Near-to-planar three-jet events at NNLL accuracy

We extend the ARES method for next-to-next-to-leading-logarithmic (NNLL) QCD resummations to three-jet event shapes in $e^+ e^-$ collisions in the near-to-planar limit. In particular, we define a NNLL radiator for three hard emitters, and discuss new features of NNLL corrections arising specifically in this case. As an example, we present predictions for the $D$-parameter, matched to exact NLO. After inclusion of hadronisation corrections in the dispersive approach, we compare our predictions with LEP1 data.


Introduction
Jet dynamics, encapsulated in event-shape distributions and jet-rates, is one of the most studied topics in QCD. These observables are designed to capture the continuous energy-momentum flow in hadronic processes, and as such offer a powerful probe of strong interactions. Jet observables probe disparate scales across the energy spectrum, starting from high scales where fixed-order perturbative calculations can be applied, and all the way down to Λ QCD where the yet unexplained phenomenon of hadronisation dominates the physics. Historically, and still up to this day, jet observables have been utilised to accurately extract the strong coupling from data as well as to test non-perturbative hadronisation models (see e.g. ref. [1] and references therein).
The study of jet observables has been pivotal in understanding the all-orders properties of QCD radiation, which subsequently lead to the discovery of non-global logarithms [2][3][4]. Distributions in jet observables can be computed at fixed-order in perturbative QCD [5], and such calculations have reached next-to-next-leading order (NNLO) accuracy for a number of relevant QCD processes. In particular, for e + e − annihilation, NNLO corrections to three-jet production have been computed in refs. [6][7][8][9].
Fixed-order calculations are reliable when the value of the jet observable is large. Nevertheless, the bulk of data lies in the region of small observable value where the cross section is dominated by large logarithms. The large logarithms emerge from the soft and/or collinear regions of phase space, due to a miss-cancellation between real and virtual corrections. Given a generic jet observable, the total (normalised) cross section is denoted by Σ(v), and represents the fraction of events where the observable takes a value less than v. In perturbation theory, Σ(v) displays logarithmic terms, L ≡ ln 1/v, and the highest power that appears at each order α n s of perturbation theory depends on the observable. For double-logarithmic observables, Σ(v) will contain logarithms as high as α n s L 2n . The fixed order approximation of the cross section becomes unreliable in the regime when α s L ∼ 1, and resummation becomes mandatory for theoretical consistency.
The primary concern of the resummation program is to reorganise the perturbative series in such a way as to allow those large logarithms to be isolated and resummed. Explicitly, the idea is to express ln Σ(v) as a series of functions with successive logarithmic accuracy. For double-logarithmic observables, we have ln Σ(v) = Lg 1 (α s L) + g 2 (α s L) + α s g 3 (α s L) + . . . , where Lg 1 (α s L) resums the so-called leading logarithmic (LL) terms, α n s L n+1 , g 2 (α s L) the NLL ones, α n s L n , g 3 (α s L) the NNLL ones, α n s L n−1 , and so on. Next-to-leading logarithmic (NLL) resummations have been available for many years for specific observables [10][11][12][13][14][15]. At the present day NLL resummation is available for all (continuously) global jet observables, which possess the property of recursive infrared and collinear (rIRC) safety [16][17][18]. The technique is based on a semi-numerical approach developed in a series of works, and the method is implemented in the computer program CAESAR [18], which automatically verifies whether or not a given observable is rIRC safe and continuously global. This paved the way to a systematic study of event shapes in hadronic di-jet production at NLL accuracy matched to next-to-leading order (NLO) results at hadron colliders [19,20].
There is no doubt that the theoretical predictions of NLL resummation are under a lot of tension due to various reasons. First, it remains true that NLL resummed predictions have a sizeable theoretical uncertainty which, when compared to current precision measurements, requires going beyond NLL. Second, recent works have started to utilise resummation results to test, and improve, the accuracy of parton shower simulations [21]. Given the absolute importance of parton showers for collider physics, it is mandatory that we push the accuracy of resummation results and aim for the widest class of observables. Third, event-shape distribution offer an important testing ground for analytical models of non-perturbative hadronisation corrections. Simultaneous fits of both the strong coupling and the parameter controlling the leading hadronisation corrections have been performed using NLL resummations for a variety of event-shapes (see [22,23] for the most accurate fits). Analogous studies using NNLL resummations exist only for a limited number of observables [24][25][26][27], it would be very interesting to have a new comprehensive picture of leading hadronisation corrections using NNLL resummations.
In the past decade, much progress has taken place in NNLL resummation. Nevertheless, most results available in the literature are performed for two-jet observables, i.e. those which vanish in the limit of two jets. Moreover, until very recently most NNLL resummations were observable specific, i.e. dependent on whether a factorisation theorem holds for the observable. Such approaches made it possible to obtain full next-to-next-to-leading logarithmic (NNLL) predictions for a number of global e + e − event shapes such as one minus the thrust 1 − T [25,28,29], heavy jet mass ρ H [30], jet broadenings B T and B W [31], C-parameter [26], energy-energy correlation (EEC) [32][33][34], heavy hemisphere groomed mass [35], and angularities [36,37].
Apart from e + e − annihilation, NNLL resummations have been performed for a number of jet observables in other QCD processes. For example, some results are available in deep inelastic scattering [38][39][40], and for hadronic collisions results were obtained when a colour singlet is produced at Born level. For instance, the transverse momentum of a colourless Boson in the final state [41,42], the variable φ * [43], the beam thrust [44,45], transverse thrust [46], and the leading jet's transverse momentum [31,[47][48][49]. We have a limited number of NNLL resummations for processes with more than two hard emitters, notably for heavy quark pair's transverse momentum [50,51] and N -jettiness [52,53]. Very recently, resummations for the Boson's transverse momenta and related observables has been pushed to N 3 LL accuracy [54,55].
Despite the advances of these approaches, many observables do not exhibit the factorisation requirements needed to carry out the resummation using, for example, the SCET framework [56]. This is particularly the case for observables which cannot be expressed as a simple analytic function of momenta, such as the thrust-major or the two-jet rate in the Durham algorithm. Very recently, the ARES (Automated Resummation for Event Shapes) approach has been completed and it is now possible to resum any rIRC safe di-jet observable, in e + e − annihilation, at NNLL accuracy. The original development of ARES focused on e + e − event shapes, but extensions thereof were presented for the two-jet rate [57]. ARES performs the resummation in direct space, i.e. without using any integral transforms, and only relies on the factorisation properties of QCD matrix elements in the soft and/or collinear limits. The fundamental ingredients of the method are as follows: • The analytic cancellation of soft and collinear divergences, which relies on the exact exponentiation of infrared poles in QCD processes with coloured particles in the final state. This exponentiation is simple to implement in the case of two, as well as three, hard legs.
• Unresolved emissions, owing to rIRC safety, yield a finite Sudakov radiator, which is analytically calculable in four dimensions. The radiator acts as a suppression factor for emissions contributing to the observable above δv, where δ defines a resolution scale.
• Resolved emissions do not generally exponentiate, but yield a class of functions starting at NLL accuracy. These functions are finite and directly calculable in four dimensions, thereby amenable to Monte Carlo integration. Each function has a distinct physical origin, and resums a given class of subleading logarithms that originate from various regions of phase space.
In this paper we extend the ARES method to resum three-jet observables, which are global and rIRC-safe, up to NNLL accuracy. Our work presents the first general NNLL resummation for three-jet observables, and paves the way to future extensions of ARES to go beyond three jets at NNLL accuracy. Given that the hard event is comprised of three partons, we are able to follow the basic constructions of the ARES approach outlined above. In particular, we derive the Sudakov radiator suitable for three-jet observables. Compared to previous two-jet results, the soft radiator in our case exhibits a richer structure in that it depends explicitly on the kinematics of the hard legs. The emission probability of soft partons, which are emitted coherently from the all the hard legs, takes the form of a sum over dipoles whose invariant masses end up appearing in the Sudakov radiator. The radiator also receives a contribution that is of pure collinear origin, and we report the full extension of this contribution in the presence of three hard legs. The dipole structure, in addition, leads to a new NNLL contribution arising due to resolved emissions. The new function is called ∆F wa and presents a new addition to the wide-angle NNLL function encountered in the ARES di-jet resummation formula. Moreover, we take full account of spin-correlations that become omnipresent due to the presence of a hard gluon in the Born configuration. The latter introduce a new ingredient in the resummation formula when the underlying gluon recoils against a hard emission. This leads to a new NNLL function, which we call ∆F rec , and we show how to compute it for a generic observable. The paper is organised as follows. In section 2 we define our notation and the various Sudakov decompositions we employ throughout the manuscript. We also give a lightening review of NLL resummation previously obtained in [58]. Section 3 contains the central piece of our work where we explicitly derive the NNLL resummation master formula. The Sudakov radiator is computed, before we move to formulate the various correction functions. We make sure to elaborate on the new ingredients that arise at NNLL, which manifest in two new correction functions. In section 4, we apply our resummation formula to the D-parameter, an example of an additive observable amenable to a fully analytic treatment in ARES. Finally, in section 5 we validate our analytic resummation against exact fixed-order calculations, and we present some simple phenomenological studies.

