Combining Higher-Order Resummation with Multiple NLO Calculations and Parton Showers in GENEVA

We extend the lowest-order matching of tree-level matrix elements with parton showers to give a complete description at the next higher perturbative accuracy in alpha_s at both small and large jet resolutions, which has not been achieved so far. This requires the combination of the higher-order resummation of large Sudakov logarithms at small values of the jet resolution variable with the full next-to-leading order (NLO) matrix-element corrections at large values. As a by-product this combination naturally leads to a smooth connection of the NLO calculations for different jet multiplicities. In this paper, we focus on the general construction of our method and discuss its application to e+e and pp collisions. We present first results of the implementation in the GENEVA Monte Carlo framework, where we employ N-jettiness as the jet resolution variable, combining its next-to-next-to-leading logarithmic resummation with fully exclusive NLO matrix elements, and PYTHIA8 as the backend for further parton showering and hadronization. For hadronic collisions, we take Drell-Yan production as an example to apply our construction. For e+e- ->jets, taking alpha_s(mZ) = 0.1135 from fits to LEP thrust data, together with the PYTHIA8 hadronization model, we obtain good agreement with LEP data for a variety of 2-jet observables.


JHEP09(2013)120
Contents Accurate and reliable theoretical predictions for measurements at collider experiments require the inclusion of QCD effects beyond the lowest perturbative accuracy in the strong coupling α s . This is especially important in the complex environment of the LHC, which requires precise predictions for a broad spectrum of observables. Higher-order corrections in α s are important to predict total cross sections and other inclusive observables. Exclusive jet observables, such as jet-vetoed cross sections, require the all-orders resummation of logarithmically enhanced contributions. For many observables, an accurate description across phase space demands a combination of both types of corrections. For experimental analyses to benefit from these advances, it is crucial to provide the best possible theoretical predictions in the context of fully exclusive Monte Carlo event generators. The goal of modern Monte Carlo programs is to provide a proper description of the physics at every jet resolution scale. This is the motivation for the by-now standard combination of matrix elements with parton showers (ME/PS). [1,2] Here, the parton shower provides the correct lowest-order description at small jet resolution scales, where the resummation of large Sudakov logarithms is needed, while at large jet resolution scales the exact tree-level matrix elements are needed to provide the correct lowest-order description. Hence, the ME/PS merging provides theoretical predictions at the formally leading O(1) accuracy relative to the lowest meaningful perturbative order. Once one has a consistent matching between these two limits of phase space, the possibility to include exact tree-level matrix elements for several jet multiplicities follows almost automatically by iteration.
Given the necessity of higher-order perturbative corrections to make accurate predictions, it is important to extend the perturbative accuracy of the Monte Carlo description to formal O(α s ) accuracy relative to the lowest order. This requires including the formally next higher-order corrections that are relevant at each scale. At small scales, i.e., small values of the jet resolution variable, this requires improving the leading-logarithm (LL) parton shower resummation with higher-order logarithmic resummation, while at large scales this requires including the fully differential next-to-leading-order (NLO) matrix elements. It is important to realize that typically a large part of phase space, often including the experimentally relevant region, is characterized by intermediate scales, i.e., by a transition from small to large scales. In the end, providing an accurate description of this transition region requires a careful combination of both types of corrections.
Such a Monte Carlo description at relative O(α s ) accuracy across phase space has never been achieved and is the subject of our paper. (We briefly summarize the existing efforts to combine NLO corrections with parton showers in section 1.1 below.) The crucial starting point in our approach is that all perturbative inputs to the Monte Carlo are formulated in terms of well-defined physical jet cross sections [3,4]. This allows us to systematically increase the perturbative accuracy by incorporating results for the relevant ingredients to the desired order in fixed-order and resummed perturbation theory.
An essential aspect of any higher-order prediction is a reliable estimate of its perturbative uncertainty due to neglected higher-order corrections. To the extent that parton shower Monte Carlos provide perturbative predictions, they should be held to the same JHEP09(2013)120 standards. An important benefit in our approach is that we have explicit control of the perturbative uncertainties and are able to estimate reliable fixed-order and resummation uncertainties. As a result, in Geneva each event comes with an estimate of its perturbative uncertainty; i.e., Geneva provides event-by-event theory uncertainties. 1 In our approach, the Monte Carlo not only benefits from the resummation, but in turn also provides important benefits to analytic resummed predictions. For one, it greatly facilitates the comparison with experimental data, as it allows easy application of arbitrary kinematic selection cuts, which can often be tedious to take into account in analytic predictions. More importantly, resummed predictions require nonperturbative inputs which can be treated as power corrections at intermediate scales but become O(1) corrections at very small scales. Here, these are provided "on-the-fly" by the nonperturbative hadronization model. In essence, we are able to combine the precision and theoretical control offered by higher-order resummed predictions with the versatility and flexibility offered by fully exclusive Monte Carlos.
In this paper, we focus on the theoretical construction. We leave a discussion of the implementation details of the Geneva Monte Carlo framework to a dedicated publication. 2 We will however highlight some of the main technical issues we had to overcome and discuss some implementation details in the application sections. In the remainder of this section, we briefly summarize the existing efforts to include NLO corrections in parton shower Monte Carlos and give a short overview of our basic construction. In section 2, we discuss in detail the requirements to obtain full α s accuracy as well as our method to achieve it. In section 3, we discuss the application to e + e − → jets, where we combine next-to-next-toleading logarithmic (NNLL) resummation with NLO matrix elements, and present results from the implementation in Geneva together with a comparison to LEP measurements. In section 4, we discuss the application to hadronic collisions and show first results for Drell-Yan production, pp → Z/γ * → + − + jets, obtained with Geneva. We conclude in section 5.
Recently, efforts have been made to extend these approaches in order to combine NLO matrix elements for several jet multiplicities with parton showers [26][27][28][29][30][31][32]. We will discuss some issues faced by some of these approaches in section 2.1.5. Here, we would like to stress that including several NLO matrix elements by itself does not provide a full extension of 1 Further uncertainties, e.g. due to nonperturpative effects such as hadronization, must be evaluated as well for a complete uncertainty analysis. 2 The current Geneva framework and implementation is new and independent of the earlier work in refs. [3,4].

JHEP09(2013)120
the lowest-order ME/PS matching to relative O(α s ) perturbative accuracy, since the fixed NLO corrections only suffice to increase the perturbative accuracy in the region of large jet resolution scales. To the same extent that the inclusion of the LL Sudakov factors in the ME/PS merging are needed to get meaningful results at intermediate and small jet scales, higher-order resummation is necessary to improve the perturbative accuracy in this region. In our approach, the full information from NLO matrix elements for several jet multiplicities is automatically included as follows: for a given Born process with N partons, a small jet scale corresponds to the exclusive N -jet region, and here the N -parton virtual NLO corrections are incorporated in conjunction with the higher-order resummation; in fact, they are a natural ingredient of it. On the other hand, a large jet scale corresponds to the inclusive (N + 1)-jet region with additional hard emissions. Here, the (N + 1)-parton virtual NLO corrections are included in the usual way by the fixed NLO calculation for N + 1 jets.

Brief overview of our construction
The starting point of our approach is the separation of the inclusive N -jet cross section into an exclusive N -jet region and an inclusive (N + 1)-jet region, Here T ≡ T (Φ N +1 ) is a suitable resolution variable, which is a function of Φ N +1 , and dσ/dΦ N +1 (T ) denotes the fully differential cross section for a given T . In ME/PS merging, this role is played by the variable that determines the merging scale. However, in our case the parameter T cut is not a jet-merging cut but instead serves as an infrared cutoff for the calculation of dσ/dΦ N +1 (T ) and ideally is taken as small as possible.
In the N -jet region at small T (both above and below T cut ), logarithms of T become large and must be resummed to maintain consistent perturbative accuracy to some order in α s . On the other hand, in the (N + 1)-jet region at large T , a fixed-order expansion in α s will suffice. To consistently match the resummed and fixed-order calculations, we use the following prescription for the jet cross sections entering in eq. (1.1): The superscript "resum" indicates an analytically resummed calculation and "FO" indicates a fixed-order calculation or expansion. This construction properly reproduces the fixed-order calculation at large T , the resummed calculation at small T , and smoothly interpolates between them. It is straightforward to extend our formulation to combine higher jet multiplicities at NLO with higher-order resummation, as we will show. This is done by replacing dσ FO /dΦ N +1 in eq. (1.2) with an inclusive (N + 1)-jet cross section separated into the exclusive (N +1)-jet and inclusive (N +2)-jet cross sections and iteratively applying eq. (1.2).

JHEP09(2013)120
The key ingredients in our approach are the higher-order resummation of the jet resolution variable, the fully differential fixed-order calculation, and the parton shower and hadronization. While each of these components is known, there is a sensitive interplay of constraints between them that must be satisfied to achieve a consistent combination. This is precisely what is accomplished in the Geneva framework and is the focus of this paper.

General construction
In this section, we derive our theoretical construction in a process-independent manner. We start in section 2.1 with a slightly simplified setup, considering the singly differential spectrum in the jet resolution variable. We use this to discuss in detail the perturbative structure and the accuracy in the different phase space regions. In section 2.2, we discuss the extension to the fully differential case and how to combine the fixed-order expansion and resummation in this situation. In section 2.3, we further generalize these results to include several jet multiplicities by iteration. Finally, in section 2.4, we discuss the Monte Carlo implementation and how to attach parton showering and hadronization.

Basic setup
The basic idea of Monte Carlo integration is to randomly generate points in phase space ("events") that are distributed according to some differential (probability) distribution. By summing over all points that satisfy certain selection criteria, we are able to perform arbitrary integrals of the distribution. In our case, that distribution is the fully differential cross section, allowing one to compute arbitrary observables. For simplicity, we will first focus on the singly differential cross section in some phase space resolution (or jet resolution) variable T of dimension one. The precise definition of T is not important at the moment, so we keep it generic for now. We use the convention that the limit T → 0 corresponds to Born kinematics, i.e., the tree-level cross section is ∼ δ(T ). We also require that T is an IR-safe observable, such that the differential cross section dσ/dT can in principle be well defined to all orders in perturbation theory and for T > 0 contains no IR divergences.
To give an example, for our application to e + e − → 2/3 jets in section 3, we will use 2-jettiness T 2 = E cm (1 − T ), where T is the usual thrust [33]. Alternatives include other 2-jet event shapes. For Drell-Yan in section 4, we will use beam thrust [34]. An alternative would be the p T of the leading jet. If the Born cross section we are interested in has N signal jets, 3 then T could be N -jettiness or the largest p T of any additional jet. The important point is that we can think of T as a resolution variable which determines the scale of additional emissions in the Φ ≥N +1 phase space, such that for T ≤ T cut there are no emissions above the scale T cut . For later convenience, we also define the dimensionless equivalent of T as

JHEP09(2013)120
Here, Q is the relevant hard-interaction scale in the Born process, e.g., Q ≡ E cm for e + e − → jets or Q ≡ m + − for Drell-Yan pp → Z/γ * → + − . In terms of τ , the limit τ 1 corresponds to the exclusive limit close to Born kinematics. For τ ∼ 1, there are additional emissions at the hard scale T ∼ Q, which means we are far away from Born kinematics and we should switch the description to consider the corresponding Born process with one additional hard jet.
To describe the differential T spectrum, we want the Monte Carlo to generate events at specific values of T , which are distributed according to the differential cross section dσ/dT . The total cross section is then simply given by summing over all events, The essential problem every Monte Carlo generator faces is that in perturbation theory the differential cross section dσ/dT contains IR divergences from real emissions for T → 0, which only cancel against the corresponding virtual IR divergences upon integration over the T → 0 region. As a result, the perturbative spectrum for T → 0 can only be defined in a distributional sense in terms of plus and delta distributions [see eq. (2.6) below]. To deal with this, we have to introduce a small cutoff T cut and define the cumulant of the spectrum as In the Monte Carlo, the total cross section is then obtained by combining the cumulant and spectrum as In practice, this is implemented by generating two distinct types of events: (i) events that have T = 0 and relative weights given by σ(T cut ), and (ii) events that have nonzero values T > T cut and relative weights given by dσ/dT . The first type of events have Born kinematics and represents the tree-level and virtual corrections together with the corresponding real emissions integrated below T cut . The second type of events contains one or more partons in the final state, since the real-emission corrections determine the shape of the spectrum for nonzero T . We now have two basic conditions: 1. From a numerical point of view, we want the value of T cut to be as small as possible, so as to describe as much differential information as possible. In practice, our ability to reliably compute the cumulant σ(T cut ) in perturbation theory sets a lower limit on the possible value of T cut few times Λ QCD .
2. Since T cut is an unphysical parameter, we want the dependence on it to drop out (to the order we are working at). In practice, this is guaranteed by including the corresponding dominant higher-order corrections in the cumulant and spectrum. JHEP09(2013)120

Perturbative expansion and order counting
In perturbation theory, the differential cross section in τ and the cumulant in τ cut have the general form where we have distinguished "singular" and "nonsingular" contributions. For τ → 0, the singular terms in dσ sing /dτ scale like 1/τ , while the nonsingular terms in dσ nons /dτ contain at most integrable singularities. For the cumulant, this means that σ sing (τ cut ) contains all terms in σ(τ cut ) enhanced by logarithms ln k (τ cut ), while σ nons (τ cut = 0) = 0. The singular part of the spectrum is given by where σ B denotes the Born cross section, and we denote the usual plus distributions as They encode the cancellation between real and virtual IR divergences. The corresponding singular contribution to the cumulant cross section integrated up to τ ≤ τ cut is At O(α k s ), only the coefficients C n (α s ) with n ≤ 2k −1 contribute, so dσ/dτ has logarithms up to α n s L 2n−1 /τ , while σ(τ cut ) has logarithms up to α n s L 2n cut , where we use the abbreviations L ≡ ln(τ ) , L cut ≡ ln(τ cut ) . (2.9) The α s expansion of the coefficients C −1 (α s ) and C n (α s ) in the singular contributions can be written as Similarly, the α s expansion of the nonsingular contributions can be written as Using eqs. (2.10) and (2.11), the spectrum and cumulant up to O(α 2 s ) are given by

