$C$-parameter hadronisation in the symmetric 3-jet limit and impact on $\alpha_s$ fit

Hadronisation corrections are crucial in extractions of the strong coupling constant ($\alpha_s$) from event-shape distributions at lepton colliders. Although their dynamics cannot be understood rigorously using perturbative methods, their dominant effect on physical observables can be estimated in singular configurations sensitive to the emission of soft radiation. The differential distributions of some event-shape variables, notably the $C$ parameter, feature two such singular points. We analytically compute the leading non-perturbative correction in the symmetric three-jet limit for the $C$ parameter, and find that it differs by more than a factor of two from the known result in the two-jet limit. We estimate the impact of this result on strong coupling extractions, considering a range of functions to interpolate the hadronisation correction in the region between the 2 and 3-jet limits. Fitting data from ALEPH and JADE, we find that most interpolation choices increase the extracted $\alpha_s$, with effects of up to $4\%$ relative to standard fits. This brings a new perspective on the long-standing discrepancy between certain event-shape $\alpha_s$ fits and the world average.


Introduction
The strong coupling constant α s is the least well known coupling in the gauge sector of the Standard Model. The latest Particle Data Group (PDG) average of α s has an uncertainty of about 1% [1,2], considerably larger than the error in the other gauge coupling determinations. Given the importance of QCD at LHC collider experiments, and rapid progress in perturbative calculations [3] and experimental accuracy, the uncertainty on α s is becoming increasingly critical for precision collider phenomenology. However, the headline figure of 1% uncertainty from the PDG average masks significant discrepancies between different extractions. In particular two of the most precise α s determinations come from event-shapes studies: α s = 0.1135 ± 0.0010 [4] from fitting thrust data and α s = 0.1123 ± 0.0015 [5] from C-parameter data. These results are several standard deviations away from the world average of 0.1179 ± 0.0010 [2,6] and from other individual precise extractions, such as 0.1185 ± 0.0008 from lattice step scaling [7] and 0.1188 ± 0.0013 from jet rates [8].
These particular event-shape and jet-rate fits are among the most precise of a wide variety of fits to e + e − hadronic final-state data [4,5,[8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23]. Many of them use a Current address: MTU Aero Engines AG Dachauer Str. 665 80995 Munich, Germany high-precision perturbative calculations, however they also all require input on non-perturbative (hadronisation) effects. These can be estimated either using Monte Carlo (MC) event generators [8,11,14,17,20] or via analytic non-perturbative models [4,5,13,18]. The use of MC event generators has long been criticised on two main grounds: they are tuned on less accurate perturbative (shower) calculations, and the separation between perturbative and non-perturbative components cannot easily be related to today's highest-accuracy perturbative calculations. Conversely the analytic models fit a non-perturbative parameter and the perturbative coupling in a single, consistent framework. The low α s values from Refs. [4,5,18] use the latter method. The price to pay in this approach is that the non-perturbative component is not controlled beyond the first order in an expansion in powers of 1/Q (the centre-of-mass energy) and furthermore only in the 2-jet limit, while fits cover both the 2 and 3-jet regions.
In this article we examine specifically the issue of going beyond the 2-jet limit for the hadronisation correction. In principle one might attempt a full calculation as carried out for top-quark production in the large-n f limit in Ref. [24]. Before embarking on such a calculation, we believe however that it is worth establishing whether 3jet hadronisation corrections bring a phenomenologically relevant effect. To do so simply, we consider the example of the C-parameter. This observable is special in that it has two singular points, one at C = 0 and the other at a Sudakov shoulder at C = 3/4 [25,26]. Existing fits calculate the hadronisation correction around the first singular point, C = 0 and extend it to the whole C-parameter spectrum. Here we point out that one can also quite straightforwardly calculate the power correction at the other singular point C = 3/4. One can then consider a range of schemes for interpolating between the two singular points and examine their impact on strong coupling fits.
This letter is structured as follows: in Section 2 we briefly recall the framework for perturbative fits with analytic hadronisation estimates. In Section 3 we then review the determination of the C-parameter hadronisation correction in the 2-jet limit and extend it to the symmetric 3-jet case (C = 3/4). Section 4 presents the results of our new fits and we then conclude in Section 5.

The C-parameter and its distribution
The C-parameter variable for a hadronic final state in e + e − annihilation is defined as follows [27], in terms of the eigenvalues λ i of the linearised momentum tensor Θ αβ [28,29], where | p i | is the modulus of the three momentum of particle i and p α i is its momentum component along spatial dimension α (α = 1, 2, 3). In events where all particles are massless, this can also be written as where Q is the centre-of-mass energy, p i denotes the fourmomentum of particle i, x i = 2(p i · Q)/Q 2 , and θ ij is the angle between particles i and j. We introduce the cumulative distribution Σ defined as The differential distribution dσ is known to next-to-nextto-leading order (NNLO) in massless QCD [30][31][32], which can be combined with the total cross section σ [33] to obtain a N 3 LO prediction for Eq. (4). The effects of heavyquark (notably the bottom quark) masses on event shape distributions [34], as well as electroweak corrections [35], are known to NLO, but we do not consider them in our study. Their omission does not affect in any way the conclusions of this article. In the following we consider the massless-QCD NNLO predictions for the differential distribution from Ref. [32].
In the two-jet region, the fixed-order expansion is spoiled by large logarithms of infrared and collinear origin, which must be consistently resummed at all perturbative orders to obtain a physical prediction. The resummation for the C-parameter distribution has been carried out in different formalisms [26,36,37] and it is known up to N 3 LL [36]. In our analysis we adopt the analytic next-to-next-leading logarithmic (NNLL) calculation from the appendix of Ref. [37], which is sufficient to illustrate our findings.
To obtain a perturbative prediction that is accurate across the whole physical spectrum, we need to match the resummed NNLL calculation to the fixed order result. This is done by combining the N 3 LO calculation Σ N 3 LO (C) with the resummed prediction Σ NNLL (C) according to the log-R scheme [38] as where Σ Exp. (C) is the fixed-order expansion of Σ NNLL (C) to O(α 3 s ). Detailed formulae are reported in Ref. [39]. Note that specific choices need to be made to limit the impact of the resummation in regions where C is not small. Our choices are discussed in Appendix D.
The hadronisation corrections to the C parameter distribution can be described in terms of an expansion in negative powers of the centre of mass energy Q. The leading correction in this sense leads to a shift of the cumulative distribution of the form where δC (C) ∝ 1/Q is the mean change in the C parameter's value due to the emission of soft non-perturbative radiation. In most work, the power correction is taken to be independent of the value of the observable, δC (C) ≡ δC (0). In this work we will be investigating the consequences of having δC (C) vary with C. We adopt the following form for δC (C), where ζ(C) encodes the C-dependence of the correction. In Eq. (7) we include the Milan factor M(n f = 3) 1.490 [40][41][42] to account for the non-inclusive correction. Eq. (7) involves the mean value of the strong coupling constant in the soft physical schemeα s (see also Appendix A) at infrared scales µ ≤ µ I , above which the prediction is assumed to be dominantly perturbative [43], Eq. (7) also includes terms to subtract the contributions already accounted for in the perturbative calculation [13,18,44]. The determination of the latter is not without subtleties, in that it assumes that non-inclusive corrections to such renormalon subtraction are described by the same multiplicative M factor as for the coefficient of α 0 (µ 2 I ). However these subtleties are numerically subdominant relative to other effects that we will be discussing here.
The core of this article relates to the coefficient ζ(C), which entirely determines the C dependence of Eq. (7). The calculation of ζ(C) at specific values of C is the subject of the next section.

Non-perturbative corrections
The calculation of δC near a singular configuration requires the amplitudes describing the emission of a "gluer" [43] k off the hard partonic system defined by the set of momenta {p i }. We can then calculate the mean change in the observable caused by the emission of k, i.e.

∆C({p
where the momenta {p i } ({p i }) describe the hard configuration before (after) the emission of the gluer. One can construct a range of prescriptions for mapping depends on the choice that is made. However, in the immediate vicinity of a singular configuration, it turns out that ∆C({p i }, {p i }; k)/k t is independent of that prescription in the limit of k t → 0, where k t is the transverse momentum of k. For the C-parameter, the singular points correspond to the 2-jet limit and the symmetric 3-jet limit. At these points the dependence on the choice of recoil scales as k 2 t and so vanishes in the k t → 0 limit of Under such conditions, one can write where [dk]M 2 (k) is the phase space and eikonal squared amplitude describing the emission of the gluer, and the {p i } are the hard momenta associated with the singular configuration. The coupling, colour factor and a dimensional factor Q are divided out, since these are included in Eq. (7). For an emission from a dipole {ij}, stretching between particles with momenta p i and p j , the matrix element and phase space are where k t , η and φ are to be understood with respect to the dipole; k t in the δ(k t − ) factor in Eq. (10) is also to be understood with respect to the emitting dipole.

Calculation of ζ(0) (2-jet limit)
Let us start by considering the two-jet limit C = 0. The shift in the C-parameter induced by a small-k t gluer is [26,38] where, for brevity, we have omitted the {p} and {p } arguments in ∆C(k). This leads us to where the rapidity limits can be taken to infinity because the integral is convergent. This coincides (to within conventions for normalisations) with the result that was given in [26,44,45].

Calculation of ζ(3/4) (Sudakov shoulder)
The leading order (LO) C-parameter distribution has an endpoint at 3/4. Just below this endpoint, the distribution tends to a non-zero constant [25], while above the endpoint the distribution is zero at order α s . This structure is known as a Sudakov shoulder. Parametrising the energies of the two quarks and the gluon as and the shoulder arises because of the absence of linear dependence on the 's. This absence of linear dependence on is also the reason that the choice of recoil prescription affects ∆C only at order k 2 t /Q 2 . Considering emission of a gluon with momentum k from an {ij} dipole (with i, j chosen among q,q, g), one can then derive in terms of the k t , η and φ of the emission with respect to the dipole (taken in the dipole's centre of mass). The O(k t /Q) contribution arises only when k is out of the 3-jet plane.
The squared matrix element times phase space can be written as a sum over three dipoles where C qg = Cq g = C A /2 and C qq = C F − C A /2. Note that for each dipole we will use the corresponding kinematic variables (k , etc.) in Eq. (16). This is equivalent to the procedure used to calculate the power correction to the D-parameter for arbitrary 3-jet configurations [46]. Integrating over η [dip] and φ [dip] , and summing over dipoles, we then obtain The functions K and E are the complete elliptic integrals of the first and second kind The numerical value of ζ(3/4) reads which provides the leading non-perturbative correction at the shoulder. 1 This simple result reveals that the leading (∼ 1/Q) hadronisation correction at the (symmetric threejet) Sudakov shoulder is less than half that in the two-jet limit (ζ(0) = 3π).

Modelling of the 0 < C < 3/4 region
Our calculations of ζ(0) and ζ(3/4) relied critically on the fact that recoil from the gluer emission had an impact that was quadratic in the gluer momentum. Away from these special points, the methods used here do not give us control over the value of the power correction, because the result depends on the prescription that we adopt for recoil (the impact of the hard parton's recoil becomes linear in the gluer momentum). One could conceivably extend the methods of Ref. [24] to attempt to determine the general dependence of ζ(C) on C, however such a calculation is highly non-trivial. So here, we want to establish whether such a calculation would be phenomenologically important. To do so, we consider a range of models that interpolate the power correction between the known values at C = 0 and C = 3/4, some of which depend on a parameter n ≥ 0. These are: where g(u) has the property that it is 0 (1) for u = 0 (1) and its first derivative is zero at u = 0, 1, The different forms for ζ(C) are shown in Fig. 1. The ζ 0 choice corresponds to using a constant shift, i.e. the standard approach for earlier studies. For both ζ a,n and ζ b,n , using n = 1 corresponds to a linear interpolation between the ζ(0) and ζ(3/4) values. For larger n, ζ a,n is flat close to C = 0, while ζ b,n is flat close to C = 3/4. Finally ζ c is flat near both C = 0 and C = 3/4. We stress that the variations in Eqs. (21) are not normally taken into account when estimating hadronisation with analytic models, which effectively all assume the ζ 0 model, corresponding to a constant shift across the whole differential distribution. In Section 4 we will see what impact this has on fits for the strong coupling from experimental data. In order to gain some insight on how ζ(C) depends on the recoil scheme, in Appendix B we carry out a fixedorder calculation of this quantity within different schemes to distribute the recoil due to the emission of the gluer among the remaining three partons. In reality, however, the behaviour that we find at fixed order in Appendix B can be substantially modified by the emission of multiple perturbative radiation (as also discussed in Appendix B). Therefore we do not rely on these calculations to assess the impact of ζ(C) on the fits, but rather use them as an insightful picture of how the leading non-perturbative correction scales across the spectrum of the event shape. We do however note that the concrete recoil schemes all yield shapes that fall below the ζ a,1 ≡ ζ b,1 line.

Fit of α s and hadronisation uncertainties
To test how our results affect the extraction of α s , we perform a simultaneous fit of the strong coupling and of the non-perturbative parameter α 0 (µ 2 I ), using data at different centre-of-mass energies from the ALEPH [48] and JADE [49] experiments, as summarised in Table 1 Ref. [5], but is largely sufficient for determining how the α s fit result depends on ζ(C).
The theory predictions are obtained using 50 bins in the 0 ≤ C ≤ 1 range, subsequently interpolated in order to be evaluated in correspondence to the experimental data bins. The fit is performed by minimising the χ 2 function defined as where V ij is the covariance matrix that encodes the correlation between the bins C i and C j . The general form of the covariance matrix is V ij = S ij + E ij , where S ij = δσ 2 stat, i δ ij is the diagonal matrix of the (uncorrelated) statistical errors in the experimental differential distribution, while E ij contains the experimental systematic covariances. The diagonal entries of E ii = δσ 2 syst,i are given by the experimental systematic uncertainty on the i-th bin. For the off-diagonal elements, which are not publicly available, a common choice (used also in Refs. [4,5,18]) is to consider a minimal-overlap model, which defines E ij as For ease of comparison, we adopt the same choice, though we note that for the normalised distributions that we fit here, the true covariance matrix would also include some degree of anti-correlation. The χ 2 minimisation is carried out with the TMinuit routine distributed with ROOT and the whole analysis was implemented in the C++ code used for a similar fit in Ref. [18]. Results with a diagonal covariance matrix, i.e. without any correlations, are given in Appendix C. They yield almost identical central results for α s and α 0 , smaller χ 2 values, and an increase in the experimental errors of O(10% − 20%), which however remain small compared to theoretical uncertainties. In order to estimate the theoretical uncertainties, we perform the following variations: • the renormalisation scale µ R is randomly varied in the range Q/2 ≤ µ R ≤ 2 Q, while the infrared scale µ I is set to 2 GeV; • for µ R = Q, the resummation scale fraction x C defined in Appendix D (default value x C = 1/2) is randomly varied by a factor 3/2 in either direction, namely in the range 1/3 ≤ x C ≤ 3/4, following the prescription of Ref. [9]; • for µ R = Q and x C = 1/2, the Milan factor M is randomly varied within 20% of its central value [40] (M 1.49) to account for non-inclusive effects in the δC shift (7) beyond O(α 2 s ); • keeping all of the above parameters at their central values, the parameter p in the modified logarithm defined in Eq. (41) of Appendix D (default value p = 6) is replaced by p = 5 and p = 7. This choice for p is discussed in Appendix D.
The theory error is defined as the envelope of all the above variations. When we quote overall results below, we add the theoretical and experimental errors in quadrature. We test several models for ζ(C) as given in Eq. (21) and shown in Fig. 1. Specifically, we consider the constant ζ 0 choice, the ζ a,n model for n = 1, 2, 3, the ζ b,n model for n = 1, 2, 3, and the ζ c model (recall ζ a,1 ≡ ζ b,1 ).
The results of the fits are given in Fig. 2 and Table 2. Fig. 2 shows results for α s and α 0 : the points give the central result for each ζ(C) choice, while the corresponding shaded areas represent the envelope of results obtained varying scales and parameters in the theoretical calculation, i.e. our overall theoretical uncertainty. Each point is accompanied by the ∆χ 2 = 1 ellipse, whose projection along each of the axes defines the 1 σ experimental uncertainty. Table 2 Table 2. Results of fits for αs and α0 using the different functional forms for ζ(C) reported in Eq. (21). The quoted uncertainties encode the total (statistical and systematic) experimental uncertainty (first number) and the total theoretical uncertainty (second number) estimated as described in the text. The The results with the ζ 0 model correspond to the standard implementation of the leading non-perturbative correction, which is assumed to amount to a constant shift across the whole C spectrum. The fit returns and agrees well with that of Ref. [5], albeit with larger uncertainties, in part due to our use of NNLL+NNLO rather than N 3 LL+NNLO theory predictions. We observe that several models lead to a χ 2 value that is the same as, or smaller than, that for the ζ 0 shape. In particular, the ζ b,2 model returns with a χ 2 that is similar to that of the ζ 0 fit. This corresponds to an increase in α s (M 2 Z ) of about 3.7%. In a number of models (ζ a,1 ≡ ζ b,1 , ζ b,2 , ζ b,3 , and ζ c ) the values of α s become compatible with the world average [2] α W.A. s = 0.1179 ± 0.0010. The result with the smallest χ 2 is the ζ a,2 model, which yields a rather small value of α s = 0.1121 +0.0024 −0.0016 . However the investigations of Appendix B, with a variety of concrete recoil-scheme prescriptions, seem to disfavour the ζ a,2 shape, suggesting that yet other factors may be relevant for maximising the fit quality.
Overall, the results suggest that one should allow for a 3−4% uncertainty in α s extractions from e + e − Cparameter data, associated with limitations in our current ability to estimate hadronisation corrections.

Conclusions
In this letter we have pointed out that the presence of a Sudakov shoulder in the differential distribution of some event-shape observables, such as the C parameter, can be exploited to gain insight on the observable dependence of the leading (∼ 1/Q) hadronisation correction to the spectrum. We found that the leading hadronisation correction at the Sudakov shoulder (C = 3/4) is over a factor of two smaller than the corresponding value in the two-jet (C = 0) limit.
In order to assess the impact of this observation on the fit of the strong coupling constant, we performed a set of fits using different assumptions on the scaling of the non-perturbative correction between the two points.
Our study is by no means exhaustive, and the inclusion of additional physical effects (such as the impact of bottom-mass effects) as well as a careful assessment of other sources of systematic uncertainty (such as the dependence on the fit range and the choice of correlation model) is necessary. However, it clearly reveals that current uncertainties in the modelling of hadronisation corrections can arguably impact the extractions of the strong coupling from event shapes at the several percent level. In particular, some of the models tested here lead to an increase in the extracted value of the strong coupling by 3% − 4%, which then becomes compatible with the world average to within uncertainties.
This necessarily raises the question of whether such observables should still be adopted for percent-accurate determinations of the strong coupling at LEP energies. Similar considerations may apply to extractions of α s obtained with jet observables, for instance those relying on accurate calculations for jet rates [31,[50][51][52][53] (e.g. the fits of Refs. [8,21]) or modifications of e + e − event shapes by means of grooming techniques [54][55][56] (an example being the analysis of Ref. [57]). Further studies are certainly warranted to investigate whether it is possible to better understand hadronisation for such observables across their whole spectrum, for example exploiting the large-n f calculational methods of Ref. [24].

A Some relevant quantities
In the present section we report the expressions for the anomalous dimensions used in the main text. The QCD β function is defined by the renormalisation group equation for the QCD coupling constant where the first two coefficients read The K (i) coefficients that appear in the non-perturbative shift (7) arise from the perturbative relation between the strong coupling in the soft physical scheme [58][59][60], denoted here byα s , and the MS coupling α s ≡ α s (µ 2 ) α s (µ 2 ) = α s 1 + α s 2π They read [58][59][60] B Fixed-order prediction for ζ(C) and recoil scheme dependence It is instructive to repeat the calculation (10) starting from a generic qqg configuration in the region 0 < C < 3/4. In this case the value of ∆C in Eq. (9) will depend on the specific scheme used to distribute the recoil due to the emission of the gluer among the remaining three partons.
Therefore, the definition of ζ(C) away from the singular points at C = 0 and C = 3/4 must be modified as follows with the normalisation N given by In the above two equations Φ qqg denotes the phase space of the qqg system and M 2 qqg ({p i }) is the corresponding squared amplitude evaluated with the unrecoiled momenta {p i } = {p q , pq, p g } prior to the emission of the gluer k. We see that as we approach one of the singular points ∆C({p i }, {p i }; k) ∆C(k) and we reproduce Eq. (10). Away from those points, the recoil will induce a linear dependence on the gluer's momentum, hence affecting the value of ζ(C) in a way that potentially depends on the specific model of recoil. The fixed-order calculation of Eq. (28) will provide some level of insight into how the leading non-perturbative correction varies across the spectrum of the event shape.
For an emission off a given dipole {ij} ({qg}, {qq} or {qg}), we express the gluer's momentum k by means of the Sudakov parametrisation where α k = (p j ·k)/(p i ·p j ), β k = (p i ·k)/(p i ·p j ), and k ⊥ = k t [n ⊥,1 cos φ + n ⊥,2 sin φ], with n 2 ⊥,m = −1, n ⊥,m ·p i/j = 0 (m = 1, 2), n ⊥,1 · n ⊥,2 = 0. We consider the following four recoil schemes 1. CS Dipole: the scheme is inspired by the Catani-Seymour map [61]. For an emission k off a dipole {ij} one identifies the emitter and spectator by considering the following quantity computed in the event centre-of-mass frame with = i, j. The emitter is then the dipole end corresponding to the smaller y k . Once the emitter (say p i ) and the spectator (p j ) are identified, the recoil is distributed as follows We also examined an alternative scheme in which the distance y k is computed in the dipole centre-of-mass frame. The two schemes produce identical results for the calculation considered in this appendix, and therefore we omit further discussion of the latter variant.

PanLocal [62] (antenna variant)
: the recoil is shared locally within the dipole ends as The quantities α i/j and β i/j in Eq. (33) are fully specified by the requirements (p ) 2 i/j = 0, (p i + p j + k) = (p i + p j ) and p i = p i for k t → 0. The PanLocal [62] rapidity-like variableη is defined as where s ij = 2p i ·p j , s i = 2p i ·Q, and Q is the total event momentum. In the event centre-of-mass frame,η = 0 corresponds to a direction equidistant in angle from p i and p j . The function f is responsible for sharing the transverse recoil among p i and p j and it is defined as Finally, we have with 3. PanGlobal [62]: the longitudinal recoil is assigned locally within the dipole as and the transverse recoil is assigned by applying a Lorentz boost and a rescaling to the full event so as to obtain final momenta {p , k } whose sum gives the original total momentum Q (see [62] for details).

FHP: the scheme is inspired by that proposed by
Forshaw-Holguin-Plätzer in Ref. [63]. It is similar to PanGlobal, with the difference that only the longitudinal recoil along the emitter, say p i , is assigned locallyp and the remaining longitudinal and transverse recoil is assigned by applying a Lorentz boost and a rescaling to the full event as in the PanGlobal scheme. Unlike the proposal in the original paper [63], we identify the emitter p i with the dipole end closer in angle to k in the event centre-of-mass frame, that is the one with the smaller y k defined in Eq. (31). For our purposes, this is physically similar to what is done in Ref. [63].
The results of the computation are reported in Figure 3, where for comparison we also report the curves corresponding to the profiles ζ a,1 ≡ ζ b,1 , ζ a,2 and ζ b,2 . We observe that the CS Dipole, PanLocal and PanGlobal schemes yield nearly identical results for ζ(C), which depart very sharply from the asymptotic value in the two-jet limit and approach the shape of the ζ b -type profiles at large values of C. Instead, the FHP scheme gives a less convex shape, close to a linear scaling in the fit region (indicated by the unshaded area in the plot). We believe that the similarity between the CS Dipole, PanLocal and PanGlobal schemes originates from the fact that, in the presence of a single perturbative gluon, ζ(C) is largely insensitive to the precise distribution of the transverse recoil among the particles in the event, which at this order and for this particular observable, is washed out by the integration over the azimuth and rapidity of the perturbative gluon p g . Conversely, the result does seem to depend on how the longitudinal recoil is assigned. The CS Dipole, PanLocal and PanGlobal schemes schemes all assign the longitudinal recoil locally within the emitting dipole, while in the FHP scheme part of the longitudinal recoil is shared among all particles in the event.
Note that all recoil schemes appear to be below the ζ a,1 ≡ ζ b,1 model. This tends to disfavour the ζ a,2 type model for interpolation of ζ(C) between C = 0 and C = 3/4, even though it gave the lowest χ 2 in the α s fits in Section 4.
A final comment concerns the limitations of the fixedorder nature of the study carried out in this appendix.
At the order at which we work, the fact that we force the perturbative qqg system to have a given value of the C-parameter causes the perturbative gluon to have a hardness comparable to CQ. Were we to go to higher orders, it would become possible for the perturbative event to contain additional, much softer gluons. Those gluons could also be involved in the recoil from the nonperturbative gluer, further altering the gluer's impact on the C-parameter. To take this into account, one would need to carry out a non-global type resummation that includes any number of perturbative soft gluons between the momentum scale set by the C-parameter value and the non-perturbative scale. First investigations in this direction (with just the PanGlobal recoil for the gluer) suggest that the impact of the resummation is non-negligible, and also continue to favour ζ(C) profiles that are below the linear ζ a,1 ≡ ζ b,1 profile.

C Fits with uncorrelated systematic experimental uncertainties
In this appendix we report the simultaneous fit of α s and α 0 obtained with the same procedure outlined in the main text, albeit replacing the model (23) with the simpler assumption of uncorrelated systematic uncertainties in the experimental data, namely The results are given in Table 3. Relative to Table 2, the absence of correlations leads to an increase in the experimental uncertainties for α s of O(10% − 20%), and the χ 2 values decrease. The central α s and α 0 results are essentially unchanged, as are the theory systematics.

D Modified logarithms and comparison to profile functions
In order to properly ensure that the resummation is turned off at the kinematic endpoint of the differential distribution, we modify the resummed logarithms by making the replacement where p denotes a positive parameter, and C max is the kinematic endpoint of the C-parameter distribution in the multi-jet regime, i.e. C max = 1. The prescription of Eq. (41) is but a possible choice and other sensible solutions can be found in the literature (see e.g. Ref. [5]). This ambiguity introduces an additional theoretical uncertainty in the calculation that must be carefully estimated. The quantity x C is of order one and its variation estimates the resummation uncertainty due to missing higher-logarithmic corrections. Specifically, the Σ N k LL (C) is defined in Ref. [5]. The centre-of-mass energy Q is set to the Z-boson mass.
resummed cross section acquires a net x C dependence such that, in the logarithmic limit as C → 0. Similarly, the parameter p determines how quickly the resummation is turned off in the region C ∼ C max . The choice of p must guarantee that the resummation does not substantially affect the prediction in regions of the spectrum dominated by hard radiation. An inspection of the first-order C-parameter distribution reveals contributions suppressed by a (linear) power of C relative to the dominant (ln C)/C dependence. Were we to take p = 1, the first-order expansion of the resummation would be associated with perturbative linear power-suppressed contributions whose coefficient would be larger than that observed in the exact fixed-order calculation. Accordingly, we believe it is sensible to apply the restriction p > 1 to avoid such contributions, and ensure that the resummation does not affect the dominant scaling at subleading power. With this constraint, we find that the extracted value of α s depends only very mildly on the choice of p and well within the quoted theoretical uncertainties. We then choose p = 6 and x C = 1/2 and vary both parameters as outlined in Section 4 in the uncertainty estimate. This specific choice is motivated by the fact that the scaling of the modified logarithm (41) in most of the fit range that we adopted in Section 4 happens to reproduce that of the profiled logarithms of the soft function of Ref. [5], which we use as a reference benchmark in our study. A comparison between the two prescriptions is shown in Figure 4 Model  Table 3. Results of fits for αs and α0 using the different functional forms for ζ(C), as in Table 2, but with an assumption of uncorrelated experimental uncertainties instead of the minimum overlap model.