Kinematics and setup
In this section we set up the kinematics and notation that will be used throughout the paper. We explain the procedure we use to select three-jet events and give a short review of how to perform NLL resummation for three-jet observables in near-to-planar kinematics.
Three-jet Born kinematics. At Born level, a three-jet event in e + e − annihilation is made up of a quark of momentum p 1 , an antiquark p 2 and a gluon p 3 : with θ 12 and θ 13 the angles between p 1 and p 2 , and p 1 and p 3 respectively. The relation between these angles and energies and the variables that are typically used to describe three-jet events is reported in appendix B. We consider event shapes that vanish in the three-jet limit, i.e. V ({p}) = 0, where {p} denotes the set {p 1 ,p 2 ,p 3 }, the actual final-state momenta, which coincide with (p 1 , p 2 , p 3 ) in eq. (2.1) at Born level. After many soft and/or collinear emissions k 1 , . . . , k n , the three hard partons will recoil, and {p} are the actual final-state momenta after recoil.
Sudakov variables. A single emission k can be decomposed along any pair of light-like momenta (p i , p j ), which constitute the (ij) dipole, as follows: We have also introduced the invariant transverse momentum with respect to the (ij) dipole where (p i p j ) is a short-hand notation for the Lorentz-invariant product p i · p j . One can also choose to decompose the emission k along a single light-like momentum p . This can be achieved by defining the light-like momentump = (E , − p ). Explicitly, ⊥ is a two-dimensional space-like vector lying in the transverse plane to p , and whose magnitude reads (2.6) Note that, if k is collinear to p , we have κ (ij) → k We now introduce the rapidity η (ij) with respect to a dipole and its counterpart, η ( ) , with respect to leg p For an emission k collinear to p i or p j , the rapidities η (ij) , η (i) , η (j) are related as follows for a light-like vector k we obtain the rapidity bounds Last, we denote with φ ( ) the azimuthal angle of k ( ) ⊥ . We adopt the convention φ (ij) = φ ( ) = 0 when an emission is in the plane formed by p 1 , p 2 , p 3 . Comparing the expressions of the component of k outside the event plane in the Sudakov decompositions in eqs. (2.2) and (2.5), we have that, for an emission k collinear to Notice that the light-like momenta we use in eqs. (2.2) and (2.5) are not necessarily the Born momenta {p 1 , p 2 , p 3 } in eq. (2.1), neither the actual final state momenta, i.e. {p}. In fact, one can choose a different Sudakov decomposition according to the kinematical limit one is interested in, see for example [59]. The choice we make will depend on the context and will be made explicit in the subsequent derivations.
Lorentz-invariant phase space. The particular Sudakov variables we employ will depend on the form of singularities in the squared matrix elements. Hence, we list here the Lorentz-invariant phase-space measure expressed in various bases. For a massless emission k = (ω, k), we have: In the above equation, when we parameterise the phase-space in terms of leg-variables k ( ) t , η ( ) , φ ( ) , we omit the integration region corresponding to the anti-collinear direction p , because this does not correspond to any collinear singularity of QCD matrix elements. Also, the boundary of the region collinear to p is conventionally chosen to be η ( ) = 0.
Selection of three-jet events. In this paper we are interested in studying event shapes in the near-to-planar limit. In order to do this, we need a procedure to select hadronic events with at least three jets. This could be, for instance, through a jet algorithm that counts the number of well separated hard jets in the event, or through a cut on some secondary, two-jet observable. This constraint is represented by H(p 1 , . . . , p n ), a function of all hadron momenta p 1 , . . . , p n that is 1 if an event passes the cut and 0 otherwise. The function H(p 1 , . . . , p n ) also embodies exact energymomentum conservation. In our case, we use the Durham algorithm [60] and we select three-jet events if the three-jet resolution variable y 3 (p 1 , . . . , p n ) is greater than y cut . Correspondingly, we have a total three-jet cross section which, in d dimensions, is given by with dΦ n the n-particle phase space. We now consider the cumulative distribution of a three-jet event shape V (p 1 , . . . , p n ), defined as: (2.13) In near-to-planar kinematics, i.e. for v 1, the cumulative distribution Σ H (v) assumes the factorised form (see e.g. [18]) where p 1 , p 2 , p 3 are now the three Born momenta in eq. (2.1) and dσ 3 /dΦ 3 is given in eq. (B.3).
When v 1, the function Σ({p 1 , p 2 , p 3 }, v) develops large logarithms of v, which we want to resum, up to a given logarithmic accuracy, to all orders in the strong coupling. Notice that eq. (2.14) we have implicitly mapped the final state momenta, {p}, to the Born level momenta, {p}. The details of this mapping are not important here. In fact, with an IRC safe three-jet selection, in the presence of infinitely soft and/or collinear emissions k 1 , . . . , k n , the fineal-state momenta {p} always reduce to {p}, and the difference between H({p}, k 1 , . . . , k n ) and H(p 1 , p 2 , p 3 ) is suppressed by powers of v.
NLL resummation. We are interested in rIRC safe three-jet observables V ({p}, k 1 , . . . k n ) in the near-to-planar kinematics, in which V ({p}, k 1 , . . . k n ) 1. We require that, for a single soft emission collinear to leg , our observables behave as follows In the above equation, Q is a typical hard scale for the process under consideration, in our case by default the centre-of-mass-energy of the e + e − collision. Note that IRC safety requires a > 0 and b > −a. The NLL resummation of near-to-planar three-jet observables can be obtained from the general procedure of refs. [17,18]. At NLL accuracy, the distribution Σ({p 1 , p 2 , p 3 }, v) introduced in eq. (2.14) reads where R NLL (v) is the NLL radiator, encoding the probability of observing no emissions k i with V ({p}, k i ) > v, and R NLL (v) is obtained from the logarithmic derivative of R NLL (v) neglecting all NNLL corrections. The NLL radiator, obtained originally in ref. [18], will be borne out as a byproduct of our formalism in section 3. We recall here its expression as a sum of contributions of the three dipoles qq, qg, gq that build up a three-jet configuration: with L ≡ ln(1/v). In the above equation, C (ij) is the colour factor associated with dipole (ij), namely C (ij) = C A for (ij) = (qg), (gq) and C (ij) = 2C F −C A for (ij) = (qq). Then, we have the following integrals over selected momentum regions: where α phys s is the soft physical coupling defined in ref. [61]. To keep strict NLL accuracy, we calculate the integrals above neglecting all subleading contributions. The quantity ln(d g ) is an azimuthal average, which for any function f (φ) is defined as The terms proportional to the coefficient B represent virtual corrections of hard collinear origin down to the scale Qv 1 a+b , the characteristic scale of hard collinear radiation. The term proportional to T (L/a) represents soft wide-angle virtual corrections down to the scale Qv 1 a , the characteristic scale of soft wide-angle radiation. In fact, at NLL accuracy, when we perform the sum over dipoles, this is the only term that does not appear as a sum of contributions of each individual leg, but rather depends on the geometry of the underlying three-jet event. Introducing C , the colour factor of leg (C F for a quark and C A for a gluon), R ,NLL (v) ≡ 2C r (L), and R ,NLL (v) ≡ 2C r (L), we can recast R NLL (v) in the form: Here we have introduced γ (0) = 2C B , which is minus the coefficient of δ(1 − x) in the splitting function P gg (x) if p is a gluon, namely At NLL accuracy, the main role played by real emissions is that of cancelling the infrared singularities of virtual corrections. Only soft and collinear emissions give a non-trivial contribution, represented by the function F NLL (R NLL ) in eq. (2.16), where To define F NLL (R NLL ) we have to parameterise the momentum of each emission k i in terms of the leg i to which the emission is collinear, the azimuthal angle φ ( i) i , and the variable ζ i ≡ V ({p}, k)/v. For event shapes only, we do not have to specify the rapidity of emission k i , since it can be suitably integrated analytically. We then obtain the compact expression where we have used the short-hand notation (2.24) In the above expression, δ is a cutoff. All emissions with ζ i < δ are unresolved, and together with virtual corrections build the factor δ R NLL . Note that F NLL (R NLL ) has to be computed using V sc ({p}, k 1 , . . . , k n ), the approximate expression of V ({p}, k 1 , . . . , k n ) when all emissions are soft and collinear, and the quantity V sc ({p}, k 1 , . . . , k n )/v is only a function of In the ARES formalism, some formal manipulations have to be performed to obtain V sc ({p}, k 1 , . . . , k n ) from V ({p}, k 1 , . . . , k n ). The CAESAR program [18], which automatically resums all rIRC finalstate observables at NLL accuracy, instead generates actual momenta k i , and with appropriate phase-space cuts forces them to be soft and collinear. In that limit where V ({p}, k 1 , . . . , k n ) is the actual observable evaluated on the soft-collinear momenta k 1 , . . . , k n . The property of rIRC safety ensures that the limit in eq. (2.25) exists, and is only a function of Given the correspondence in eq. (2.25), we have V sc ({p}, k i ) = vζ i . A last remark is in order. When dealing with multiple soft and/or collinear partons, an emission might be collinear to one of the final state partons {p} without being collinear to the corresponding Born momentum that initiated the three-jet event. This needs caution in the parametrisation of the soft and/or collinear phase space. This issue has been discussed in detail in ref. [18] and recalled in ref. [61]. The outcome is that the light-like momenta required to perform a Sudakov decomposition according to eq. (2.5), which we will refer to as the "emitters", might need to be redefined after each emission. In general, the emitters do not coincide with the Born momenta in eq. (2.1), but are related to those via a mapping, whose details can be found in ref. [18]. Note that the emitters coincide with the Born momenta in the limit where all emissions are infinitely soft and/or collinear.