JHEP09(2013)120
Note that the c k,−1 constant term in the singular corrections, which contains the finite virtual corrections to the Born process, only appears in the cumulant. We now distinguish three parametrically different regions in τ , which are illustrated in figure 1: • Resummation ("peak") region τ 1: in this limit, the logarithms in the singular contributions are large, such that parametrically one has to count 4 This means one has to resum the towers of logarithms (α s L 2 ) n in the spectrum and (α s L 2 cut ) n in the cumulant in eq. (2.12) to all orders in α s to obtain a meaningful perturbative approximation at some order. At the same time, the nonsingular corrections can be regarded as power suppressed, since they are of relative O(τ ).
• Fixed-order ("tail") region τ ∼ 1: in this limit, the logarithms are not enhanced, and a fixed-order expansion in α s is applicable. The singular and nonsingular contributions are equally important and both must be included at the same order in α s . In particular, there are typically large cancellations between these for τ ∼ 1, so it is actually crucial not to resum the singular contributions in this region, since otherwise this cancellation would be spoiled.
• Transition region: the transition between the resummation and fixed-order regions.
There are of course no strict boundaries between the different regions. This is why it is important to have a proper description not just in the two limits but also in the transition region, which connects the resummation and fixed-order regions. In fact, in practice the experimentally relevant region is often somewhere in the transition region, where both types of perturbative corrections can be important.

Lowest perturbative accuracy
For the Monte Carlo to provide a proper description at all values of T , it has to include at least the lowest-order terms relevant for each region. Keeping only these, and dropping all other terms, the spectrum and the cumulant are given by  where the functions f 0,1 and F 0,1 are given in terms of the coefficients c ij in eq. (2.10) as The f 0 and F 0 resum the leading-logarithmic series in the cross section, which we denote as LL σ . The functions f 1 and F 1 resum the next-to-leading-logarithmic series in the cross section, which we denote as NLL σ . In the resummation region at τ 1, the LL σ terms in the spectrum scale as L ∼ 1/ √ α s (relative to the overall α s /τ scaling) and provide the lowest level of approximation. The NLL σ terms scale as ∼ 1, and one can argue about whether they are needed as well in order to get a meaningful lowest-order prediction. Formally, they are necessary to obtain the spectrum at ∼ α s /τ , which one might consider the natural leading-order scaling of the spectrum (or equivalently if one does not want to rely on the ∼ 1/ √ α s enhancement of the LL series). Experience shows that the NLL terms are indeed numerically important. For example, in analytic resummations, one rarely gets a sensible prediction without going at least to NLL. Similarly, to obtain sensible predictions from a parton shower, it is almost mandatory to include important physical effects such as momentum conservation in the parton splitting and the choice of α s scale [35]. In the cumulant, the LL σ series in F 0 scales as ∼ 1 and must be included. The NLL σ series in F 1 scales as ∼ √ α s and, for consistency, should be included in the cumulant if it is included in the spectrum.
In the fixed-order region at τ ∼ 1, the lowest meaningful order in the spectrum is given by the complete O(α s ) terms, requiring one to include the c 11 and c 10 terms, which are part of the f 0 and f 1 functions, as well as the nonsingular corrections f nons 1 (τ ). Since we take τ cut to be small, the cumulant is always in the resummation region. Hence, its nonsingular corrections F nons 1 (τ ) [see eq. (2.12)] are suppressed by O(α s τ cut ) and can be safely neglected.
The leading level of accuracy in eq. (2.14) closely corresponds to what is achieved in the standard ME/PS matching. In this case, the LL resummation is provided by the JHEP09(2013)120 parton shower Sudakov factors (either generated by the shower or multiplied by hand), where the jet resolution variable corresponds to the shower evolution variable, since that is the variable for which the shower directly resums the correct LL σ series. The LL σ series has a well-known and simple exponential structure, The resummation exponent at LL σ is given by the integral over the leading c 11 α s ln(τ )/τ term in the spectrum. This is precisely what the standard parton shower veto algorithm exploits to generate the resummation exponent. The analogous structure does not hold at NLL σ , which is why the parton shower cannot resum the NLL σ series by exponentiating the integral of the c 10 α s /τ term. As already mentioned, in practice, parton showers include important partial NLL effects, so practically this provides a numerically close approximation to the correct NLL σ series. The nonsingular corrections in the spectrum, f nons 1 (τ ), are obtained by including the full tree-level matrix element for one additional emission. Since the full matrix element also includes the c 11 and c 10 terms, this requires a proper matching procedure to avoid double counting these terms. At LL σ , a simple way to do this is to multiply the full fixed-order result from the matrix element with the shower's LL σ resummation exponent, which corresponds to the CKKW-L [1,2,36,37] procedure. The reason this gives the spectrum correctly at LL σ is the simple structure in eq. (2.17), where the LL σ exponent multiplies the c 11 term in the spectrum. 5 At large τ ∼ 1, the exponent in eq. (2.18) can be expanded as 1 + O(α s ), so eq. (2.18) gives the correct leading fixed-order result. Compared to eq. (2.14), the NLO matching performed in MC@NLO and Powheg amounts to adding to the cumulant the c 1,−1 singular constant, containing the O(α s ) virtual corrections, as well as the nonsingular contributions F nons 1 (τ ). Assuming the same set of NLL terms are included in the cumulant and spectrum, this achieves that inclusive quantities that are integrated over a large range of τ , such as the total cross section, are correctly reproduced at fixed NLO, which provides them with O(α s ) accuracy. In these approaches, the goal is not to improve the perturbative accuracy of the spectrum (or the cumulant at small τ cut ), which has the same leading accuracy as in eq. (2.18).

Next-to-lowest perturbative accuracy
We now want to improve the Monte Carlo description in eq. (2.4) from the lowest-order accuracy, given by eq. (2.14), to the next-to-lowest perturbative accuracy in α s . This requires us to include the appropriate higher-order corrections in each region, which gives where we denote the series of logarithms resummed by the functions f 2 and F 2 by NLL σ and the series resummed by f 3 and F 3 by NNLL σ . They can again be written in terms of the c ij coefficients in eq. (2.10) as In the resummation region, the NLL σ series in the spectrum scales as ∼ α 3/2 s and thus provides the ∼ α s correction to the LL σ series in f 0 . Similarly, the NNLL σ series scales as ∼ α 2 s providing the ∼ α s correction to the NLL σ series in f 1 . They can again be obtained by performing the standard resummation in the exponent of the cross section to NLL and NNLL respectively. (Here, NLL refers to those parts of the full NNLL resummation that arise from the combination of the one-loop matching corrections with the NLL resummation, see section 3.1.1 and table 2.) In the fixed-order region, increasing the perturbative accuracy by ∼ α s requires the complete O(α 2 s ) corrections, including the f nons 2 (τ ) nonsingular corrections. Similarly, for the cumulant, F 2 and F 3 resum the NLL σ and NNLL σ series of logarithms, which scale as ∼ α s and ∼ α 3/2 s , respectively, and provide the ∼ α s improvement over the LL σ and NLL σ series in F 0 and F 1 . In addition, going to the next higher order in the cumulant requires including the full singular constant c 1,−1 , 6 as well as the nonsingular corrections F nons 1 (τ ), which both scale as ∼ α s .
It is instructive to see where the information from the virtual NLO matrix elements enters in eq. (2.19). As already mentioned, the virtual NLO corrections to the Born process are given by c 1,−1 . In addition, by multiplying the LL series it contributes part of f 2 and F 2 . Hence, consistently combining the virtual corrections with the resummation requires one to go to at least NLL . The virtual NLO corrections with one extra emission (plus the JHEP09(2013)120 integral over the two-emission tree-level matrix element) yield the full O(α 2 s ) corrections in the spectrum, i.e., both the singular c 2k terms as well as the nonsingular f nons 2 terms in eq. (2.12). Adding these corrections again requires one to avoid double counting the singular c 2k terms that are already included in the resummation. In analytic resummation, it is well known how to do this, namely by simply adding the nonsingular corrections. These are obtained by taking the difference of the full NLO corrections and the singular NLO corrections, where the latter are given by expanding the resummed result to fixed order. Since this construction involves the virtual contribution to both the Born process and the process with one extra emission, we see that going consistently to higher order in both the resummation and fixed-order regions naturally leads to a combination of the information from two successive NLO matrix elements.

Merging NLO matrix elements with parton shower resummation only
We stress that, for a description at the next-higher perturbative accuracy across the whole range in τ , it is not sufficient to include the fixed NLO corrections to the spectrum and take care of the double counting with the parton shower resummation. This only provides the proper NLO description in the fixed-order region at large τ . In the transition and resummation regions, a proper higher-order description necessitates higher-order resummation. Of course, this is not a problem if the only goal is to improve the fixed-order region at large τ , as is the case for example in a recent MC@NLO publication [30].
However, including the fixed NLO corrections outside the fixed-order region, as is done in Sherpa's recent NLO merging [28,29], can actually make things worse in two respects: first, numerically this will typically force the spectrum to shift toward the fixedorder result and away from the resummed one. Since this can shift the spectrum in the wrong direction, it can potentially make the result less accurate. 7 At the same time, the perturbative uncertainties from fixed-order scale variation decrease, which only aggravates this problem. Multiplying the NLO corrections to the spectrum with LL parton shower Sudakov factors (see, e.g., ref. [38]) can mitigate this to some extent but does not solve the problem. The only consistent way to include the fixed NLO corrections to the spectrum outside the fixed-order region, and in particular obtain reliable perturbative uncertainties, is to properly combine them with a higher-order resummation.
Second, this explicitly spoils the formal O(α s ) accuracy of the inclusive cross section. To see this, consider adding the fixed NLO corrections to the lowest-order spectrum and cumulant in eq. (2.14), properly taking care of the double counting at O(α 2 s ), which gives One can see this for example in the case of 2-jettiness in figure 3 in section 3. Here, the NLL +LO3 result is much closer to the slightly higher NNLL +NLO3 best prediction than the fixed NLO3 result. We have checked that in this case, adding the NLO3 to the NLL +LO3 by expanding it to O(α 2 s ) forces the result to move in the wrong direction toward the lower NLO3.

JHEP09(2013)120
Using these expressions yields for the inclusive cross section While the first two terms give the correct NLO inclusive cross section, the O(α 2 s ) terms induced by the fixed NLO corrections in the spectrum formally scale as α s and α 3/2 s and therefore spoil the formal O(α s ) perturbative accuracy for the inclusive cross section and in fact for any inclusive observable. This directly contradicts the claim in refs. [28,29] that this description maintains the higher-order accuracy of the underlying matrix elements in their respective phase space range. It only preserves the fixed O(α s ) terms, which in the context of combining fixed-order corrections with a logarithmic resummation is necessary but not sufficient to preserve the higher perturbative accuracy.
This problem cannot be avoided by multiplying the α 2 s corrections in the spectrum with the LL parton shower Sudakov factors, since this does not provide the proper NLL σ and NNLL σ series. Note also that we have already assumed in eq. (2.22) that the full NLL σ series is included in the spectrum and cumulant. In general, the parton shower cannot provide this, which means there will be even α 2 s L 3 cut ∼ √ α s terms induced in eq. (2.22).
Pragmatically, the inclusive cross section can be restored to formal O(α s ) accuracy by either explicitly including the corresponding α 2 s corrections in the cumulant to cancel these terms, where numerical methods to do so have been described very recently in refs. [31,32,39], or alternatively by explicitly restricting the fixed NLO corrections in the spectrum to the fixed-order region at large τ , such that the induced O(α 2 s ) terms in the total cross section are not logarithmically enhanced and are formally O(α 2 s ). This is essentially the approach taken in ref. [30]. However, neither of these approaches improves the perturbative accuracy in the spectrum outside the fixed-order region.