NNLL resummation
In this section, we explain how to extend the ARES formalism to three-jet observables. For the most part we rely on the previous results obtained in ref. [61], while making sure to stress the new features that arise in three-jet events. The quantity of interest is the cumulative cross-section which reads where dΦ3 +n is the phase spece for the final-state momenta {p} and the secondary emissions k 1 , . . . , k n . Here we remind the reader that the jet-selection function H implicitly contains a delta function for conservation of total four-momentum. In the near-to-planar limit, i.e. v 1, all emissions (k 1 , ..., k n ) are either soft and/or collinear. In this region of phase space it is always possible to provide an explicit mapping between the final state, (p 1 ,p 2 ,p 3 ), and the Born momenta (p 1 , p 2 , p 3 ) in eq. (2.1), see e.g. ref. [61]. Once the Born momenta have been identified, in the limit v → 0 the cumulative cross-section Σ H (v) assumes the factorised form [18] where S(n) is a symmetry factor, e.g. 1/n! for n identical gluons. In the above equation V({p}) represents virtual corrections to the Born process, an explicit function of the Born momenta {p} ≡ {p 1 , p 2 , p 3 }, normalised by Born matrix element squared. Moreover, M 2 represents the real corrections. Both M 2 and V ({p}) are divergent in four-dimensions, we therefore consider all our expressions to be regularised in some way. Specifically, we adopt dimensional regularisation, and all the quantities are computed in d = 4 − 2 dimensions.

Soft and collinear factorisation of matrix elements
A basic ingredient of ARES is the factorisation properties of QCD squared matrix elements, real and virtual. This makes the structure of divergences manifest and allows for the sought after cancellation of infrared poles. In this subsection, we recall the form of factorised amplitudes which are needed to build the NNLL resumed cumulative distribution.
Soft limit of squared matrix elements. The fundamental premise of ARES is the analytic cancellation of infrared singularities. To this aim, we start with the soft limit of real radiation, where we first notice the following factorisation 1 Since the final state contains only three coloured legs, it remains true that the soft limit of the squared matrix elements is factorised in terms of "soft correlated blocks", similar to the two-jet case (see e.g. [62]) As we shall explain below, eachM 2 s (k 1 , . . . , k n ) then takes the form of a sum over dipoles. The above decomposition in eq. (3.4) is a property of QCD. Nevertheless, one has to realise that beyond three jets an equivalent expression is very difficult to obtain and will involve complicated colour correlations [62].
Hard-collinear radiation. The second ingredient is radiation in the hard-collinear region of phase space. At a fixed logarithmic accuracy, we only need to consider a fixed number of hardcollinear emissions. The corresponding NNLL contributions are generated by a single hard-collinear parton k hc plus an ensemble of soft and collinear emissions k 1 , . . . , k n , all emitted independently. In this region of phase space, we have the following factorisation and the hard-collinear matrix element explicitly depends on the Born momenta through the mapping mentioned above. We will write where f = q, g is the flavour of leg . For a quark (or anti-quark) leg, we have where k t is the emission's transverse momentum with respect to its emitter, defined according to the procedure explained in ref. [18] and recalled in ref. [61]. We have the angular measure in the 2 − 2 -dimensional transverse plane The function is an azimuthally averaged splitting function where, to avoid double-counting the soft-collinear region, we have appropriately eliminated the divergent part of the full splitting function for z → 0.
For a gluon we have where the gluon averaged splitting function, with the soft divergence subtracted, is given by For collinear splittings of a gluon, we need to keep track of spin correlations with the hard event. These are accounted for by the un-averaged splitting function which has the property Spin correlations do not simply factorise from the Born amplitude, therefore, we have to introduce a new function T ({p}) of the Born momenta. In our case, we have [5] T where x 1 , x 2 are the invariants of the Born event defined in appendix B.

Virtual corrections. The last ingredient is the virtual corrections
where H({p}, α s (Q)) is a finite hard function, S({p}, α s (Q)) is a soft function [62] containing all soft singularities, while finally J (α s (Q)) encapsulates all hard-collinear singularities. We choose to incorporate a dependence on the Born momenta in the soft function, which we can always perform provided we appropriately adjust the hard function at each fixed order in the strong coupling.
For the sake of clarity, let us pause and further discuss the function J . The latter admits the expression [63] if leg is a quark or antiquark, and of P where β (d=4) is the beta function in four dimensions, given by the following expansion In anticipation of our next steps, we separate the collinear jet function as follows where the neglected terms, when combined with the corresponding real corrections, give rise to N 3 LL contributions. Last, the hard-collinear radiator R hc (v) can be defined at all logarithmic orders as follows 2 : Last we have the hard function, H({p}, α s (Q)), which is quite involved as it captures all the finite terms in the one-loop corrections to the Born event. Explicitly, in our case we have is given in ref. [5] and σ 0 is the Born quark-antiquark total cross section. Now that we have the various ingredients of eq. (3.2), we can define two cumulants, a soft and a hardcollinear cumulant, each encoding a separate non-overlapping portion of phase space.

Soft cumulative distribution
The soft cumulant Σ soft (v) is defined as which, as it stands, holds up to any logarithmic accuracy.
Exponentiation. The cancellation of the soft divergences is performed in two steps. First, we define a resolution variable, and divide the soft emissions into resolved and unresolved clusters according to the value of the resolution variable. The clustering algorithm is explained in ref. [61]. The resolution variable is a fake observable designed to cancel the soft divergences, with the only condition being that it has to share the same leading logs as the full observable. A natural choice is the soft-collinear limit of the observable in the presence of a single emission given in eq. (2.15). The second step is to express the soft function, S, in terms of the soft blocks in eq. (3.4) where the function W is called a web, whose properties will be discussed later on. Note that, the fact that we have written the soft function in term of an integral over real emission matrix elements implies definite kinematical boundaries for the k-integration. Note that our representation is a choice. Other representations, provided they correctly incorporate the soft singularities, lead to a different hard function H({p}, α s (Q)).
The unresolved clusters drop from the theta function in eq. (3.22) and thus the unresolved soft blocks exponentiate trivially 3 . Hence, the soft cumulant becomes where all soft divergences have been cancelled leaving a manifestly finite Sudakov radiator Note that in eq. (3.25), the soft radiator is still a function of the Born momenta and therefore still appears inside the integral over the Born phase space. To simplify notation we will drop the explicit dependence on the Born phase space in the remainder of the paper. In the above we have the resolution parameter δ, upon which the phase space is clustered. In particular, the phase space of the resolved clusters is bounded from below by making the corresponding resolution variable of each cluster bigger than δv. All expressions in eq. (3.25) are in four dimensions, and the dependence on δ cancels out in all of the final expressions.

Soft Radiator
It is in fact quite straightforward to compute the soft radiator in eq. (3.26) by directly utilising the two-jet results in ref. [61]. First, we realise that the soft blocks naturally take the form of a sum over dipoles. For example, the single-emission soft block at tree level reads where in our case This allows us to express the web function as a sum over dipole webs as follows (3.29) The notation k (ij) stresses the fact that each dipole web w({p}, α s (Q), k (ij) ) enjoys the property hence it is natural to express the momentum of the web in terms of the dipole variables in eq. (2.2). In order to compute the soft radiator we need to express also the resolution observable, V sc (k), in terms of same variables. This gives and d (ij) is defined in such a way that V (ij) sc (k) reduces to the expression in eq. (2.15) when k is collinear to leg . This gives the relation between d (ij) and the coefficient d introduced in eq. (2.15): This allows us to define a soft radiator for each dipole as follows in terms of which the total soft radiator becomes Notice that in eq. (3.34) the sole dependence on the Born kinematics is due to d (ij) . Furthermore, we realise that the phase space measure in eq. (3.34) contains the rapidity of a massive web momentum. Explicitly, where the rapidity is bounded as follows Since the web is uniform is rapidity, ref. [61] presented a simple strategy to extract the mass dependence of the web. Hence, each dipole contributes a radiator where the subscript '0' means that eq. (3.34) is to be evaluated with massless rapidity bounds, while δR (ij) (v) is a mass correction that accounts for the correct rapidity boundary of eq. (3.37).
At NNLL accuracy, this is given by To simplify the calculation of R (ij) ,0 (v) further we separate out the dependence on the Born momenta and on the azimuthal angle by expanding the step functions as follows where we truncated appropriately for NNLL resummation. Using eq. (3.40), we write and With this setup, it is straightforward to see that the dipole radiator in Eq. (3.42) is identical to the soft radiator of two-jet observables, which has been obtained in ref. [61]. The only difference is that now α s is a function of the invariant mass of each dipole, instead of the hard scale Q as in [61]. Nevertheless, we can re-expand the coupling to cast our results as a function of the resummation variable Implementing these steps we obtain and K (1) and K (2) are the coefficients of the soft physical coupling defined in ref. [61]. Finally, the mass correction reads , (3.53) which, up to NNLL accuracy, has no dependence on the dipole kinematics.

Hard-collinear radiator
In ARES, the Sudakov radiator receives a contribution from the virtual hard-collinear region truncated at the collinear scale 4 which is completely independent of the dipole structure, i.e. it only knows about the emitting leg, and hence the dependence on b . Following ref. [61], we write The various functions are expressed in terms of the coefficient of δ(1−x) in the regularised Altarelli-Parisi as follows This concludes the analytic construction of the Sudakov radiator.

Soft resolved clusters and correction functions
The contribution of resolved clusters to the cumulant, the second line of eq. (3.22), arranges itself in the form of a correction function, where contributions with successive logarithmic accuracy can be systematically extracted. This correction function is (3.62) The above expression starts at NLL accuracy. The next-to-leading logarithms emerge only from soft and collinear emissions widely separated in angle, which build the function F NLL . The remaining contributions emerge from relaxing, one at a time, the approximations performed to obtain F NLL . This systematic procedure gives rise to various NNLL functions, which we analyse in the following. To avoid clutter of notation, we will drop the dependence on {p} in the soft radiator R s .
The first thing to notice is that the observable in eq. (3.62) is not constrained to be evaluated in the soft limit. In fact, one of the gluons k 1 , . . . , k n can be hard and collinear, although emitted with the soft matrix element. This introduces an overlap between the soft and hard collinear regions, which we need to isolate. Let us denote this hard-collinear gluon as k hc . With this gluon, we have Therefore, subtracting the double counting with the soft-collinear region, we obtain the following NNLL contribution (3.64) This contribution will be incorporated in the function δF rec , to be discussed later on in section 3.3 along with the other NNLL functions of hard-collinear origin. Given this logic, another contribution naturally arises when one of the emissions is soft, but emitted at large angle. This gives the following correction where, in eq. (3.65), V wa means that we need to probe the observable in the limit when a single soft gluon k is emitted at large angles. In the second step function, the observable is evaluated as if k were soft and collinear. Without this subtraction, F wa would contain NLL terms that are part of the function F NLL defined in eq. (2.23). Notice that we dropped the resolution parameter, δ, from the integral over the extra emission, k, since the difference of theta functions renders the result finite in the limit δ → 0. Notice importantly that F wa , as it stands, contain contributions beyond NNLL accuracy. We will show below how to isolate the NNLL contributions. When this is done we get F wa (α s /π)δF wa (λ). The explicit expression of δF wa (λ) will be discussed later. After extracting F wa from eq. (3.62), we observe that we still have logarithms of arbitrary accuracy. At NLL, it suffices to treat all soft emissions as independent, while starting at NNLL we need to take into account the correlated portion of the double-soft squared matrix element, as follows 5 : The above form is valid up to NNLL accuracy. We now aim to re-arrange it for manifestly finite integrals that produce exact logarithmic accuracy. We concentrate first on the part that contains the correlated matrix elementM 2 s (k a , k b ). The dependence on δ in the integration over (k a , k b ) can be eliminated in the most elegant way by replacing the strong coupling for each soft emissioñ M 2 s (k i ), eq. (3.27), by the physical coupling that defines the soft radiator [61] . (3.67) With this replacement, we can isolate the function F correl (v) that starts at NNLL accuracy, and is given by In eq. (3.68), m 2 is the invariant mass of the correlated pair, i.e. m 2 ≡ (k a + k b ) 2 . The second step function precisely encapsulates the inclusive limit of the double emission that allowed us to introduce the physical coupling in the soft-collinear matrix element squared. Notice in particular that we dropped δ from the double emission phase space because, once again, the integral is manifestly finite in the limit δ → 0.
The previous steps leave us with the following expression The last step is to extract F NLL from the second line of the above equation, and isolate the remaining NNLL soft contributions. First, we expand the exponential prefactor in eq. (3.69), up to NNLL, as follows where ∆R s,NNLL denotes NNLL contributions to the first derivative of the full soft radiator, R s (v). Now, using eq. (3.67), for each leg in a certain dipole we implement the following transformation of variables In particular, ξ ( ) represents the rapidity fraction of the emission, i.e. the ratio of the emission's rapidity with respect to leg to the maximum available rapidity at fixed observable value vζ. The essence of the above transformation is that the observables we are interested in are event shapes and therefore do not depend on the rapidity fraction ξ ( ) . This allows us to integrate out this variable for each emission, and we reconstruct the logarithmic derivative of the massless radiator. In terms of these new variable, for each emission, we can write Now it is straightforward to extract various contributions to F s (v) by expanding (R ,0 ) (v) as follows: ζ .
(3.73) We can then write (3.74) In the above equation, F NLL (λ) is given by (3.75) To avoid confusion, we introduced the notation i in the above equation to denote the legs in a fixed dipole to which the i th soft emission belongs. This expression can be further simplified by observing that, when an emission k is soft and collinear to leg , the azimuthal angle φ (ij) reduces to the azimuthal angle with respect to the leg, i.e. φ ( ) . Then, we can also eliminate the sum over dipoles by defining a leg-dependent quantity We first eliminate as much as possible the dependence on the cutoff δ. The procedure, introduced in ref. [64], consists in writing logarithms of δ as integrals over an auxiliary variable ζ. Using we obtain (3.80) We can further simplify the above expression by extracting a piece that is purely soft and collinear, and that can be seen as the generalisation of δF sc introduced for two-jet observables [64]. However, we anticipate that, in the current case, ∆F s contains a term that is manifestly of wide-angle origin, and therefore is more naturally associated with F wa . First, similar to eq. (3.76) we can define another leg-dependent function ,0,NNLL ) (v) .
Using the explicit expression of the full radiator given in section 3.2.1, we have, to NNLL accuracy, We stress that, at NNLL accuracy, (R (ij) ,0,NNLL ) does not depend on the dipole kinematics, but only on the leg contained in the dipole (ij). Combining all terms that depend on Q ij we obtain, to NNLL accuracy which corresponds clearly to a term of soft wide-angle origin. Last, we define ∆F s = α s (Q) π (δF sc (λ) + ∆F wa (λ)) .