What Monte Carlo can do for resummation
For T being the resolution variable between N and more than N jets, we showed in the previous subsection that combining the NLO matrix-element corrections for N and N + 1 partons at the level of the singly differential T spectrum is equivalent to combining the NNLL resummation of the singular contributions with the higher-order nonsingular contributions. Our goal now is to extend this singly differential description to the fully differential case, in order to use the full N -parton and (N + 1)-parton information of the matrix elements. We will use the notation (N)LO N or (N)LO N +1 to indicate up to which fixed order in α s the N -parton or (N + 1)-parton matrix elements are included.
To start with, it is straightforward to generalize the jet resolution spectrum dσ/dT and its cumulant σ(T cut ) to include the full dependence on the N -body Born phase space, Here, dσ incl /dΦ N is the inclusive N -jet cross section. The discussion in section 2.1 can be precisely repeated in this case, since the perturbative structure of the differential spectrum dσ/dΦ N dT with respect to T is precisely the same as in eqs. (2.5) and (2.6). Namely, we can write it as the sum of singular and nonsingular contributions, The nonsingular contributions are general functions of Φ N and T , but as before are integrable in T for T → 0. The singular contributions have the structure where dσ B /dΦ N is now the fully differential Born cross section. Since the singular contributions arise from the cancellation of virtual and real IR singularities, which only know about Φ N , their T dependence naturally factorizes from the Φ N kinematics of the underlying hard process. This is what allows the resummation of the singular terms to higher orders for a given point in Φ N . At LL, the entire Φ N dependence is that of the Born cross section. At higher logarithmic orders, this is not the case anymore, since the coefficients C n can have nontrivial Φ N dependence. In addition, the precise definition of T also becomes important. Depending on its definition, the higher-order singular coefficients can depend on clustering effects or other types of nonglobal logarithms [40][41][42][43][44][45], which can be difficult to resum to high enough order with currently available methods. Therefore, it is important to choose a resolution variable with simple resummation properties. An example is N -jettiness, for which the complete NNLL resummation for arbitrary N is known [46,47]. For the purpose of our discussion below, we will assume that a resummed result for the spectrum and its cumulant in eq. (2.24) at sufficiently high order is available to us. We can think of the cumulant dσ/dΦ N (T cut ) in eq. (2.24) as the exclusive N -jet cross section with no additional emissions (jets) above the scale T cut , while the spectrum dσ/dΦ N dT for T > T cut is the corresponding inclusive (N + 1)-jet cross section. While the cumulant dσ/dΦ N (T cut ) is differential in dΦ N and thus already as differential as it can be, the spectrum contains a projection from the full dΦ ≥N +1 phase space down to dΦ N dT . To also be fully differential in the (N + 1)-jet phase space, we can generalize eq. (2.24) to where dσ/dΦ N +1 (T ) denotes the fully differential spectrum for a given T ≡ T (Φ N +1 ). We explicitly denote the dependence on T and T cut to clearly distinguish the spectrum from the cumulant. We have also used the shorthand notation

JHEP09(2013)120
where Φ N (Φ N +1 ) denotes a projection from an (N + 1)-body phase space point to an Nbody phase space point. This projection defines what we mean by N jets at higher orders in perturbation theory. Note that beyond LO, both the cumulant dσ/dΦ N (T cut ) and spectrum dσ/dΦ N +1 (T ) must be well-defined jet cross sections; i.e., they require a specific IR-safe projection from Φ ≥k+1 to Φ k for both k = N and k = N +1. We will see below where this definition enters. Using eq. (2.27) at the next-higher perturbative accuracy requires us to combine the higher-order resummation in T for the cumulant and spectrum with the fully exclusive N -jet and (N + 1)-jet fixed-order calculations at NLO N and NLO N +1 . To achieve this, we have to construct appropriate expressions for the cumulant dσ/dΦ N (T cut ) and the spectrum dσ/dΦ N +1 (T ), which we do in the next two subsections.

Matched cumulant
We start by discussing the cumulant in eq. (2.27). Since the resummation is naturally differential in the dΦ N of the underlying Born process, we can combine the resummed result with the fixed-order one by adding the fixed-order nonsingular contributions to it, (2.29) The first term contains the resummed contributions, while the difference of the two terms in square brackets provides the remaining nonsingular corrections that have not already been included in the resummation. The NLO N fixed-order result is given by where B N and B N +1 are the N -parton and (N + 1)-parton tree-level (Born) contributions, V N is the N -parton one-loop virtual correction, and we abbreviated Here, T (Φ N +1 ) implements the definition of T . The NLO N result also depends on the projection from Φ N +1 to Φ N , i.e., the precise NLO definition of Φ N . However, this dependence only appears in the nonsingular corrections. For a given definition of T , the singular NLO corrections do not depend on how the remaining Φ N +1 phase space is projected onto Φ N , since they arise from the IR limit in which all (IR-safe) definitions agree. In eq. (2.29), the singular contributions inside the full fixed-order cumulant, dσ FO /dΦ N (T cut ) are canceled by the NLO expansion of the resummed result at NLL σ or higher, leaving only the nonsingular fixed-order contributions in square brackets.

Matched spectrum
To properly combine the higher-order resummation in T with the fully differential (N + 1)jet fixed-order calculation, the inclusive (N + 1)-jet spectrum dσ/dΦ N +1 (T ) in eq. (2.27)

JHEP09(2013)120
has to fulfill two basic matching conditions, The first condition states that integrating the fully differential spectrum over the additional radiative phase space has to reproduce the correct spectrum in T including the desired resummation and fixed-order nonsingular corrections, such that eq. (2.27) reproduces eq. (2.24). The second condition states that the fixed-order expansion of the fully differential spectrum has to reproduce the full (N + 1)-jet fixed-order calculation, where at NLO N +1 , Here, B N +1 and B N +2 are the (N + 1)-parton and (N + 2)-parton tree-level (Born) contributions, and V N +1 is the (N + 1)-parton one-loop virtual correction. Integrating over dΦ N +2 in the last term now requires a projection from Φ N +2 to Φ N +1 , analogous to eq. (2.28), which now defines precisely what we mean by N + 1 jets at NLO. In principle, there is some freedom to construct an expression for dσ/dΦ N +1 (T ) that satisfies both conditions to the order one is working. Our master formula to combine the resummed spectrum dσ resum /dΦ N dT with the fully differential dσ FO /dΦ N +1 is given by Expanding the right-hand side to a given fixed order, we can see immediately that Condition 2 is satisfied by construction. Imposing Condition 1 yields the consistency (or "matching") condition If the resummed result already has the nonsingular contributions at the desired fixed order added in, then the term in brackets is by construction equal to unity for any value of T . Otherwise, the expansion of the resummed result reproduces the singular terms of the full fixed-order result, leaving the nonsingular fixed-order contributions, such that we get Here, dσ sing,resum denotes the pure resummed result only containing the resummation of the singular contributions. Hence, eq. (2.36) not only multiplies in the additional dependence JHEP09(2013)120 on Φ N +1 /Φ N at fixed order, but if needed also adds the nonsingular corrections to the spectrum multiplied by the higher-order resummation factor. (Note that for the expansion of the resummed result to indeed reproduce all the singular terms at the desired fixed order, the resummation has to be carried out to sufficiently high order, which we have already seen in section 2.1.) To apply Condition 1, we have to integrate eq. (2.34) using the projection onto Φ N and T in eq. (2.31). Therefore, to get the correct T spectrum at NLO N +1 , the projection in eq. (2.35) has to satisfy i.e., it has to preserve the value of T when constructing the projected Φ N +1 point. Usually, the simplest way to handle this would be to use the left-hand side to define T (Φ N +2 ). However, in our case, eq. (2.39) provides a very nontrivial condition on the projection since T (Φ N +2 ) is already defined by our choice of jet resolution variable, which in particular has to be resummable. This turns out to be a nontrivial technical challenge one has to overcome to be able to satisfy Condition 1. We will see where this enters in section 3.1.2 and section 4.1.2.
Note that to ensure that the resummation factor in square brackets in eq. (2.36) is well behaved in the fixed-order region at large T , it is important to turn off the resummation such that the ratio of the resummed spectrum and its expansion becomes O(1) up to higher fixed-order corrections. In principle, the fixed-order result in the denominator can also become negative at very small values of T . This is not a problem in practice, since this region is explicitly avoided by imposing the cut T > T cut .

Perturbative accuracy and order counting
The appropriate order counting in the resummation and fixed-order regions is precisely the same as in section 2.1, so there is no need to repeat it here. Applying eq. (2.36) at the very lowest order, namely LL σ resummation with LO N +1 fixed-order corrections, we get at small T , and we used that at LL σ the ratio in brackets in eq. (2.36) is just the resummation exponent. This directly corresponds to the CKKW-L procedure [1,2,36,37], which multiplies the tree-level matrix elements with the shower Sudakov factors. Hence, we can think of our master formula eq. (2.36) as a consistent extension of this to higher orders. As demonstrated in section 2.1, going to the next higher perturbative accuracy in all phase space regions requires the NLL σ and NNLL σ series of logarithms. We obtain these by performing the full NLL and NNLL resummation in the exponent, as well as the fixed NLO N and NLO N +1 corrections in the cumulant and spectrum, respectively. The resummation naturally connects both jet multiplicities, since the NLO N corrections are included in the cumulant and are part of the resummation for the spectrum starting JHEP09(2013)120 inclusive N -jet exclusive N -jet inclusive (N +1)-jet notation fixed order accuracy log. order accuracy fixed order accuracy at NLL , where they effectively predict the singular NLO N +1 contributions, and the full NLO N +1 corrections are obtained by adding the nonsingular corrections to the spectrum. In the following, we will use the notation (N)NLL T +(N)LO N +1 to indicate the resummation order for the employed jet resolution variable together with the (N + 1)-jet fixed order. For simplicity, we do not explicitly denote the N -jet fixed order and keep it implicit in the resummation order, i.e., LO N at (N)LL and NLO N at NLL and above. This is summarized in table 1.
An immediate and important question to ask is to what accuracy resummed spectra for jet resolution variables other than T are predicted in our approach. A detailed theoretical investigation of the formal resummation order one attains for other variables would be very interesting but is beyond the scope of the present work. What is certainly clear is that other variables will not be resummed at the same formal level as the primary jet resolution variable T itself. However, we know that other variables are correct to NLO N +1 , while at the same time, the inclusive cross section is not changed, as it is independent of which variable one integrates over. This implies that the NLO N +1 corrections for other variables do not induce uncanceled, higher-order logarithmic terms as in eq. (2.22), and hence, some higher-order resummation must be partially retained for other observables as well. Numerically, the higher-order resummation in T provides an improved weighting of the IR region of phase space, from which other variables are expected to benefit as well. We can validate to what accuracy other variables are obtained by comparing predictions from our highest order to the analytically resummed results for other observables, which we do in section 3.3.

Extension to more jet multiplicities
The method proposed in this paper is completely general and can be extended to more jet multiplicities essentially by iterating the procedure discussed in section 2.2. We start by introducing separate jet resolution variables T N to distinguish N from N + 1 jets, T N +1 to distinguish N + 1 from N + 2 jets, and so on. One can choose any IR-safe observable that goes to zero in the limit of N pencil-like jets. For each N , the inclusive N -jet cross section JHEP09(2013)120 is obtained by combining the cumulant and spectrum for T N as in eq. (2.27), The exception is the highest jet multiplicity, N max , for which T cut Nmax = ∞, corresponding to the fact that no additional jets are resolved.
For the cumulants in eq. (2.41), the discussion in section 2.2.1 applies separately for each N , so the cumulants matched to higher resummed and fixed order are given, as in eq. (2.29), by (2.42) The fully differential T N spectra dσ/dΦ N +1 (T N ) are now obtained recursively as follows. We start with the highest jet multiplicity, N max , for which no resummation is needed since T cut Nmax is essentially removed. Furthermore, the highest jet multiplicity is, by construction, only required at leading order, where the result is simply given by the Born contribution, For each N < N max , we apply the discussion in section 2.2.2. To combine the resummation in T N with the (N + 1)-jet fixed-order calculation, the fully differential T N spectrum dσ/dΦ N +1 (T N ) must satisfy the matching conditions as in eqs. (2.32) and (2.33), These can be satisfied by a straightforward generalization of eq. (2.36), The prefactor on the right-hand side is now the inclusive (N + 1)-jet cross section from eq. (2.41). This is what ties together the different jet multiplicities. The condition in eq. (2.45) now leads to the consistency condition

JHEP09(2013)120
which states that for each N , the cumulant and spectrum in T N +1 must be included to sufficiently high order so as to reproduce the (N + 1)-jet fixed order that is required by the T N spectrum. Imposing the condition in eq. (2.44) yields the consistency condition for the T N spectrum, which is the generalization of eq. (2.37). To satisfy eq. (2.47) at NLO N +1 , it requires that as in eq. (2.39). That is, for each N , the projection from Φ N +2 to Φ N +1 which defines the (N + 1)-jet cross section at NLO has to preserve the value of T N . In addition, eq. (2.48) requires that, upon integration, the T N +1 resummation contained in dσ incl /dΦ N +1 does not interfere with the T N resummation, e.g., by inducing higher-order logarithms in T N . Since eq. (2.41) relates the dσ incl /dΦ N to dσ/dΦ N +1 (T N ), the relationship in eq. (2.46) gives rise to a recursive definition, which when combined with the result for the highest jet multiplicity in eq. (2.43) determines dσ/dΦ N +1 (T N ) for all N .
In the Monte Carlo implementation, the phase space is split up recursively as . . . , where in each step, the total cross section for N or more jets is separated into an exclusive N -jet cross section, which is assigned to partonic events with N final-state partons, and the integral over the remaining cross section for N +1 or more jets. For the highest multiplicity, N max , the remaining cross section for N max or more jets is represented by events with N max final-state partons.
Note that the structure of eq. (2.50) is very similar to eq. (2.41). The crucial difference is that in eq. (2.50), each inclusive cross section on the left-hand side is the same that appears under the integral on the right-hand side in the line above. By comparing eq. (2.50) with eq. (2.41) and repeatedly inserting eq. (2.46), we obtain the higher-order, "fully resummed," exclusive N -jet cross sections that serve as inputs to the Monte Carlo. Abbreviating the resummation factor in eq. (2.46) as The careful reader will have noticed that the above is in one-to-one correspondence to the structure generated by a parton shower with up to N max emissions. The crucial difference is that, in our case, all ingredients are well-defined physical jet cross sections defined in terms of a global jet resolution variable. This allows us to systematically increase the perturbative accuracy by computing the relevant ingredients to higher order in resummed and fixed-order perturbation theory as well as to systematically estimate the perturbative uncertainties. The analogous parton-shower-like structure underlies the CKKW-L ME/PS merging, which replaces the splitting functions in the shower with the full tree-level matrix elements. Restricting eq. (2.52) to the lowest order as in eq. (2.40), it reduces to the ME/PS merging as a special case. In principle, the above construction allows us to go to even higher fixed and resummation order, as long as the fixed-order ingredients are available and the resummation is known to a correspondingly high enough order. It also lets us combine as many jet multiplicities as we like at the order they are available. In particular, it is straightforward to add additional multiplicities at the lowest accuracy in a CKKW-L-like fashion.

Attaching parton showering and hadronization
In the Monte Carlo, a point in Φ N is represented by N (massless) four-vectors together with the appropriate flavor information. We then generate events with N to N max partons and assign the N -parton events the weight dσ mc N /dΦ N (T cut N ), the (N + 1)-parton events the weight dσ mc N +1 /dΦ N +1 (T cut N +1 ), and so on. The events with N max partons are assigned the weight dσ mc (2.50) are included in the weight, which means that all events with ≥ N + 1 partons that have T N < T cut N get zero weight. 8 In this way, by summing up the weights of all events, we can integrate up the cross sections in eq. (2.50), including arbitrary kinematic cuts in Φ N , Φ N +1 , etc. What is important is that, although the events contain massless partons, they represent the exclusive jet cross sections of eq. (2.52). (From the resummation point of view, the massless partons represent the kinematics of the hard function.) In the next step, the events are given as a starting point to a parton shower, whose purpose it is to fill up the jets with additional emissions inside the jets without changing 8 Technically, the split up of phase space is usually flavor-aware. This means that an event with TN < T cut N is only set to zero if the closest two partons produce a QCD singularity.

JHEP09(2013)120
the weight of the event. Formally, this means that the shower should not be allowed to change the underlying distribution in the jet resolution variable, since this has already been computed at the higher perturbative accuracy. For example, starting from an event with N + 1 partons with kinematics Φ N +1 and weight dσ mc N +1 /dΦ N +1 (T cut N +1 ), the fully showered event should have the same jet kinematics Φ N +1 as the unshowered event from which it originated. Most importantly, the showered event should have the same value of T N (Φ N +1 ) and should have T N +1 < T cut N +1 so it still has the correct weight dσ mc N +1 /dΦ N +1 (T cut N +1 ). In the cumulant N -jet bin, the shower is allowed to fill out the phase space from T N = 0 to T cut N . Since for the highest jet multiplicity, T cut Nmax → ∞, the shower fills out the remaining phase space. In practice, these are quite nontrivial constraints on the shower. The easiest way to enforce them is to repeatedly run the shower on the same event until it produces an acceptable showered event, where we allow the value of T N to be changed at most by a numerically small amount consistent with a power correction. This method is of course computationally intensive (though it is not computationally prohibitive), since one may have to rerun the shower many times, and it would be interesting to develop a more efficient way of constraining the shower for this purpose. Notice that in this procedure no events are discarded, so the cross section is not changed.
In the final step, the showered event is passed to the hadronization routine. In this case, there are no constraints on the kinematics of the hadronized event; i.e., the hadronization is allowed to smear out the T N spectrum. The reason is that our perturbative calculation does not take into account nonperturbative effects, which are instead supplied by the hadronization. This is discussed in more detail in section 3.1.3.

Application to e + e − collisions
In this section, we apply the framework described in section 2 to e + e − → 2/3 jets, implemented in the Geneva Monte Carlo. The higher-order resummation for 2-jet event shapes in e + e − collisions is very well understood and many precise measurements from LEP exist, which are used, for example, for precise determinations of the strong coupling constant α s [48][49][50][51][52][53][54].
In this context, one important aspect is the interplay between both resummed and fixed-order perturbative contributions with the nonperturbative corrections. Here, the Geneva framework provides an important development by being able to combine the perturbative higher-order resummation with the nonperturbative information provided by Pythia's hadronization model [55,56]. For example, this allows us to use a common theoretical framework to make predictions for different phase space regions and different observables.
The e + e − implementation also provides an important and powerful validation of our approach and its practical feasibility, while avoiding the additional complications arising for hadronic collisions, such as initial-state radiation and parton distribution functions (PDFs). The implementation and first results for pp collisions are presented in section 4.
In our e + e − implementation, we use 2-jettiness, T 2 , as the 2-jet resolution variable,

JHEP09(2013)120
which is defined as [46] and is simply related to thrust T [33] by T 2 = E cm (1 − T ). Its kinematic limits are 0 ≤ T 2 ≤ E cm /2. In the limit T 2 → 0, there are precisely 2 pencil-like jets in the final state, while for T 2 ∼ E cm , there are 3 or more jets. We perform the resummation in T 2 to NNLL and include the full NLO 2 , NLO 3 , and LO 4 fixed-order matrix elements, i.e., we obtain NNLL T +NLO 3 predictions. The default running parameters for our e + e − studies are E cm = 91.2 GeV, α s (m Z ) = 0.1135, and Pythia 8.170 with e + e − tune 1. 9 Using this value of α s (m Z ) is motivated by the fact that it was obtained from fits to the thrust spectrum using N 3 LL resummation. These fits were performed in a region (corresponding to 6 GeV ≤ T 2 ≤ 30 GeV for our E cm ) where the nonperturbative corrections due to hadronization are power suppressed and can be described by a single nonperturbative parameter, which leads to a shift in the spectrum and is included in the fit in ref. [52]. We find that this value of α s (m Z ), in conjunction with Pythia's tune 1, provides overall the best description of the data, including the peak region below T 2 ≤ 6 GeV and other 2-jet event shapes. For comparison, we show results using the world average α s (m Z ) = 0.1184 [57] as well as from using Pythia tune 3.
In the next subsection, we summarize the various ingredients that go into the master formula, with the intention of giving a concise and informative overview, while leaving a detailed discussion of our implementation to a separate publication. In section 3.2, we discuss the T 2 spectrum, validating our implementation using analytic predictions as well as comparing our results to LEP data. In section 3.3, we present our results for other 2-jet variables, namely C-parameter, heavy jet mass, and jet broadening, comparing Geneva's predictions at NNLL T +NLO 3 to the analytic higher-order resummation for each variable as well as to the experimental measurements. In all cases, we find good consistency and agreement with the data.

Ingredients
The master formula is given by

JHEP09(2013)120
Its three key ingredients are the higher-order resummation of 2-jettiness, which we include at NNLL T +LO 3 , the full fixed-order matrix elements at NLO 2 , NLO 3 , and LO 4 , and the interface to parton showering and hadronization, for which we use Pythia 8. Following the construction in section 2.3 with N max = 4, the inclusive 3-jet cross section is separated into 3 and 4 or more jet contributions using 3-jettiness, T 3 , as our 3-jet resolution variable, For e + e − collisions, N -jettiness is defined by [46] where i = 1, · · · , N andn i is a unit vector along the direction of the ith jet, where the jet directions can be determined by a jet algorithm or by directly minimizing T N . 10 There are N pencil-like jets in the limit T N → 0 and N or more jets in the limit T N ∼ E cm . As discussed in section 2.3, the master formula naturally incorporates the resummation of the 3-jet resolution variable in eq. (3.4) and extends to higher jet multiplicities, i.e., N max > 4. However, since our current focus is on the main conceptual development of combining the higher-order resummation with the fixed NLO matrix elements for 2 and 3 jets, we leave these extensions to future work. As we will not be interested in the T 3 spectrum or other exclusive 3-jet observables, it is sufficient for our purposes to calculate the two terms on the right-hand side of eq. (3.4) at fixed order (i.e., we do not include resummation for T 3 ). Thus, we use In the results that follow, we use T cut 2 value between 0.5−1 GeV, which is selected randomly from a flat distribution. This smoothing out of T cut 2 avoids small numerical discontinuities that can arise with a sharp cutoff. For T cut 3 , we use T cut 3 = 2 GeV. This value is chosen small enough that the NLO 3 calculation is fully exclusive and our results are insensitive to scales below T cut 3 . Changing T cut 3 by a factor of two up and down, the results remain unchanged, with any variations well within our perturbative uncertainties.

Resummation
Our jet resolution variable, T 2 , has the important property that it can be factorized. The factorization theorem for the T 2 spectrum provides the resummed prediction that is one of the primary inputs to our master formula in eq. (3.3). It is obtained by using the framework of Soft Collinear Effective Theory (SCET) [60][61][62][63] and allows the resummation to be JHEP09(2013)120

Fixed-order corrections
Resummation input Table 2. Perturbative inputs included at a given order in resummed and fixed-order perturbation theory. The columns in the resummation input refer to the noncusp anomalous dimension (γ x ), the cusp anomalous dimension (Γ cusp ), and the QCD beta function (β).
systematically carried out to higher orders and combined with the nonsingular fixed-order result. Our highest-order resummed input to the master formula has NNLL resummation. We use the standard resummation formalism, where the large logarithms are resummed in the exponent of the cross section, with the corresponding resummation orders summarized in table 2.
We write the jet resolution distribution in T 2 as where the separation into singular and nonsingular contributions was discussed in section 2.2 [see eq. (2.25)]. The singular contribution is given by [64,65] dσ sing Here, dΦ 2 = dΩ 2 = d cos θdφ is the angular phase space for the orientation of the thrust axis with respect to the beam, and dσ B /dΩ 2 is the tree-level 2-parton cross section. Note that the overall dependence on Ω 2 here is that of the Born cross section, which is correct in the limit T 2 → 0 in which eq. (3.8) is obtained. The hard function H 2 in eq. (3.8) contains the fixed-order 2-parton matrix elements, which describe the short-distance corrections at the scale E cm . The jet functions J 1 and J 2 describe the back-to-back collinear final-state radiation along the thrust axis, and the soft function S 2 describes the soft radiation between the jets. The soft function contains perturbative and nonperturbative components, which can be separated as [66][67][68] where S pert 2 (T 2 − k, µ) is the perturbative soft function, while the shape function f (k, µ) describes the nonperturbative hadronization corrections. For T 2 ∼ Λ QCD , the shape function JHEP09(2013)120 gives an O(1) contribution to the cross section, while for T 2 Λ QCD it can be expanded, and only the leading O(Λ QCD /T 2 ) nonperturbative power correction is relevant. For further discussion and the derivation of the factorization theorem, see refs. [64,65]. The resummed prediction used in Geneva only includes the perturbative soft function, while the nonperturbative corrections are provided by the hadronization in Pythia.
The nonsingular contribution in eq. (3.7) is given by the spectrum at fixed order with the singular terms subtracted. It includes all O(T 2 /E cm ) corrections to the singular distribution to a given order in α s . The O(α s ) nonsingular corrections in T 2 are known analytically and can be taken from ref. [52], so we include them in our resummed result. Each function in eq. (3.8) depends on the renormalization scale µ and the characteristic scale of the physics it describes. These are µ H ∼ E cm , µ J ∼ √ T 2 E cm , and µ S ∼ T 2 for the hard, jet, and soft functions, respectively. Renormalization group evolution (RGE) between the soft, collinear, and hard scales resums the logarithms of the form ln µ S /µ H ∼ ln T 2 /E cm and ln µ 2 J /µ 2 H ∼ ln T 2 /E cm in the factorized singular distribution in eq. (3.8). The anomalous dimensions and singular fixed-order corrections required at a given resummation order are summarized in table 2.
The resummed cumulant in eq. (3.3) is obtained in an analogous way to the resummed T 2 distribution. It is given by a singular and nonsingular component, where the singular contribution is obtained by integrating eq. (3.8) over T 2 from 0 to T cut 2 . The nonsingular contribution to the cumulant is given by the difference between the fixed-order result and the resummed singular terms expanded to fixed order.
The perturbative uncertainties in the resummed spectrum are estimated by scale variation and receive a contribution from two distinct sources, the fixed-order corrections and the higher-order logarithmic resummation. The fixed-order uncertainties are estimated by a correlated overall variation of all scales by factors of two. The resummation uncertainties are instead estimated by varying the lower scales µ J (T 2 ) and µ S (T 2 ), which are functions of T 2 , and are referred to as profile scales [52,68,69]. The profile scales satisfy the criteria that, in the resummation region, µ J,S (T 2 ) have their canonical scaling (given above) and in the fixed-order region, µ J,S (T 2 ) ∼ µ H , which turns off the resummation. In the transition region, the profile scales provide a smooth interpolation between the resummation and fixed-order regions. These three regions are determined based on where the fixed-order singular contributions dominate over the nonsingular ones. The variations in the profile scales subject to the above constraints determine the resummation uncertainty, where we take the largest absolute variation from the central scale. The resummation uncertainties are combined in quadrature with the fixed-order uncertainties to generate our theory uncertainty estimate. For a given partonic event in Geneva, each profile scale variation gives rise to a different event weight, which is computed analytically. Hence, we can provide each event with its own perturbative uncertainty estimate by assigning it several weights from the profile scale variation in addition to its central weight.