(3.87)
We now briefly derive the form of each NNLL correction that is suitable for numerical integration.
Soft-collinear NNLL correction. We collect from eq. (3.82) all the terms that depend explicitly on each leg, and not on the event geometry, and we make use of the new function ∆R ,NNLL defined in eq. (3.85). This gives the generalisation of the soft-collinear function δF sc introduce for the two-jet case in ref. [64]: Soft wide-angle NNLL correction. Let us move to eq. (3.65) and extract the NNLL contribution. Since the emission k is at largest angle with respect to all the others, to this aim all soft emissions are independent, hence where, once again, the single-emission soft block is defined with the physical coupling. For soft and collinear emissions, we can introduce the soft-collinear measure dZ[R i,NLL , {k i }] following the same steps as for F NLL . Furthermore, for the soft wide-angle emission k, we perform a change of variables that reflects the dependence of F wa on the dipole kinematics. We then use the Sudakov variables of eq. (2.2), and for each dipole (ij) we further introduce (3.90) Following the same steps as in ref. [64], at NNLL accuracy, we obtain F wa (v) = (α s (Q)/π) δF wa (λ), where (3.91) We then collect from eq. (3.82) all terms that contain the ratios Q ij /Q. This gives the new NNLL function

(3.92)
Soft correlated NNLL correction. Eq. (3.68) can be simplified further to extract the NNLL contributions. First, we write the correlated portion of the double-emission tree-level matrix element in terms of the variables introduced in appendix A: where and S, H g , H q can be found in appendix A. Note that the variables κ (ij) , η (ij) , φ (ij) refer to the Sudakov decomposition of the parent momentum k = k a + k b , and the construction is explained in Appendix A. Now in eq. (3.68) we change variables in a similar fashion to eq. (3.71) Owing to the fact that the observable does not depend on ξ ( ) , we can integrate it out analytically and find (3.97) Using eq. (3.81), and the fact that k a , k b are soft and collinear to the same leg , we can approximate φ (ij) φ ( ) , and finally obtain F correl (v) = (α s (Q))/π)δF correl (λ), where (3.98) The master NNLL formula for the soft cumulative distribution. We now put together all the ingredients to write down our master formula for the soft cumulant, valid up to NNLL accuracy × F NLL (λ) + α s π (δF sc (λ) + δF wa (λ) + δF correl (λ) + ∆F wa (λ)) . (3.99)

Hard-collinear cumulative distribution
Up to NNLL, the hard-collinear cumulant reads Our first task is to cancel the collinear divergence in the above expression. To this aim we notice that the singularity is encoded solely in the portion of the hard-collinear matrix element proportional to the Born amplitude, namely the pieces having the averaged splitting functions in eqs. (3.7) and (3.10). We will show below that the extra piece in eq. (3.10), proportional to the un-averaged splitting function, produces a finite term due to the vanishing of the azimuthal average, eq. (3.13), as k t → 0. We can partition the above expression into various pieces in order to arrange for manifestly finite expressions that could then be evaluated in 4 dimensions. The steps follow ref. [64], albeit with a new contribution arising from the spin-correlations of the gluons. We have comprises a constant leftover after cancelling the collinear divergence. We can easily evaluate eq. (3.102) using the splitting functions given in eqs. (3.7) and (3.10). For a (anti)-quark we get 6 This result coincides with that of ref. [61] for two hard legs. For a gluon we have Moreover, we have two correction functions. The first arises solely due to our choice of regularisation in eqs. (3.102) and (3.108), and reads We can further simplify the above expressions by introducing the phase-space measure over softcollinear emissions. We first introduce the observable fraction of the hard emission and using eq. (3.71) for the soft-collinear emissions we find F hc (v) = (α s /π)δF hc (λ)

(3.107)
Notice that in the above expression we can send the upper limit of the ζ integral to infinity, with corrections suppressed by powers of v. The second correction incorporates the recoil of the event shape due to the hard-collinear emission where the extra emission, k, is treated as soft and collinear in the second step function. This correction can be conveniently combined with the overlap function F s/hc introduced in eq. (3.64).
Performing the same formal manipulations that lead to eq. (3.107) (see also ref. [64] for details) we obtain that, at NNLL accuracy, F rec + F s/hc = (α s /π)δF rec , where where C q = Cq = C F and C g = C A . Note that, in eqs. (3.105) and (3.109), we obtain identical results if we use the following alternative definition for ζ In fact, what really matters is only the scaling of the hard-collinear emission k with respect to k t . It is crucial that we pause here to address an intricate point in the derivation of eq. (3.109). The squared matrix elements in eqs. (3.7) and (3.10) are expressed in terms of the transverse momentum with respect to the emitter. Consequently, we have to express the observable in terms of the same set of variables utilised in the matrix elements, which might turn out to be non-trivial depending on the observable. Any specific event shape will either use an axis in its definition, e.g. the thrust axis, or will depend on the relative transverse momentum between the particles of the final state. For soft emissions, the situation is simple because the direction of the emitter is the same as the direction of the final state hard momenta, up to terms that vanish as k 2 t → 0. For a hard-collinear emission, extra care must be taken. The transverse momentum appearing in the observable defines the integration variable ζ, and its precise relation to k t in the emission probability, i.e. eqs. (3.7) and (3.10), must be explicitly worked out.
Finally we have a new correction which is absent in the case of di-jet observables, and is due to the spin-correlation in the final state. Explicitly, we have (3.111) The above integral is indeed finite, but requires extra care. The collinear divergence is regulated because as k t → 0, the azimuthal average of ∆P g (z, φ; ) vanishes identically in d = 4 − 2 . The finite contribution that arises can be isolated. Let us change variables according to eq. (3.110) and extract the NNLL contribution 7 Now we can use (see for example [65]) where the plus distribution is defined as follows Hence, as promised the pole term disappears because the azimuthal average vanishes in 4 − 2 dimensions, cf. eq. (3.13). Applying the plus prescription yields a finite result Although the second step function vanishes because of the azimuthal average, it is still quite important to keep it in order for numerical integration to be feasible. The goal is to utilise the second step function as a regulator in a Monte Carlo integration. To simplify the implementation one ideally wants to push the limit of the ζ integral to infinity. For most observables, the first integral is effectively cut off by the observables constraint, so we can push the limit of the ζ integration to infinity. There are however a number of observables, especially those who are affected by cancellations of the contribution of emissions with comparable values of ζ, for which the integral is damped by the result of the integration over the soft-collinear measure dZ[R i,NLL , {k i }] (see e.g. appendix H ref. [18]). For those observables, the integral in ζ converges for R NLL lower that a certain critical value, which is the region in which our resummation is valid. Note that this consideration applies to F NLL and to all NNLL corrections, and we recall it here for completeness. We can also split the second integral at ζ = 1, and the contribution from 1 to 1/v vanishes identically upon azimthal integration. Hence, one can identically recast the above expression as follows: Eq. (3.116) is one of the main results of this paper, and is suitably defined for numerical evaluation.

Additive observables
For additive observables, such as the D-parameter, we can obtain closed form expressions for all NNLL functions. Additivity implies that the observable can be written as the sum of contributions of individual emissions. For soft and collinear emissions k 1 , . . . , k n , this means The evaluation of eq. (3.75) becomes simple and yields the well-known result (see e.g. [18]) , (3.118) where R NLL ≡ R s,NLL . Using eq. (3.117) we can compute δF sc , ∆F wa and δF correl using the procedure described in appendix C of ref. [64], and we get Proceeding to the wide-angle correction, all what we really need is to probe the observable, when a single soft emission k is emitted at wide angle. If we parametrise k using the Sudakov decomposition in eq. (2.2), and for an additive observable, we obtain where ζ is defined in eq. (3.90), and Using the above relations in eq. (3.91), we find . (3.124) Here, it is important to notice that the wide-angle correction is sensitive to the invariant mass of the dipole in contrast to the soft-collinear correction in eq. (3.119). For δF correl , we follow the procedure of ref. [61]. In particular, for k a , k b collinear to leg , we can write where ζ is defined in eq. (3.96). After some formal manipulations, we obtain Now we discuss NNLL contributions induced by hard-collinear radiation. For additive observables, following appendix C of [64], for the hard-collinear correction δF hc , we find where γ (0) arises due to the integral over the splitting function. Now we move to computing the function δF rec . First, we write where ζ is now given by eq. (3.110), and thus we have Now we follow almost identical steps to the treatment of the soft wide-angle correction and we get .

(3.131)
Finally, we have the new function ∆F rec . Instead of starting from eq. (3.116), we show here that this function can be computed directly using eq. (3.112) and use additivity, as follows: Owing to rIRC safety, we can rescale the ζ i 's, cf. ref. [64], to construct F NLL Now this integral is well defined in 4 − 2 dimensions, and therefore we can rescale ζ → ζ/f where now the ζ integral can be trivially performed and yields (3.135) Finally, we recall eq. (3.13) and expand the above equation around = 0 to find our final expression

NNLL resummation of the D-parameter in the near-to-planar limits
As a proof of concept, in this article we concentrate on a specific three-jet event shape, the Dparameter. This is defined in terms of the determinant of the spherocity tensor [5] where the sum runs over all hadron momenta p i and Q is the centre-of-mass energy of e + e − annihilation. The spherocity tensor has three eigenvalues λ 1 , λ 2 , λ 3 satisfying λ 1 + λ 2 + λ 3 = Tr Θ = 1. Using these eigenvalues we construct the C-parameter and the D-parameter D = 27 det Θ = 27λ 1 λ 2 λ 3 . For an isotropic event all eigenvalues are equal to 1/3, and hence both the C-and the D-parameter are equal to 1. Another useful form of the D-parameter is [66] which is very convenient to obtain analytic expressions for the D-parameter in the soft and collinear limits, as needed to compute the various components of our resummation master formula. In particular, in the presence of multiple soft emissions k 1 , . . . , k n , eq. (4.4) can be approximated as follows: where k i = (ω i , k). Note that, in the presence of soft emissions, the final-state hard momentã p 1 ,p 2 ,p 3 can be approximated by their Born counterparts p 1 , p 2 , p 3 . Therefore where k ix is the component of k i in the direction of p j × p k , i.e. out of the plane formed by the Born momenta p 1 , p 2 , p 3 . Using the fact that, for three particles (see e.g. [66]) we obtain the final expression for the D-parameter in the presence of soft emissions: where λ 1 λ 2 has to be computed using Born momenta.
NLL resummation. To compute the NLL resummation of the D-parameter we consider its behaviour after a single soft emission, collinear to leg . Using eq. (4.8) and the Sudakov parametrisation in eq. (2.5) we obtain Comparing the above expression with eq. (2.15) we get: This information is enough to compute the resummed cumulant at NLL. We first note that the D-parameter is additive, i.e. obeys eq. (3.117), which is clear from eq. (4.8). The parameters in eq. (4.10) allows us to directly compute the NLL radiator using the following relation and plugging a = b = 1 in eq. (2.17). Finally, F NLL is given by eq. (3.118), where one merely computes the logarithmic derivative R NLL .
NNLL resummation. In order to use our prescription for the NNLL radiator, we first construct d (ij) for each dipole by combining eq. (4.10) with eq. (3.33). Once we have d (ij) , we can compute the soft NNLL radiator using the formulae of section 3.2.1. In particular, we utilise the following relation The hard-collinear coefficients C hc, can be computed by replacing ln(d g ) in eqs. (3.103) and (3.104) with the appropriate expression in eq. (4.11).
We now consider the various real-emission NNLL corrections. The function δF sc is the one for additive observables, and is given by eq. (3.119). Furthermore, due to additivity both the wide angle, ∆F wa , and the hard-collinear, δF hc , functions are given by eqs. (3.120) and (3.127).
To compute recoil corrections, δF rec and ∆F rec , we need to obtain the expression for the Dparameter after a single hard splitting of leg . This produces an emission k with a fraction z ( ) of the energy E of the parent momentum p , and a final-state momentump carrying the remaining energy fraction 1−z ( ) . The important point to notice is that both k andp carry equal and opposite out of plane momenta,p ,x = −k x . From eq. (4.4), labelling the remaining two hard partons with the indexes 1 , 2 = , we have to consider four terms (4.13) If we add an arbitrary number of soft and collinear emissions k 1 , . . . , k n , their transverse momenta are much smaller than that of the hard collinear emission, which is the only one that effectively recoils against the hard partonp . In particular, the soft emissions do not change the direction of the emitter, up to non-singular corrections. Therefore, k x is the out-of-event-plane component of the transverse momentum with respect to the emitter p . Therefore, k x coincides with the emission's transverse momentum with respect to the emitter p , and we get This means that the D-parameter is additive also in the presence of an extra hard and collinear emission. For z → 0 we have Using eqs. (4.15) and (4.16), as well as the additivity of the D-parameter, we can compute δF rec using eq. (3.109) as follows: -30 -Using eq. (4.15) we can also compute ∆F rec from eq. (3.133), as follows: The next NNLL correction we need to compute is δF wa . According to eq. (4.8), for soft emissions the D-parameter is additive, so we can make use of the general expression in eq. (3.124). To achieve this we need to recast the expression of the D-parameter, with a single soft wide-angle emission, in the form of eq. (3.123). Using the Sudakov decomposition in eq. (2.2), for the dipole (ij) we obtain where . (4.20) Using eq. (3.122) and eq. (4.10), we find (4.21) Inserting the above expressions in eq. (3.124) we obtain where In fig. 1 we provide a plot of the integral in eq. (4.23) as a function of the three-parton variables (x 1 , x 2 ) defined in Appendix B. The plot shows the explicit result only for the qq dipole, and we choose not to explicitly show the similar plots for either the qg or theqg dipoles. The only difference is simply that the contours rotate in the (x 1 , x 2 ) plane. The last contribution we need to compute is δF correl . Since the D-parameter is additive, we can again use the general formula for additive observables in eq. (3.126). We recast the D-parameter, with two soft-collinear emissions, in the form of eq. (3.125). This gives which is the same for all three legs. This gives 25) and the various integrals are computed via a Monte Carlo routine H q ln f correl (µ, z, φ, φ ( ) ) = 1.1562 .