Fixed order
As we can see from eqs. (3.4) and (3.6), we need the 3-jet cumulant dσ/dΦ 3 (T cut 3 ) as well as the Born 4-parton cross section B 4 (Φ 4 ). The Born 4-parton cross section is trivial and requires no further discussion. To calculate the 3-jet cumulant at NLO 3 , we use the generic formula given in eq. (2.34), (3.12) The projection Φ T 3 (Φ 4 ) defines what we mean by Φ 3 at NLO 3 . It implicitly depends on our choice of resolution variable since eq. (2.39) requires it to satisfy To deal with the IR singularities that are present both in V 3 and in the integral of B 4 over Φ 4 , we use the FKS subtraction method [70]. We introduce a set of projecting functions, θ T m (Φ 4 ), that partition the phase space into nonoverlapping regions, such that m θ T m (Φ 4 ) = 1. In our case, this partition is effectively determined by the resolution variable, which is indicated by the superscript. The resulting partition must be such that each region m contains at most one collinear and one soft singularity. Then, we can write where Φ m rad denotes the radiative phase space describing a 1 → 2 splitting in each region. For each region m, we define a mapping that identifies which particle in Φ 3 is undergoing the 1 → 2 splitting, which generates the phase space point. It also unequivocally defines how the recoil is shared amongst the remaining particles in the event, which is needed to enforce total momentum conservation. Notice that our definition of the θ T m -functions leaves us the freedom to include phase space regions in the partition which do not contain any IR singularity. This freedom is in fact essential to be able to satisfy Condition 1 in eq. (2.32), namely, to ensure that in each phase space region the correct functional form for the resolution variable T 2 (Φ 4 ) is used. For example, in e + e − → qqg production, the region in which the quark and antiquark are closest to each other does not contain any QCD singularity. Nevertheless, it must be treated as a separate region in the phase space partition, since in this region the invariant mass between the quark and antiquark determines the value of T 2 .
With this notation, one can write

JHEP09(2013)120
If the region m contains an IR divergence, the FKS subtraction requires one to define the soft, collinear, and soft-collinear limits of Φ m rad , which we denote as Φ m,s rad , Φ m,c rad , and Φ m,cs rad , respectively, together with the resulting points in the 4-body phase space Φ m,s 4 , Φ m,c 4 , and Φ m,cs 4 . We can then write where θ s m encodes the soft limit of the θ T m -functions and we have used the fact that in the collinear and soft-collinear limits the θ T m -functions are trivially satisfied. Also, since in each of these limits T 3 ≡ 0, the θ(T 3 < T cut 3 ) functions are satisfied by construction. If m is not singular, in principle, no such subtraction is needed, and one could simply evaluate the 4-parton tree-level matrix element B 4 (Φ m 4 ). However, given that the integral of the subtraction counterterms over the whole phase space is known analytically for both massless and massive partons [17,70,71], The crucial point, discussed in section 2.2.2, is that our construction requires the phase space map that generates Φ m 4 (Φ 3 , Φ m rad ) to preserve the value of T 2 ; i.e., Comparing this to eq. (3.13), we see that the map Φ m 4 (Φ 3 , Φ m rad ) must be precisely the inverse of Φ T 3 (Φ 4 ) in the region m. In principle, this condition can be relaxed to only hold up to power corrections. Additionally, the map can fail to preserve T 2 in a region of phase space that gives a power-suppressed contribution to the cross section. The phase space maps used in the standard FKS implementations [16,70] were not designed to preserve the value of T 2 , and thus, they change its value by an O(T 3 /T 2 ) amount over a large region of phase space. 12 Since 4-parton events with T 3 < T cut 3 and T 2 > T cut 2 are the 11 These integrals can also easily be defined by restricting the integration of the FKS variable ξ up to some ξcut value. These are, however, not in direct correspondence to the partition of phase space we are considering. 12 Generically, an emission that takes a 3-parton event to a 4-parton event will change T2 by the scaling T2(Φ4) − T2(Φ3) ∼ T3.

JHEP09(2013)120
only real emission contributions included in the NLO calculation for the 3-jet cumulant, eq. (3.11), one can impose the restriction T cut 3 T cut 2 and use the standard FKS phase space maps. However, this hierarchy strongly restricts T cut 3 , and it is preferable to define a map that is specifically designed for our goals. We have constructed such a map, which preserves the exact value of T 2 up to power corrections, except in a region of phase space whose contribution to the cross section scales as O(T 2 /E cm ). In this region, the value of T 2 is altered by an O(T 3 /T 2 ) amount, meaning the net correction scales as O(T 3 /E cm ). Therefore, enforcing the much looser constraint T cut 3 E cm is sufficient to achieve our purposes. We postpone the detailed discussion of this map to a dedicated publication describing the implementation of Geneva.

Parton shower and hadronization
The phase space points Φ 2 and Φ 3 represent jet kinematics, which are defined by the jet resolution variable 2-jettiness. As discussed in section 2.4, we require that the parton shower does not change the underlying hard jet distribution so that the higher-order weights we calculate, dσ/dΦ 2 (T cut 2 ) and dσ/dΦ 3 (T 2 ), are correctly assigned. Without any constraints the parton shower will not preserve the value of T 2 . We address this problem in our current implementation in a physically motivated way. For small T 2 , the resummed singular jet resolution spectrum dominates and is determined up to power corrections of order λ ∼ T 2 /E cm . We require that for Φ ≥3 events, the change in T 2 due to showering, ∆T 2 , satisfies ∆T 2 /T 2 < λ. This represents a power correction to the 2-jettiness spectrum, which scales as 1/T 2 for small T 2 . For the 2-jet cumulant bin, we require that Φ 2 events, which have T 2 = 0 when unshowered, remain in the 2-jet bin after showering up to a power suppressed correction, with T 2 < T cut 2 (1+λ ). Here, the effect of a small nonzero λ induces a change to the shape of the distribution generated by Pythia that scales as a power correction and does not affect the formal accuracy of the spectrum. Formally, we work in the limit where λ and λ are effectively taken to zero. For λ = 0, the shower would be required to exactly preserve T 2 , making it maximally inefficient. Therefore, for the results shown in this section, we use small nonzero values λ = 2 λ = 0.05. We have checked that these are small enough to be in the asymptotic region where the results become independent of the precise value.
Furthermore, the shower must also be restricted to not change the NLO 3 result. This requires that, for 3-parton events, we effectively only allow showers to start from T cut 3 . Likewise, we limit the showering of 4-partons events down from their T 3 (Φ 4 ) value. This can be seen as a proxy for what would be the correct approach in a T N -ordered shower.
We use Pythia 8.170 with e + e − tune 1 for showering and hadronization. The choice of tune for e + e − data in Pythia affects both the time-like showering and hadronization model. However, since in our implementation we restrict the shower from changing the T 2 spectrum, the effect of changing the tune in Pythia primarily reflects the uncertainty from hadronization in Geneva. We have checked that this is also the case for observables other than T 2 by verifying that the effect of the tune on the showered Geneva predictions is very small compared to the change due to hadronization. The uncertainty from hadronization is associated with the nonperturbative contribution to the soft function in eq. (3.9) in our JHEP09(2013)120 framework and is not included in our event-by-event perturbative uncertainties, which are derived from the analytical resummation and fixed-order matching. A complete uncertainty analysis should also include uncertainties due to hadronization as well as due to the remaining amount of parton showering. As an indication of the size of the uncertainty from hadronization, we also show Geneva hadronized results using e + e − tune 3. The shift from the partonic to the showered results could be taken as a conservative upper limit on the remaining showering uncertainty.
It is important to note that we use the standard tunes in Pythia, without changing any internal parameters. Since in our approach the shower evolution in standalone Pythia 8 is substituted with higher-order resummation above the 2-jet resolution scale, we advocate that a separate tuning of Geneva + Pythia 8 should be employed to obtain the best results. This would also allow one to obtain meaningful estimates of hadronization and remaining showering uncertainties.

Validation using the jet resolution spectrum
Before comparing the Geneva prediction for various e + e − spectra to analytic predictions and LEP data, we first validate the implementation of our procedure to combine higherorder resummation and full NLO matrix-element corrections by using the jet resolution spectrum. At the level of the singly differential T 2 spectrum only, the standard approach to resummation achieves the same matching between resummation and fixed order by adding the nonsingular contribution to the resummed result. This provides a nontrivial crosscheck of the master formula and in particular validates the event-by-event theory uncertainties generated by Geneva. For each comparison in this section, we show the peak, transition, and tail regions, described in section 2.1, at the LEP center-of-mass energy E cm = 91.2 GeV. In all cases, the error bars or bands on the Geneva histograms are built from its event-by-event perturbative uncertainties. The statistical uncertainties from Monte Carlo integration are much smaller and are not shown.

Partonic results
The analytic resummed T 2 spectrum is shown in figure 2 at successively higher orders: NLL, NLL +LO 3 , and NNLL +NLO 3 (see table 2 for the order-counting definitions). The perturbative uncertainties are generated by using the same profile scale variations employed in Geneva and discussed in section 3.1.1. The theory uncertainties decrease at increasing order and demonstrate excellent convergence at all values of T 2 . Below T 2 < 0.5 GeV, we enter the purely nonperturbative region, and the scale uncertainties diverge since even resummed perturbation theory breaks down. In the far tail, the scale uncertainties also grow rapidly, which reflects missing higher fixed-order corrections. The uncertainties in the NNLL +NLO 3 prediction diverge past the 3-parton endpoint at E cm /3, where the fixedorder prediction is only correct at leading order for 4 partons. In the transition region, there is a smooth interpolation between the resummation and fixed-order regions.
In figure 3, we compare the partonic T 2 spectrum from Geneva with T cut 2 = 1 GeV to the analytic resummed results from figure 2. To illustrate the interpolation between resummed and fixed-order results, we also show the pure resummed results at NLL and JHEP09(2013)120  NNLL and the pure fixed-order contribution at LO 3 and NLO 3 . The latter are calculated using Event2 [72,73], which serves as an independent crosscheck of our NLO 3 implementation. Using the NLL +LO 3 resummation of T 2 and the LO 3 fixed-order contribution as inputs to our master formula for the spectrum in eq. (3.3), the dσ FO /dΦ 3 and dσ resum /dΦ 2 dT 2 FO contributions exactly cancel for the T 2 spectrum. As a result, we see precise agreement between Geneva and the analytic NLL +LO 3 result in figure 3(a)-3(c) in the peak, transition, and tail regions. This result agrees well with the pure NLL resummed contribution in the peak, while in the tail, it is consistent within uncertainties with the LO 3 result, where the latter clearly underestimates the full perturbative uncertainties.
At next higher order, using as inputs to the master formula the NNLL +LO 3 resummation of T 2 and the NLO 3 fixed-order calculation, we see that the central value and eventby-event uncertainties in Geneva agree very well with the full analytic NNLL +NLO 3 resummed prediction in the peak and transition regions, as shown in figures 3(d) and 3(e). In the tail region, figure 3(f), Geneva has significantly smaller uncertainties of the same size as the pure fixed-order contribution. This is because there is a substantial cancellation between singular and nonsingular contributions in this region, which is incorporated JHEP09(2013)120  differently in the analytic resummation and the master formula at NNLL +NLO 3 . For the former, the nonsingular α 2 s contributions are added. This preserves the absolute size of residual resummation uncertainties, which are very small relative to the singular contributions but large relative to the total result after cancellation. In the master formula in eq. (3.3), the nonsingular contributions are incorporated multiplicatively through the ratio of dσ FO /dΦ 3 and dσ resum /dΦ 2 dT 2 FO . This preserves the relative size of residual resummation uncertainties, thus leading to much smaller absolute variations when compared to the final result. Comparing the Geneva prediction with the pure NNLL resummed and NLO 3 fixed-order results, we see that the master formula precisely interpolates as expected between the fixed-order and resummation regions, with the transition region properly describing the transition between the two, including uncertainties.
Combining the exclusive 2-jet cross section with the integral of the inclusive 3-jet cross section, the Geneva prediction at NNLL +NLO 3 formally reproduces the total inclusive cross section at NLO. Numerically, we have σ NLO tot = 44.1 ± 0.2 nb. With T cut 2 smeared between 0.5 − 1 GeV, the total inclusive cross section in Geneva is σ Geneva tot = 42.5 ± 1.6 nb, where the uncertainties are given by integrating over the different profile scale variations. The central value is 3.8% low and agrees within the uncertainties of ±3.8%. The uncertainty in Geneva that comes from integrating the spectrum over T 2 > T cut 2 , as in eq. (2.4), is much larger than the fixed-order uncertainty. The reason is that, at any given point in the JHEP09(2013)120 spectrum, but especially in the peak region, the relative uncertainties, reflecting both shape and normalization, are larger than in the total cross section. Hence, when integrating the spectrum to obtain the total cross section, the uncertainties in the spectrum must cancel each other, meaning there is a negative correlation in the uncertainties between different regions in the spectrum. When the resummation and matching to fixed order is performed for the spectrum, this correlation and cancellation is numerically not exact for the total cross section. This is a well-known limitation of analytic resummation [52]. In fact, the result from Geneva is completely consistent with the inclusive cross section obtained using the analytic resummed result in eq. (2.4) with σ(T cut 2 ) calculated at NNLL +LO 3 and dσ/dT 2 calculated at NNLL +NLO 3 . In this case, with T cut 2 = 1 GeV, the central value is 3.5% low with uncertainties of ±3.7%. One way to solve this problem would be to enforce a (highly nontrivial) constraint on the profile scale variation to reproduce the required correlation exactly, in which case the total cross section would come out exactly right. In practice, a simpler way to enforce this is to compute the result for the resummed cumulant as the difference between the total cross section and the integrated resummed spectrum. (This is similar in spirit to the method proposed in refs. [31,32,39].) Since our focus in this paper is the differential spectrum, which serves as the primary input to the Monte Carlo, rather than the total cross section, we leave this for future improvement.

Showered results
Next, we validate our interface with the parton shower. In figure 4, we compare the NNLL +NLO 3 partonic and showered Geneva predictions with T cut 2 smeared between 0.5 − 1 GeV. We also show the analytic resummed NNLL +NLO 3 and pure fixed-order NLO 3 spectra for comparison. Before showering, the cumulant dσ/dΦ 2 (T cut 2 ) is in the T 2 = 0 GeV bin, and we see the effect of the smeared T cut 2 on the spectrum in the Geneva partonic histogram in figure 4(a). The parton shower generates emissions inside the 2jet bin, which fills out and determines the shape of the Geneva showered result in the region below T cut 2 and agrees remarkably well with the analytic resummed spectrum below the cut. While the shape of the spectrum here is determined only by Pythia, the cross section below T cut 2 is still accurate to NNLL +LO 3 . We can see this explicitly in figure 5 from the separate contribution of 2-, 3-, and 4-parton events before and after showering for the central value in the peak region. The shape of the 2-parton showered histogram is determined by Pythia, and the area under the histogram is the cumulant dσ/dΦ 2 (T cut 2 ) calculated at NNLL +LO 3 . The relative contribution of 3-parton and 4-parton events is determined by T cut 3 = 2 GeV, for which the 4-parton contribution is well behaved, giving 15% of the total cross section and no large cancellation with 3-parton events. These contributions all combine smoothly to generate the total Geneva showered result.
The action of the shower on 3-parton and 4-parton events, which make up the spectrum above T cut 2 , is restricted to not change T 2 by more than a power suppressed amount λ T 2 , as discussed in section 3.1.3. This controls the allowed shift from the Geneva partonic to showered histograms in figure 4. We can see that there is excellent agreement, including uncertainties, between the two in the peak and transition regions. This validates that, with our choice of λ, the higher-order accuracy of the resummed T 2 spectrum is not com-JHEP09(2013)120

JHEP09(2013)120
promised by the shower. (Increasing λ, we do observe, at some point, a shift of showered results away from partonic.) The showering does shift the T 2 spectrum in the far tail away from the partonic result, which matches the NLO 3 curve, as can be seen in figure 4(c). This is allowed, since our partonic prediction in this region becomes only leading order for 4 partons.

Hadronized results and comparison to data
The full prediction for the jet resolution spectrum is obtained by turning on the hadronization in Pythia. This gives rise to a shift in the T 2 spectrum, shown in figure 6, where "default" refers to the default running parameters α s (m Z ) = 0.1135 and Pythia e + e − tune 1. As discussed in section 3.1.3, we use the standard Pythia 8 tunes without modifying any internal parameters. For comparison, we show the Geneva hadronized result for tune 3 with our default α s value as well as for tune 1 with the world average value α s (m Z ) = 0.1184. We also show a comparison to experimental data from ALEPH [74] and OPAL [75]. We only show ALEPH data in the tail since the OPAL data in this region is sparse. These measurements are fully corrected to the particle level, allowing us to directly compare to our hadronized predictions. Since the data are normalized to the total cross section, we rescale them to the total NNLO cross section and convert from thrust T to . This allows us to directly compare the data to the absolute cross section predictions in Geneva, unlike a comparison between normalized spectra, which would only test the shape. The Geneva prediction at the default values agrees impressively well with the data within uncertainties across the peak and transition regions and into the tail. The difference in the far tail is expected since here fixed-order contributions beyond LO 4 are important and are not yet included in our results. The partonic Geneva prediction does not include nonperturbative effects in the soft function of O(Λ QCD /T 2 ), nor power corrections of the form O(Λ QCD /E cm ). Since we strongly constrain the action of the Pythia parton shower to not change the analytic resummed NNLL +NLO 3 result, as discussed in section 3.1.3 and demonstrated in figure 4, we expect the hadronization in Pythia to supply these missing nonperturbative effects. In effect, Pythia provides a well-tested model of the nonperturbative soft function in eq. (3.9). We show the hadronized Geneva result with Pythia e + e − tune 3 at the central scale in figure 6 as a measure of the uncertainty from hadronization. Tune 3 turns out to give a smaller shift due to hadronization than tune 1, which makes a significant difference in the peak below 3 GeV, where nonperturbative corrections are O(1) and depend on the details of the hadronization model. In the transition and tail regions, we see a smaller difference, with tune 3 being systematically lower than tune 1. This is consistent with the fact that the transition and tail regions are sensitive only to the first nonperturbative power correction in the soft function of O(Λ QCD /T 2 ).
There is an important interplay between the effect of hadronization and the value of α s (m Z ), as discussed in ref. [52], where a simultaneous fit to α s (m Z ) and the first nonperturbative correction to the soft function of O(Λ QCD /T 2 ) was carried out. Generically, larger nonperturbative corrections shift the partonic spectrum to larger values of T 2 , while a smaller value of α s (m Z ) shifts the 2-jettiness spectrum downward. This gives rise to JHEP09(2013)120 (d) Ratio of Geneva to data Figure 6. The showered NNLL +NLO 3 Geneva prediction with and without hadronization using the default values Pythia 8 e + e − tune 1 and α s (m Z ) = 0.1135 compared to data from ALEPH [74] in the (a) peak, (b) transition, and (c) tail regions and to OPAL [75] in the peak and transition regions. The ratio of Geneva predictions to the ALEPH data is shown in (d). Also shown is the Geneva prediction at the central scale with α s (m Z ) = 0.1184 and e + e − tune 3. compensating effects. Since tune 3 gives a smaller shift due to hadronization than tune 1, the combination of tune 3 and α s (m Z ) = 0.1135 gives an estimate of the lower bound on the combined uncertainty of these two effects in the transition and tail regions, while the combination of tune 1 and α s (m Z ) = 0.1184 gives an estimate of the upper bound. This is illustrated very well in the ratio of Geneva to ALEPH data in figure 6(d). Both are, however, still within the perturbative uncertainties from Geneva across most of the transition and tail regions.
We have also checked that the nonperturbative shift from Pythia tune 1 is of similar size as expected from the fit results in ref. [52]. This is consistent with the fact that it gives a good description of the data when used together with their fitted value of α s (m Z ). Hence, we use tune 1 with α s (m Z ) = 0.1135 as the default since it agrees best with the data in the peak and provides a consistent description of the data across larger values of T 2 .

Predictions for other event shapes
In this section, we present Geneva's predictions for a variety of dijet event shape variables. Examining observables other than the jet resolution variable we use as input serves to validate our master formula at the fully differential Φ 2,3 level [see eq. (3.2)] rather than its projection onto the T 2 spectrum [see eq. (2.24)]. Event shapes are particularly useful to consider because there exist both higher-order resummed results and precision LEP data with which we can compare.
By construction, Geneva correctly predicts other observables at NLO 3 , while maintaining the correct inclusive cross section. However, as discussed in section 2.2.3, it is an important open question to what extent the NNLL T +NLO 3 resummation of the T 2 spectrum increases the accuracy of resummed predictions for the other observables (beyond the partial NLL order naively expected by interfacing with the parton shower). While other observables will not be predicted at the same resummed order as T 2 , the accuracy of the predictions for event shapes is expected to increase as a function of their correlation with 2-jettiness. The comparison of Geneva to the higher-order analytic resummation of event shapes plays a crucial role in numerically testing the accuracy achieved in our approach and validating the event-by-event perturbative uncertainties.
We present results for the C-parameter [76][77][78], heavy jet mass (ρ) [79,80], and jet broadening (B) [81,82] event shapes. These are defined as follows: wheren T is the thrust axis and is used in heavy jet mass to divide the event into two hemispheres, hemi 1,2 , with respect to which the masses M 1,2 are measured. C, ρ, and B provide a useful range of event shapes to compare to since their resummation structure is increasingly different from that of T 2 . The resummation of C-parameter is precisely the same as T 2 to NLL and has the same convolution structure as eq. (3.8) beyond. Heavy jet mass has a different convolution structure from T 2 . Both ρ and T 2 are projections of the same doubly differential spectrum dσ/dM 2 1 dM 2 2 , where T 2 is related to the sum and heavy jet mass to the maximum of the hemisphere masses. Of the event shapes we consider, jet broadening is most different from T 2 ; it measures momentum transverse to the thrust axis and, in the dijet limit, is sensitive to the recoil of the thrust axis due to soft emissions [83], unlike T 2 . This complicates the higher-order resummation of jet broadening, which was only recently extended to NNLL B [84] and gives a logarithmic structure that is very different from T 2 . As a result, jet broadening provides a highly nontrivial test of the accuracy and theory uncertainties of the Geneva prediction.
For each of these observables, we compare to analytic resummed predictions as well as the NLO 3 fixed-order contribution from Event2. We present new results for the analytic JHEP09(2013)120 resummation of C-parameter at NNLL C +NLO 3 , extending the previous NLL C resummation of [85]. 13 Note the subscript on the order of resummation indicates the observable for which analytic resummation was carried out. Since resummed results for jet broadening do not exist at NNLL B , we compare to the highest available resummation NNLL B +LO 3 , where we use the results of [84], which we extend to include fixed-order matching that is necessary to describe the tail and transition regions. Finally, for heavy jet mass, N 3 LL ρ resummed results exist [51]; however, we show the NNLL ρ +NLO 3 resummation since this is consistent with the highest T 2 resummation we use.
It is important to note that all running parameters were set based on the T 2 spectrum alone, and no further optimization was carried out for other observables. This ensures that our results for other observables are true predictions of the Geneva framework.

C-parameter
In figure 7, we show the Geneva prediction for C-parameter both at the partonic level and showered, using NLL T +LO 3 resummation as input to our master formula in figures 7(a) and 7(b) and at next higher order NNLL T +NLO 3 in figures 7(c) and 7(d). We compare this to the analytic resummed C-parameter prediction at the same order as the T 2 resummation we input, as well as one order lower. The comparison of the Geneva prediction at different orders in the peak and transition regions is useful because it highlights the features of resummation that are consistently captured by our implementation. In the tail region, figure 7(e), where the comparison to the NLO 3 fixed-order result is most relevant, we only show our highest order NNLL T +NLO 3 Geneva result.
We see the effect of the cut on 2-jettiness of T cut 2 = 0.5 − 1 GeV up to C = 0.066 in the partonic prediction from Geneva in figures 7(a) and 7(c) since C ≤ 6T 2 /Q [85]. Interfacing with the shower generates emissions inside the jets and fills out the region below C = 0.066. The action of the parton shower is restricted based on the constraints on T 2 discussed in section 3.1.3. This effectively constrains the C-parameter distribution as well, giving very little change from the partonic to showered predictions at both resummation orders except in the multijet region of the far tail. Here, the constraints on the shower are looser, reflecting the fact that our prediction is correct at LO 4 . The size of the shift from the partonic to the showered result in the peak and transition regions is a measure of the correlation between the C and T 2 event shapes, where, although the two differ beyond NLL, their logarithmic structure is the same. It is worth noting however, that despite the similarity of the resummation structure between C and T 2 , the shape of the C-parameter spectrum is very different, with the singular terms dominating the nonsingular for a much larger region of the spectrum.
One might naively expect that the accuracy of the resummation achieved in Geneva for any observable other than T 2 would only be the partial NLL of the parton shower. However, it is clear that the Geneva prediction at NLL T +LO 3 in the peak and into the transition region, figures 7(a) and 7(b), is much more consistent with NLL C +LO 3 than 13 We thank Vicent Mateu and Iain Stewart for pointing out to us the relationship between thrust and C-parameter in SCET.
(e) Geneva NNLL T +NLO 3 tail region Figure 7. The C-parameter partonic and showered Geneva predictions are shown compared to the analytic resummation of C-parameter at different orders. The Geneva result at NLL T +LO 3 is compared to NLL C and NLL C +LO 3 in (a) and (b). In (c) and (d), the Geneva prediction at one order higher, NNLL T +NLO 3 , is compared to NLL C +LO 3 and NNLL C +NLO 3 , while in the tail (e), we also show the fixed-order NLO 3 prediction from Event2.
NLL C resummation, both in its central value and also in the size of the perturbative uncertainties it predicts. This appears to hold even in the peak region below C ∼ 0.05, where the parton shower determines the shape of the spectrum. Going to one higher order JHEP09(2013)120 in figures 7(c) and 7(d), we see that the same pattern holds: the Geneva prediction is consistent with the higher-order NNLL C +NLO 3 resummation rather than NLL C +LO 3 , including uncertainties. This is particularly clear in the peak region where the central values of the two analytic resummation orders are significantly different and Geneva tracks the NNLL C +NLO 3 prediction. The convergence of the Geneva result for C-parameter from NLL T +LO 3 to NNLL T +NLO 3 demonstrates the consistency of the Geneva implementation including the event-by-event uncertainties for this observable. Although the accuracy of the Geneva prediction for C-parameter is not formally of the same order as the T 2 resummation we used as input to the master formula, the fact that it matches the analytic C-parameter resummation remarkably well both at NLL C +LO 3 and at NNLL C +NLO 3 shows that numerically the accuracy achieved is very close.
The Geneva uncertainties in the transition region start to shrink relative to the analytic resummation as we interpolate to the fixed-order NLO 3 result. In the tail region, the partonic Geneva prediction matches smoothly to the fixed-order NLO 3 result past the Sudakov shoulder at C = 0.75 [86], demonstrating the validity of the multiplicative implementation of dσ/dΦ 3 (T 2 ) in eq. (3.3) in this limit.
The Geneva prediction including hadronization with the default running values of Pythia e + e − tune 1 and α s (m Z ) = 0.1135 is shown in figure 8 compared to ALEPH and OPAL data rescaled to the NNLO inclusive cross section. Geneva agrees with the data remarkably well across the entire distribution up to the multijet region in the tail. We show the effect of α s (m Z ) = 0.1135 with tune 3, which gives a smaller correction from hadronization than tune 1, as seen from the size of the shift from the Geneva unhadronized result to the hadronized in figure 8. We also show the Geneva prediction at the central scale using the world average α s (m Z ) = 0.1184 and tune 1. These two combinations provide an estimate of the upper and lower bounds on the combined uncertainty of the nonperturbative effect and α s (m Z ) value in the transition and tail regions, as discussed in section 3.2. The ratio of the Monte Carlo to data in figure 8(d) shows that they are both largely within the perturbative uncertainties from Geneva in these regions. Of the values we consider, the default tune 1 with α s (m Z ) = 0.1135 gives the best agreement with the data across the C spectrum and is consistent with our findings for the T 2 distribution.