Validation and phenomenology
In this section we validate the analytic results of section 4, match them to fixed order at NLO and finally compare our matched distribution to LEP1 data. We do so in two steps. First, we use the Monte Carlo event generator EVENT2 to check most, but not all, of the pieces in the expansion of the resummation. For the D-parameter, EVENT2 provides results at LO, i.e. O(α 2 s ), and thus we will validate all terms in the expansion at O(α 2 s ) up to NNLL. Second, we use NLOJet++ to match the resummation at NLO, i.e. O(α 3 s ). Below we explain our choice of the matching scheme and point out interesting aspects of the resulting phenomenology.

Partial validation using EVENT2
In this subsection we start with EVENT2 to validate various ingredients in the resummed cumulative distribution of section 4. Given that the Born event is already at O(α s ), all the pieces in the expansion of the resummation that starts at O(α s ) can then be checked against EVENT2. Table 1 lists these various terms, which contribute at successive logarithmic accuracy. Moreover, the results of EVENT2 are sufficient to validate the geometry dependence in the radiator at NLL, i.e. the d -dependent term in eq. (3.42).
In fig. 2 we expand the resummation of the cumulative cross section and subtract the result from EVENT2. We see indeed that the difference is consistent with zero, and displays an asymptotic behaviour for the entire resummation region.
δF wa (λ) δF rec (λ) ∆F rec (λ)  Moreover, we can isolate and validate an extra NNLL function using EVENT2 results, i.e. δF correl (λ). Indeed, we can not achieve this directly because the fixed-order expansion of δF correl (λ) starts at O(α 2 s ). Nevertheless, we can look for an event shape whose behaviour in the presence of correlated emissions is identical to that of the D parameter. We shall call this observable D 2−jet , and it reads wheren beam is a unit vector along the electron beam direction. The design of D 2−jet is meant to capture the behaviour of δF correl (λ) for the actual D-parameter. This observable does not possess a lot of the structure witnessed in three-jet observables, nevertheless, it allows us to check our Monte Carlo implementation for obtaining the results in eq. (4.26). We have computed the analytic resummation of D 2−jet , although we do not quote the results here. Fig. 3 shows the result of subtracting the resummed differential distribution from that of EVENT2. Admittedly, the plot in fig. 3 does not exhibit a satisfactory asymptotic behaviour, however, it is quite suggestive. By making use of quadruple precision one should be able to probe sufficiently small values of D 2−jet . Finally, we have also tried to fully validate our resummation of the D parameter, at O(α 3 s ), using NLOJet++. Unfortunately, double precision did not allow us to reach values of ln 1/D larger than 10−12, which are not asymptotic enough to provide a reliable validation of our NNLL resummation.

Matching to fixed order at NLO
In order to provide a suitable cumulative distribution that paves the way for phenomenological studies, one has to match the resummation to fixed order. Matching is required to provide results across all values of the observable. The basic idea is to combine the results of both the resummation and fixed order, while making sure to get rid of contributions that are double counted. There are two generic conditions that the matching procedure must satisfy. First, based on physical grounds the matched total cross section should go to zero as v → 0. Second, the matched distribution, or likewise the total cross section, must reproduce the fixed order at the kinematic endpoint v → v max .
The two most popular matching schemes for e + e − annihilation are the R and log-R schemes [12,13]. In other contexts, multiplicative matching schemes are used [67,68]. Adopting one over the other is a choice that depends on the problem at hand. In our case, we could not achieve a stable matched distribution using the R and log-R schemes, given that the available data sets for the D parameter forces us to use low values of y cut . The problem with both schemes is that, given that the various components of the resummed cross section contain powers of ln y cut , the resummation does not switch off quickly enough and ends up substantially contributing to the tail of the matched distribution. This situation might be expected given that the K-factor NLO/LO is very large, approximately 100 %. This is similar to the case of resuming the distribution in the Higgs transverse momentum p t,H where the K factor is also known to be large [54]. We therefore need to supplement our matching scheme with a factor that effectively damps the resummation at large values of D.
Based on the above discussion, we use the multiplicative matching scheme designed in ref. [54]. The goal of that scheme is precisely to suppress the large terms, present in the resummation, which emerge outside the resummation region. This enables us to control the tail of the distribution and achieve a stable matching. In this scheme, matching is performed on the level of the total cross section. Given that NLOJet++ simulates the inclusive cross section, i.e. integrated over Born kinematics with the three-jet selection cut, we have to match on the same level. Explicitly, we have where H denotes the expansion of the resummation to NLO. The details of the matching can be found in appendix C and the expanded version of eq. (5.2) is given in full in eq. (C.5). The presence of the step function in eq. (5.3) suggests that the transition region, between the resummation and fixed order, will not be smooth enough. In fact, we verified that even if the step function is removed from the definition of Z, the resummation still shuts down smoothly well before reaching the kinematic endpoint v max . We carry out the matching using the values u = 1, h = 3 and v 0 = 1/2. As is customary in resummed calculations, we need to probe the size of subleading logarithmic terms. This is done using two simultaneous variations. The first introduces a rescaling x V as follows In the above, X is a variable choice to define the resummation scale, i.e. the logarithms being practically resummed, while X V controls the scale variation. We expand the total cross section around ln x V /v neglecting subleading terms. Furthermore, the resummed logarithm, ln x V /v, must be modified in order to impose that the total cross section is reproduced at the kinematical endpoint v max [68] ln where p denotes a positive number that controls how quickly the logarithms are switched off close to the endpoint. The parameter p is free, but is only constrained by the behaviour of the fixed order distribution near the endpoint [68]. In our case, we set p = 1.
Another estimate for the uncertainty in our matched distribution comes from varying the renormalisation scale, µ R , around a central scale that we take to be the centre-of-mass energy of the hard scattering, Q. For LEP1 energies, Q = M Z corresponding to α s (M Z ) = 0.118 while for FCCee energies, Q = 500 GeV corresponding to α s (500 GeV) = 0.094.
We implement two different choices for X in eq. (5.4). The first is referred to as the X const scheme which corresponds to setting X = 1 in eq. (5.4), while the second is the X prod scheme which corresponds to setting which is a function of Born kinematics. Finally, we construct the uncertainty bands by varying µ R by a factor of two in either direction and X V by a factor of three-halves in either direction. In figs. 4 and 5 we plot the matched distribution for Q = M Z using the two resummation schemes and for two different values of y cut , namely y cut = 0.1 and y cut = 0.05. We immediately notice the following features: • The uncertainty bands are not drastically reduced when increasing the logarithmic accuracy of the resummation, at least when compared to the typical situation with two-jet observables.
• The position of the peak is stable under varying y cut .
• For NLL, the uncertainty bands remain almost unchanged with decreasing y cut . In contrast, the uncertainty bands for NNLL are noticeably enhanced as we increase y cut .
The above features can be traced to the fact that jet selection generates terms that go as ln 2 y cut , for each power of α s relative to the Born cross section. The largest transverse momentum of softcollinear emissions, at fixed value of D, is of the order of √ DQ. Our resummation is strictly defined when the largest momentum is much smaller than the largest transverse momentum available, the latter being of the order of √ y cut Q. Essentially, our resummation is formally correct, as D 1, but phenomenologically viable only in the limit D y cut 1. Inspection of figs. 4 and 5 shows that the most probable value of D, which corresponds to the position of the peak of differential distributions, is of the same order as y cut and that is why we see the features described above. This is also reflected in the sensitivity of the uncertainty bands of the NNLL distribution to the variation of y cut in comparison to NLL. Simply, the NNLL pieces in the cross section, e.g. ∆F wa , contain extra powers of ln y cut compared to NLL. These logarithms are large, for y cut = 0.05-0.1, and thus we observe this behaviour of the uncertainty bands. The situation becomes better at FCC-ee energies, as we see clearly in figs. 6 and 7. Noticeably the position of the peak tends towards smaller values of D, and we start approaching the strict resummation regime D y cut 1. Simultaneously we see a reduction in the uncertainty by almost 50%. To conclude, for this observable, and depending on the value of y cut , we expect large subleading corrections that are not under control in any resummation formalism. This calls for a joint resummation of both types of logarithms, the observable and y cut , along the line of the presented resummation of both p t,H and the transverse momentum of the leading jet [69].
Leaving these caveats aside, we note that NNLL corrections generically yield harder D-parameter distributions. The effect is larger using the X prod scheme, because the resummed logarithms in the latter scheme are typically larger than the X const scheme. Indeed, this is one of the reasons why the NNLL uncertainty bands get larger when we use the X prod scheme, while their counterparts at NLL remain virtually the same. Note that in the X prod scheme the resummation scale is effectively of the order √ y cut Q which is the appropriate upper bound for transverse momenta. Therefore, this scheme automatically captures some of the terms which are enhanced by logarithms of y cut . Last, we compare our predictions to existing LEP1 data [70]. In order to do so, we need to supplement our perturbative resummation with some estimate of non-perturbative hadronisation corrections. Before we do this, we need to choose whether to use X const or X prod as our default choice for the resummation scale. We have observed that NNLL distributions obtained with X prod are not   Figure 6. The matched distribution for ycut = 0.1 and Q = 500GeV. The left plot using the Xconst scheme and the right using the X prod scheme.  Figure 7. The matched distribution for ycut = 0.05 and Q = 500GeV. The left plot using the Xconst scheme and the right using the X prod scheme.
very stable with respect to the choice of the matching parameter v 0 , which points to the fact that such a choice brings in numerically large subleading corrections, which we cannot control within our framework. Therefore, we decide to present non-perturbative plots using X const as our resummation scale, and v 0 = 1/2. We have checked that using other values of v 0 does not change considerably our findings. We include hadronisation corrections in the dispersive approach of ref. [71], where leading hadronisation corrections result in a shift of the corresponding perturbative distributions. In our case, we use the non-perturbative shift computed in ref. [72], and define where In the above equation, the geometry dependent functions g ij are the ones of ref. [72], which we rewrite using our own notation and conventions as follows: The non-perturbative parameter a NP is given by where α 0 (µ I ) = and α s (k) is the dispersive coupling defined in ref. [71]. As in previous non-perturbative studies, we set µ I = 2 GeV. In eq. (5.7), we have also introduced the factor (5.12) that ensures that the shift vanishes at the endpoint of the distribution. Also, to ensure that the distribution vanishes at its endpoint, we replaceL defined in eq. (5.5) with  Figure 8. The matched distribution, including the non-perturbative corrections, is compared to data from LEP1 for the two values of ycut we adopt in this article.
plots for non-perturbative matched distributions, with central scales, corresponding to NLL and NNLL accuracy. The non-perturbative shift corresponds to a value of α 0 (2 GeV) that give us a good agreement with data, and is within the range favoured by existing fits to event-shape data [23]. We see that NNLL resummation has a shape that resembles data more closely than NLL resummation, whilst favouring a similar value of α 0 , namely α 0 = 0.5. This trend persists irrespective of the value of y cut . This value of α 0 is similar to the central value of a fit obtained with the NNLL thrust distribution [24]. Finally, increasing the value of v 0 up to D max does not change the distributions close to the peak, but gives a better agreement with data in the tails.