Heavy jet mass
The Geneva prediction for heavy jet mass is shown in figure 9, where we compare the partonic and showered results using NLL T +LO 3 resummation in the master formula in figures 9(a) and 9(b) to the analytic resummation of ρ at NLL ρ and NLL ρ +LO 3 . In figures 9(c) and 9(d), we show the same results at one order higher, comparing NNLL T +NLO 3 Geneva results to NLL ρ +LO 3 and NNLL ρ +NLO 3 analytic ρ resummation. In the tail, we show only the highest-order Geneva and resummed results along with the pure fixed-order NLO 3 contribution, since this is sufficient to demonstrate the behavior in this region.
In the peak region, figures 9(a) and 9(c), we see the effect of T cut 2 on the partonic ρ spectrum up to ρ = 1 GeV/E cm = 0.011, which is smoothly filled out by interfacing with the parton shower. The Geneva showered prediction in figure 9(a) shows impressive agreement with the NLL ρ +LO 3 resummed result in the peak region, including uncertainties. The     improvement in accuracy of the Geneva prediction for heavy jet mass over the partial NLL resummation provided by the parton shower is clear. Going to one higher order in figure 9(c), the Geneva prediction is more consistent with the NNLL ρ +NLO 3 , with which it agrees within uncertainties, rather than the NLL ρ +LO 3 result.
The perturbative uncertainties of the Geneva showered prediction are larger than those at NNLL ρ +NLO 3 and smaller than at NLL ρ +LO 3 . This is consistent with the fact that, while the Geneva prediction for ρ is not formally of the same order as the T 2 resummation that is input into the master formula, there is a gain in accuracy going from Geneva at NLL T +LO 3 to NNLL T +NLO 3 that is transferred to the ρ prediction.
In the transition region, both at NLL T +LO 3 and NNLL T +NLO 3 , adding the parton shower gives rise to a larger shift from the partonic spectrum than for the C-parameter, because heavy jet mass is less correlated with T 2 than C. This shift is necessary to obtain agreement with the NNLL ρ +NLO 3 resummation within uncertainties in figure 9(d). The partonic Geneva prediction in this region is more consistent with NLL ρ +LO 3 analytic re-   (e) Geneva NNLL T +NLO 3 tail region Figure 9. The heavy jet mass partonic and showered Geneva predictions are shown compared to the analytic resummation of ρ at different orders. The Geneva result at NLL T +LO 3 is compared to NLL ρ and NLL ρ +LO 3 in (a) and (b). In (c) and (d), the Geneva prediction at one order higher, NNLL T +NLO 3 , is compared to NLL ρ +LO 3 and NNLL ρ +NLO 3 , while in the tail (e), we also show the fixed-order NLO 3 prediction from Event2. summation, which is higher than the NNLL ρ +NLO 3 result. By restricting the change in T 2 due to the shower, we are constraining the sum of the hemisphere masses, M 2 1,2 in eq. (3.20). For a given value of T 2 , ρ is largest when either hemisphere mass is zero, and so ρ = T 2 . Adding the parton shower tends to make this mass nonzero (while constraining the sum)

JHEP09(2013)120
and therefore gives an overall shift of the spectrum to lower values of ρ. In the tail region, the partonic spectrum interpolates to the fixed-order NLO 3 result, as expected, with the shower giving rise to a larger shift in the multijet region where our constraints are looser.
In figure 10, we show the showered Geneva prediction with and without hadronization, with our default parameters. As before, we compare to ALEPH and OPAL data, which shows impressive agreement with the data within uncertainties across all three regions of the ρ spectrum, with the expected deviation in the multijet region of the far tail. As discussed previously, tune 3 with α s (m Z ) = 0.1135 and tune 1 with α s (m Z ) = 0.1184 provide bounds on the estimate of the combined uncertainty from these two inputs. It is interesting to note that heavy jet mass is relatively insensitive to hadronization in the transition and tail regions. This is demonstrated by the shift from the shower-only to the fully hadronized result in figures 10(b) and 10(c), as well as the small change in the default central value from using tune 1 to tune 3 above ρ ∼ 0.1 seen in figure 10(d). This breaks the coupling in some respect between the nonperturbative effects and the value of α s in this region and suggests that α s (m Z ) = 0.1135 provides better agreement with the data.

Jet broadening
Finally, we turn to the results of Geneva for jet broadening, which is the most orthogonal event shape to our jet resolution variable that we consider. In figure 11, we show the Geneva partonic and showered results using NLL T +LO 3 resummation in figures 11(a) and 11(b) and NNLL T +NLO 3 resummation in figures 11(c) and 11(d). We compare these to the analytic NLL B and the best available NNLL B +LO 3 resummed prediction. Note that we would like to compare the NLL T +LO 3 resummation in Geneva to the resummation of B at the same order. However, since going from NLL B to NNLL B (which incorporates α 2 s ln B terms into the resummation) is a comparatively small effect in this case, we will find it useful to compare the NLL T +LO 3 prediction with the analytic NNLL B +LO 3 result. In the tail, we compare to the fixed-order NLO 3 result from Event2, which is the more relevant comparison in this region.
The effect of the cut on T 2 extends up to B 0.055 in the peak region and is smoothly removed by attaching the parton shower. Both at NLL T +LO 3 and NNLL T +NLO 3 , there is a significant shift induced by the parton shower across the jet broadening spectrum toward larger values of B. The size of this shift is a measure of the lack of correlation between B and T 2 , the variable used to constrain the parton shower, and is therefore progressively larger for C, ρ, and B, as we have seen.
As discussed in section 2.2.3, in the IR limit where both T 2 , B → 0, one might expect to see some improved accuracy of the Geneva prediction over the partial NLL of the parton shower, since the higher-order resummation of T 2 provides a better description in this region. This is consistent with figure 11(a), where the showered Geneva prediction agrees well with the NNLL B +LO 3 resummed result, including uncertainties in the peak region. In the transition region in figure 11(b), the central value of the Geneva showered prediction agrees with the NNLL B +LO 3 resummed result within uncertainties; however, the uncertainties from Geneva are smaller than the corresponding analytic ones, which suggests that in this region they may be underestimated.   In the far transition region and into the tail, the uncertainties in Geneva generically are smaller than the corresponding analytic resummation and of order the NLO 3 scale variation, as seen for example in the T 2 , C, and ρ spectra. This difference arises because Geneva multiplicatively interpolates to the fixed-order result, while the analytic resummation does an additive matching, as discussed in section 3.2.1. For an observable such as jet broadening, the lack of correlation with T 2 means that larger values of the T 2 spectrum contribute at smaller values of B. This can lead to an underestimate of the uncertainties from Geneva at intermediate values of B where the resummation is still important.
Going to higher order in Geneva in figures 11(c) and 11(d), the uncertainties of the showered prediction decrease and overlap with the NNLL B +LO 3 uncertainties over much of the peak and transition regions. It would be interesting to compare the Geneva prediction to the next higher-order analytic resummed jet broadening prediction to numerically test the accuracy achieved; however, this is not yet available. Determining the formal accuracy of the Geneva prediction for a given observable and systematically including the  uncertainty associated to the lack of correlation with T 2 are next steps that we leave for future work. As mentioned in section 3.1.3, one possibility would be to include the size of the shift from the partonic to the showered results as a conservative estimate of the uncertainty due to the remaining showering.   The partonic Geneva jet broadening prediction interpolates smoothly to the fixedorder NLO 3 result in the tail, with uncertainties that match the fixed-order result in this region. As for the other observables, this validates the behavior of the fully differential master formula in eq. (3.3) in this limit.
In figure 12, we show the hadronized Geneva prediction for jet broadening compared to data from ALEPH and OPAL, which shows good agreement within uncertainties across the peak and transition regions and is low as expected in the far tail. The uncertainty from the Pythia tune and value of α s are indicated by the central values of the tune 3 and α s (m Z ) = 0.1184 histograms, which agree within the perturbative uncertainties of the default Geneva prediction across most of the transition and tail regions. As for the other observables, we see better agreement in the peak (below B ∼ 0.1) with Pythia tune 1 and α s (m Z ) = 0.1135.

JHEP09(2013)120 4 Application to hadronic collisions
In this section, we apply the framework developed in section 2 to hadronic collisions and present first results from the implementation in the Geneva Monte Carlo. To accommodate hadrons in the initial state, special consideration is required for each component of our master formula. Our goal is to demonstrate that the same methods can be applied in a hadronic environment to obtain a consistent description at the next higher perturbative accuracy. We use Drell-Yan production pp → Z/γ * → + − as a concrete example to study the framework, deferring a detailed comparison with Tevatron and LHC data to later work.
In hadronic collisions, N -jettiness can be used as a jet resolution variable, with the observable taking the initial states into account. The theoretical framework exists to resum N -jettiness at hadron colliders, and this resummation has been applied to several processes [69,87,88]. Similarly, the techniques required to perform the next-to-leadingorder calculations at hadron colliders are known. A phenomenological study additionally requires Geneva to be interfaced with a parton shower and hadronization model that includes multiple parton interactions.
In the Drell-Yan example, the 0-jet resolution variable is beam thrust, T 0 , which is the analog to 2-jettiness, T 2 , used in the e + e − application. 14 The resummation of beam thrust is performed to NNLL, and the jet multiplicities at fixed order are calculated at NLO 0 and LO 1 . The prediction of Geneva is compared to the analytic resummation of T 0 at NNLL+LO 1 as a demonstration that the program correctly describes the matching between 0-and 1-jet multiplicities. Finally, we discuss how the accuracy of these ingredients can be improved and the challenges present in applying Geneva to hadron collisions.
Beam thrust is defined as a sum of contributions from particles in the final state [46], where the observable is evaluated in the center-of-mass frame of the hard partonic collision. The n a,b are light cone vectors along the beam (ẑ) axis, with n a = (1,ẑ) and n b = (1, −ẑ). Beam thrust can be evaluated in any frame by performing a longitudinal boost on eq. (4.1).
With the addition of more final-state jets, the N -jettiness definition can be generalized from the 0-jet case, As for beam thrust, this observable is evaluated in the partonic center-of-mass frame. The n i = (1,n i ) for i = 1, . . . , N are light cone vectors along the jet directions. Note that the above definition of T N is a simple extension of the observable for e + e − collisions, which has no contribution from the beam directions but is otherwise identical. JHEP09(2013)120

Master formula and ingredients for hadronic collisions
As in the e + e − case, the master formula for the cross section in Geneva is given by eq. (2.27), eq. (2.29), and eq. (2.36). To match the 0-and 1-jet multiplicities for a general process, the master formula is where the 0-jet cumulant, dσ/dΦ 0 (T cut 0 ), and the 1-jet spectrum, dσ/dΦ 1 (T 0 ), are Φ 0 is the phase space for the hard scattering that produces the 0-jet final state, and Φ 1 includes the additional phase space for the final-state jet.
In Geneva, the contributions to these cross sections are calculated separately for each parton subprocess. Because the fixed-order matching in the 0-jet cumulant is performed additively, the net weight in the 0-jet cumulant after summing over events is the same as the flavor-summed cumulant. In the 1-jet spectrum, the fixed-order matching is performed multiplicatively, meaning the sum over all events for a given T 0 has a different cross section than if we used flavor-summed components for the different terms in the matching formula. The two approaches agree up to higher-order corrections, but the former approach is natural in the Monte Carlo. In the following subsections, we will discuss how the resummed and fixed-order contributions to the master formula are obtained.

Resummation
Like T 2 for e + e − collisions, beam thrust can be factorized in SCET and the resummation can be carried out systematically to higher orders. The factorized beam thrust spectrum for Drell-Yan is given by [34] dσ where Q is the dilepton invariant mass, and Φ 0 is the phase space for the qq → + − hard scattering. The momentum fractions x a,b are defined in terms of the total rapidity Y of the final state from the hard scattering, Comparing eq. (4.5) to the e + e − analog, eq. (3.8), it is clear that the factorization theorems are structurally similar. The chief difference is that, while the jet functions in eq. (3.8) parametrize the collinear evolution of final-state jets, the beam functions in eq. (4.5) parametrize the collinear evolution of the incoming partons as well as the nonperturbative process of extracting high-energy partons from the proton. The beam functions can JHEP09(2013)120 be further factorized into a convolution between the parton distribution functions f j and perturbatively calculable Wilson coefficients I ij [34,89], Note that due to initial-state radiation, the x a,b are distinct from the Bjorken variable ξ appearing in the convolution above that gives the momentum fraction of the energetic partons that are liberated from the proton. The sum over partons i, j in eq. (4.5) is a sum over flavor singlet quark-antiquark combinations, such as uū orbb. For each flavor, the beam functions are different, and the hard function H ij differs for up-and down-type quarks due to the different electroweak couplings to the gauge boson. This flavor sum is an important consideration when implementing the master formula in Geneva, since the Monte Carlo generates events for each flavor combination, and the flavor sum is performed in the sum over events.
For processes with final-state jets, the extension of the beam thrust factorization theorem to N -jettiness is known and has the schematic form [46,47] The trace is over the nontrivial color structures that can exist in the hard and soft functions. Additionally, there is a sum over parton channels for the hard scattering that are labeled by the index κ. The additional jets are associated with additional collinear sectors in SCET, and the factorization theorem reflects this by including additional jet functions. The soft function also changes to account for the soft radiation between the final-state jets and the initial-state radiation from the colliding partons. The factorization theorems in eqs. (4.5) and (4.8) can be used to perform the resummation for both the spectrum and cumulant. Although these factorization theorems directly describe the spectrum in T 0 or T N , they can be integrated over the observable to obtain the cumulant. The perturbative part of each function in the factorization theorem is calculable, and for many processes the functions are known to high order. Each function is associated with a scale that is connected to the physical degrees of freedom that the function describes. As in the e + e − case, renormalization group evolution resums the large logarithms of ratios of these scales (see table 2).

Fixed order
Following eq. (4.3), we need to define the 0-jet cumulant dσ/dΦ 0 (T cut 0 ) and the 1-jet spectrum dσ/dΦ 1 (T 0 ). At the order we are interested in, the 1-jet spectrum will be given by the tree-level cross section B 1 (Φ 1 ) for the process pp → Z/γ * j → + − j. The 0-jet cumulant is given by

JHEP09(2013)120
where the corresponding parton distribution functions have been included into the definitions of the Born, virtual, and real emission cross sections, In addition, in order to account for the incomplete cancellations of initial-state collinear singularities, we have included one collinear counterterm G a,b for each initial-state parton, Assuming the UV divergences of V 0 have already been taken care of by a proper renormalization procedure, the remaining divergences present in B 1 , V 0 , and G a,b are of IR origin. We handle these divergences with the FKS subtraction procedure. After having partitioned the phase space into nonoverlapping regions m, which contain at most one collinear and one soft singularity, by means of a set of θ T m -functions, the final formula, including the subtraction counterterms, is As mentioned in section 3.1.2, we choose to partition the phase space by means of θ T m -functions that depend on the jet resolution variable. It is therefore crucial to evaluate the jet resolution variable and the subtraction in the same frame. This ensures the proper cancellation of IR singularities by subtraction counterterms. The preferred frame for the fixed-order calculations is the partonic center-of-mass frame, since the subtraction is most naturally expressed in terms of variables defined in that frame. Also, the jet resolution variable, eq. (4.1), is defined in this frame and resummation can be performed in it. Therefore, our approach is to perform the entire calculation in the partonic center-ofmass frame. At this point, all the ingredients of eq. (4.12) are known and available in the literature [17,70]. Note that additional complications arise when extending this construction to higher multiplicities because one must use a map that preserves the value of T 0 (up to power corrections), just as for T 2 in the e + e − → 3 jet case. For our pp → Z/γ * → + − study, this problem can be avoided since we currently include only up to one extra parton and, consequently, the T 0 (Φ 1 ) value is well defined. In order to obtain T 0 (Φ N ) with N > 1 in general, this issue may be addressed in a similar fashion to what has been done for e + e − .

Application to Drell-Yan production
We study Drell-Yan production in pp collisions at E cm = 8 TeV in Geneva, sampling the invariant mass Q of the + − pair around the Z pole between M Z − 10Γ Z and M Z + 10Γ Z , JHEP09(2013)120  where M Z = 91.1876 GeV is the mass of the Z and its width is Γ Z = 2.4952 GeV [57]. The dominant contribution in this range of Q comes from Z exchange, although the photon does contribute. Profile scales identical to those used in the e + e − T 2 resummation are used, which is justified since the logarithmic structure is the same between the two observables. The resummation is turned off just above T 0 ∼ M Z /2, and for greater T 0 , the spectrum reproduces the fixed-order distribution.
In figure 13, we show the analytic beam thrust resummation at NLL and NNLL+LO 1 in the peak, transition, and tail regions. In the peak and transition regions, the resummed result converges well. At the end of the transition region and in the tail region, the pure NLL resummed distribution goes to 0 as the resummation is turned off, but the NNLL+LO 1 distribution moves into the fixed-order result.
Implementing the Drell-Yan process in Geneva allows us to study the feasibility of the multiplicative matching for the spectrum in eq. (4.4) and compare with the analytic resummed distribution. We show this comparison in figure 14, where the analytic curve is evaluated at NNLL+LO 1 . Additionally, we show the NNLL and LO 1 distributions separately. Overall, the partonic Geneva distribution is quite close to the NNLL+LO 1 distri-JHEP09(2013)120  bution, both in terms of the central values and the size of uncertainties. In the peak region, the spectra match closely and agree well with the pure NNLL resummed distribution. At the low end of the transition region, the resummed spectra are still in fair agreement while moving to higher T 0 values the Geneva partonic prediction and the NNLL+LO 1 distributions begin to systematically deviate from the NNLL distribution. This deviation arises from the LO 1 nonsingular terms that are present in the matched spectrum but absent in the pure resummed one. In the tail region, the Geneva partonic and NNLL+LO 1 predictions move into the LO 1 spectrum. After the resummation has been turned off, these spectra match the LO 1 precisely in central value and uncertainties, as expected. Note that in the corresponding comparison in the e + e − case, shown in figure 3(a)-3(c), the analytic and the partonic Geneva distributions are in closer agreement because the resummed components of the multiplicative matching in eq. (2.36) include the nonsingular terms at LO 3 , which are known analytically. These are not included in the Drell-Yan case, and so the difference between the analytic and Geneva distributions in figure 14 is more sensitive to the subleading corrections that the nonsingular terms generate.

JHEP09(2013)120
As in the e + e − case, the uncertainty bands for the resummed curves and Geneva predictions in figures 13 and 14 are obtained by adding the fixed-order and resummation uncertainties in quadrature. In the peak and transition regions of the distribution, the resummation uncertainties dominate, while the fixed-order uncertainty dominates as the resummation is being turned off in the tail region. Comparing the uncertainty of the resummed distributions with that of the LO 1 distribution, which is much smaller, one can see that the fixed-order uncertainty is an underestimate of the missing higher-order terms. The reason for this is twofold: the missing large logarithmic corrections at higher orders, whose associated uncertainties are instead included in the resummed distribution, and the partial cancellation between the renormalization and factorization scale dependence, whose variations are correlated in the results we show.
In both the e + e − and Drell-Yan processes, the partonic Geneva spectrum is determined by eq. (2.46), which for an event multiplies the fully exclusive fixed-order cross section by the ratio of the resummed cross section for the jet resolution variable divided by the fixed-order expansion of that resummation. Compared to the e + e − case, where each subprocess contributing to the cross section is trivially proportional, in Drell-Yan, the convolution with the PDFs requires treating every possible qq initial state separately, in both the fixed-order and the resummed cross sections. In Geneva, the flavor sum is performed in the Monte Carlo sense, since every event has a definite flavor for the initial-state quarks and the correct flavor-summed cross section is obtained after a sum over all events. This means that the separate factors in eq. (2.46) are evaluated for an individual flavor, and the entire expression is summed over flavors. In the analytic resummation, since the matching between the resummed and fixed-order cross sections is additive, there is instead only one way to perform the sum over flavors.
A version of the master formula where both the resummed and the resummed-expanded are separately flavor summed before entering eq. (4.4) would be equally valid. We have checked, however, that this is a very minor effect and is not the main contribution to the apparent differences of the Geneva predictions with the analytic resummed cross section. In fact, the reason for the discrepancy is in the difference of higher-order terms that are included in the Geneva multiplicative approach with respect to the additive matching used in the analytical calculation. This can also be seen as an indication of the relative freedom in implementing the master formula in eq. (2.46) at a given perturbative accuracy. As one can evince from figures 3(d)-3(f), should the resummation and fixed-order calculations be evaluated at the next order in perturbative accuracy, the size of the yet-missing terms would decrease, and consequently, the difference between Geneva and the analytic results would be reduced.
The Geneva implementation in this example can be extended to higher accuracy in terms of both fixed-order matrix elements and resummation. An equivalent accuracy to the e + e − results shown in section 3 can be achieved if the fixed-order matrix elements for the 1-jet multiplicity are calculated at NLO, the 2-jet multiplicity are calculated at LO, and the resummation of beam thrust is continued to NNLL . Although this is beyond the scope of this work, we nonetheless demonstrate that the Geneva framework is capable of merging matrix elements of different jet multiplicities beyond the lowest order. As in JHEP09(2013)120 the e + e − case, jet multiplicities are defined using a physical jet resolution variable, which allows for a consistent extension of the entire framework to O(α s ) perturbative accuracy.
In section 3, we found that NNLL resummation of the jet resolution variable T 2 , when combined with the parton shower, was capable of describing the spectrum in other 2-jet observables with an accuracy that clearly exceeded NLL, the naive accuracy of the parton shower. The resummation of T 2 captures an important set of logarithms that are correlated with other 2-jet observables, and, when combined with the fully exclusive, allorders description of the parton shower, the accuracy of other observables can be improved beyond NLL. At a hadron collider, the effective dynamic range of observables is much larger, meaning the correlation between the jet resolution variable and another observable of interest may be small. In this case, the parton shower may play a greater role in determining the spectrum, and hence the accuracy, of other observables. We will investigate these features in a future work.

Conclusions
In this paper, we have shown how to combine higher-order resummation of a jet resolution variable with fully differential next-to-leading-order calculations to extend the perturbative accuracy of cross sections beyond the lowest order for all values of the jet resolution scale. This framework has been interfaced with a parton shower and hadronization to give the Geneva Monte Carlo program.
Our framework provides both the versatility of fixed-order calculations and the accuracy of higher-order analytic resummation. From the point of view of Monte Carlo generators, the Geneva approach allows the combination of higher-order resummations with higher fixed-order calculations. From the point of view of resummed calculations, it allows one to obtain a fully differential cross section that correctly resums the jet resolution variable to higher logarithmic accuracy. Since this construction maintains the higher perturbative accuracy for all values of the jet resolution scale, it naturally allows NLO calculations of different jet multiplicities to be combined with one another.
The higher logarithmic resummation of the jet resolution scale allows us to use a low cut on the jet resolution scale, much lower than the point where fixed-order perturbation theory breaks down but still above the nonperturbative regime. This is a major difference to other approaches [28][29][30], where the jet-merging scale has to be chosen much larger, such that α s ln 2 τ cut 1. In this paper, we have concentrated on the theoretical construction, which is valid for any number of jets, and for both e + e − and hadron colliders. We have shown that one has to carefully choose a jet resolution variable that is resummable to higher logarithmic accuracy. In our approach, the N -jet and (N + 1)-jet regions are described by the same fully differential calculation without the need of an explicit jet-merging scale. The cut on the jet resolution variable is only needed due to the presence of IR divergences. We have given expressions for both the integrated cross section below the IR cutoff and the differential cross section above that properly combine the higher logarithmic resummation with a higher fixed-order calculation.

JHEP09(2013)120
This approach has been implemented in the Geneva Monte Carlo. As a first application, we have presented results for e + e − collisions. The jet resolution variable for this case was chosen to be 2-jettiness, which is directly related to thrust, and we combined its NNLL resummation with the fully differential 3-jet rate at NLO 3 . Varying the profile scales and the renormalization scale has allowed us to obtain event-by-event uncertainties. As a final step, we have interfaced our perturbative result with Pythia 8, which added a parton shower and hadronization to our results. The parton shower adds additional radiation beyond the highest jet multiplicity in Geneva and has been restricted to only fill out the jets of the exclusive jet cross sections at lower jet multiplicity. Since the cut on the jet resolution variable could be chosen to be very small, the effect of the perturbative shower (without hadronization) is rather small, and different tunes in Pythia do not affect the resulting distributions significantly. Hadronization has a significant effect, and the difference in final results due to different hadronization parameters is more manifest.
We have shown that Geneva correctly reproduces the higher-order resummation of the thrust spectrum, even after showering, which serves as a nontrivial validation of our approach. Using α s (m Z ) = 0.1135, as obtained in ref. [52] from fits to the thrust spectrum using higher-order resummation, together with tune 1 of Pythia 8, we obtain an excellent description of ALEPH and OPAL data. The same setup was then used to predict other event shape variables, namely C-parameter, heavy jet mass, and jet broadening. In all cases, our results agree remarkably well with the explicit analytical resummations, even though only the thrust resummation was used as an input. This comparison shows numerically that we achieve a higher resummation accuracy than NLL, which is what one would naively expect to obtain from the parton shower. This is especially remarkable for jet broadening, where the resummation formula has a completely different structure from the thrust resummation. Comparing our results after hadronization to the data, we again find excellent agreement for these other observables.
Finally, we have presented first results toward an implementation for hadron colliders in Geneva. Choosing the Drell-Yan process at the LHC with beam thrust as the jet resolution variable, we combined the resummation of beam thrust at NNLL with the leading-order matrix element for the emission of an extra jet. The results from Geneva agree well with analytical results, which shows the applicability of the framework to hadron colliders. As we have shown, our theoretical framework to combine higher-order resummation with fixed-order matrix elements and parton shower Monte Carlos is very general, and there are many avenues to pursue in the future. Obvious next steps for e + e − collisions are to include NLO calculations for 4 jets, which would require including the logarithmic resummation for 3-jettiness [90] as well as additional tree-level matrix elements. In addition, one can consider the resummation for other jet resolution variables. For hadronic collisions, next steps are to include the resummation and NLO calculations for higher jet multiplicities, as well as adding parton showering and hadronization using the different available programs.