Conclusions
This article presents a general method to compute the NNLL resummation of rIRC safe observables for processes characterised by the presence of three hard emitters. The method is a generalisation of the ARES approach to NNLL resummations, and paves the way to a general NNLL resummation with an arbitrary number of hard emitters. Although we concentrate on three-jet events in e + e − annihilation, our treatment of NNLL contributions induced by final-state radiation is completely general. Similar to the two-jet case, we are able to combine unresolved real radiation and virtual corrections to the Born process into an analytically computable NNLL radiator. The remaining corrections are all induced by real radiation, and can be computed for a general observable using suitable Monte-Carlo procedures.
Two new functions appear in the three-jet case. First, a new NNLL correction of soft wideangle origin appears, due to the fact that now we have three-hard emitters with non-trivial colour correlations. Second, since we have a hard gluon initiating a three-jet event, we need to take into account non-trivial spin correlations in hard-collinear splittings. These are embedded in a new NNLL correction that adds to those of hard-collinear origin.
As an example, we have applied our method to the D-parameter. Since this is an additive observable, we are able to compute most NNLL functions analytically, with a couple of integrals to be computed numerically. Then, we have performed phenomenological studies by matching our resummation to exact fixed-order and presenting predictions for LEP1 and future colliders. Both validation of resummation and phenomenology is tricky for three-jet observables. First, while it is possible to easily check NLL contributions against exact fixed-order, it is impossible to check NNLL ones without resorting to quadruple precision. With the aid of a fake two-jet observable that resembles the D-parameter, we have been able to check some NNLL contributions using the NLO code EVENT2. For what concerns the actual phenomenology, current cuts to select three-jet events give rise to large subleading effects at LEP1 energies. The situation is a bit better at FCC-ee. Nevertheless, we envisage that, to improve phenomenological studies of the D-parameter, one should attempt a joint resummation of logarithms of the D-parameter and of the variable determining the three-jet selection, with a similar procedure to that for angularites, or for the transverse momentum of a colour singlet and an accompanying leading jet.
A comparison with LEP1 data requires the inclusion of non-perturbative hadronisation corrections. We have added to our NNLL resummation the leading hadronisation corrections evaluated in the dispersive model. In general, the NNLL distribution has a shape that is similar to data, so hadronisation corrections are compatible with a rigid shift of perturbative distributions. We find that the value of the non-perturbative parameter α 0 determining the size of hadronisation corrections is comparable to the one obtained with NLL distributions. It might be very interesting at this stage to perform a comprehensive simultaneous fit of α s and α 0 using NNLL resummations for different event shapes.
In conclusion, our study sets the main building blocks for a general NNLL resummation of rIRC safe final-state observables with an arbitrary number of hard emitting legs. The only missing ingredient is a general treatment of both initial-state radiation and soft wide-angle corrections for a system with more than three hard emitting legs. Despite the technical difficulties, the philosophy of our method stays unchanged. In particular, ARES does not depend on the specific factorisation properties of an observable, and gives promise to achieve a fully general solution to the problem of NNLL resummation in the near future.
This makes it possible to write the angles between pairs of momenta in terms of the x i 's as follows The three-parton cross section, differential in x 1 and x 2 , in four dimensions reads dσ dx 1 dx 2 = σ 0 C F α s 2π with σ 0 the Born cross section for producing a quark-antiquark pair in e + e − annihilation.
To obtain the Born three-jet cross section σ 0 (y cut ) with the Durham algorithm [60] we need to integrate the differential cross section in eq. (B.3) with the constraint y 3 (p 1 , p 2 , p 3 ) = min{y 12 , y 13 , y 23 } > y cut , where y ij is the "distance" between pairs of partons defined by The Durham algorithm defines a six-sided region in the (x 1 , x 2 ) plane, as shown in Fig. 9 for the three different values of y cut we consider here. The corresponding Born cross section σ This cross section can be computed analytically. Its expression, not particularly illuminating, can be found in [73].

C Full matching formulae
In our matching formulae Σ Mat. H (v) we normalise all of the distributions to the total cross section σ H . However this is not what is provided by NLOjet++, instead it provides the un-normalised differential distribution for the D-parameter. We can transform the output of NLOjet++ into our conventions as follows. First we compute the un-normalised, barred, total cross section Exp. (v) . (C.5)