Designing Gapped Soft Functions for Jet Production

Distributions in jet production often depend on a soft function, S, which describes hadronic radiation between the jets. Near kinematic thresholds S encodes nonperturbative information, while far from thresholds S can be computed with an operator product expansion (OPE). We design soft functions for jets that serve this dual purpose, reducing to the perturbative result in the OPE region and to a consistent model in the nonperturbative region. We use the MSbar scheme, and in both regions S displays the appropriate renormalization group scale dependence. We point out that viable soft function models should have a gap associated with the minimum hadronic energy deposit. This gap is connected to the leading O(Lambda_QCD) renormalon ambiguity in jet event shapes. By defining the gap in a suitable scheme we demonstrate that the leading renormalon can be eliminated. This improves the convergence of perturbative results, and also the stability by which non-perturbative parameters encode the underlying soft physics.

Soft functions play an important role in the study of cross sections close to kinematic thresholds, characterized by jets of collimated hadrons with small invariant mass. These cross sections are frequently described by factorization theorems involving hard Wilson coefficients, jet functions describing the jets of hadrons, and a soft function S. The hard coefficients and the jet functions are perturbative, while S encodes universal nonperturbative information on soft radiation between the jets. The prototype examples are event-shape distributions in e + e − annihilation for large c.m. energies Q [1,2,3], such as the thrust T [4,5,6], where T ≡ maxn i | p i ·n|/ i | p i | [7] and the kinematically allowed range is 1/2 < T < 1. In the threshold "dijet" region of large thrust, T ∼ 1, the events are characterized by two back-to-back jets, and at leading order in 1/Q the factorization theorem has two jet functions and one soft function [1,2,3]. Other examples include distributions for jet broadening [8], the heavy jet mass [9], and their generalization to angularities [10]. The dijet region also plays a crucial role in event shapes for massive particles, such as the invariant mass distribution of jets from top-quarks [11]. For applications at hadron colliders soft functions which account for initial state radiation are important [12]. Finally, for studies of weak B-meson decays to jets, soft functions involving the initial state B play a crucial role. Examples are B → X s γ and B → X u eν [13,14,15], as well as B → X s ℓ + ℓ − [16,17]. Here phase space cuts enhance the region where the soft function has a large effect.
Near threshold one can distinguish two regions. Very close to threshold the distribution typically shows an enhanced peaked structure, and nonperturbative information in the soft function is important for determining the shape and the maximum of the distribution. The size of this "peak region" is set by the hadronic scale Λ QCD . Next to the peak region the distribution typically falls off and shows a tail-behavior but is not yet highly suppressed. The dynamics is still dominated by jets and soft radiation, but in this "tail region" the leading soft function can be computed perturbatively since it is probed at scales larger than Λ QCD . In the tail region operators sensitive to nonperturbative physics are power-suppressed. Computations of moments involving integrations over both the peak and tail regions can be done with this same power expansion.
Since S encodes different types of physics in the peak and in the tail region, one possibility is to make separate predictions for the corresponding cross-sections. However, phenomenologically it is often desired to treat both regions coherently. In a pioneering analysis of e + e − → jets [18] this was handled by implementing a "hard" IR cutoff on the event shape variable "e". Perturbation theory was used above the cutoff and the perturbative corrections were frozen below it, with where R NLL PT contained perturbative results up to two-loop order with next-to-next-to-leading log resummation (NLL). The function R PT (e, Λ IR ) was then convoluted with a normalized soft function model S mod as dictated by the factorization theorem. With a simple choice for S mod good agreement with LEP data was found for several event shapes. This cutoff procedure does not attempt to treat explicitly the renormalization scale dependence in the region where the soft function is non-perturbative, nor does it systematically implement the perturbative corrections in this peak region.
The multi-region issue has also been analyzed in the context of B-meson decays. In Ref. [19] a perturbative tail was glued to the soft-function model, where S part is the "partonic" soft function obtained from perturbation theory and where the argument in the θ-function was chosen such that the tail turns on without discontinuity, using the condition S PT (ω, √ e(ω − Λ)) = 0. This method provides the correct renormalization group behavior for the treatment of the tail region at leading order, and is an improvement because it allows the perturbative jet function corrections to be incorporated systematically in the peak region. Shortfalls are that in the peak region it still hides the dependence on the renormalization scale µ in model parameters, and that the perturbative tail is turned on by hand at a particular point, rather than allowing it to appear once it dominates the non-perturbative corrections.
In this paper we develop a procedure for constructing soft function models for jets that i) reduce to the perturbative result in the OPE region and a consistent model in the nonperturbative region, ii) exhibit the proper renormalization group scale dependence in the MS scheme, iii) have a gap associated with the minimum hadronic energy deposit, and iv) are stabilized to perturbative corrections by being free from the leading O(Λ QCD ) renormalon. We show that the soft function gap parameter is essential for removing the renormalon ambiguity of the partonic threshold energy order-by-order in perturbation theory.
Although our procedure is quite general, in order to make all the steps explicit we will carry it out in the context of a specific example. We consider event shapes for top-quark jets produced in e + e − → tt at c.m. energies Q ≫ m t . The soft function we construct applies equally well for massless event shapes in the dijet region, that is, very little of our discussion depends on the presence of the top-quark mass or width. We consider the double differential top-antitop invariant mass distribution, d 2 σ/dM t dMt, where M t,t are either in the peak or the tail region. In the peak region near the top mass resonance, s t,t ≡ (M 2 t,t − m 2 t ) ∼ m t Γ t where Γ t is the top-quark width, and we have the factorization theorem [11] dσ peak which is valid at leading order in m t /Q and Γ t /m t . Here H is a calculable hard coefficient and B ± are calculable jet functions, whereas S np (ℓ ± , µ) is a nonperturbative soft function which peaks for ℓ ± ∼ Λ QCD when µ ∼ ℓ ± . In general the convolution probes momenta ℓ ± ∼ s t,t /Q in the soft function, and large logs in S are avoided by taking µ ∼ ℓ ± and summing large logs in the jet and hard functions. In the peak region s t,t ∼ QΛ QCD + m t Γ t , so the nonperturbative distribution described by S np (ℓ ± , µ) directly effects the differential cross section. On the other hand, in the tail region, s t,t ≫ QΛ QCD + m t Γ t , and the dominant momenta in the soft function are ℓ ± ∼ s t,t /Q. In the interesting region this is a perturbative scale of ℓ ± ≃ 3−30 GeV or larger, depending on the size of Q. The leading order factorization theorem in this tail region is which is valid to leading order in s t,t /Q 2 , m t /Q, and Λ QCD Q/s t,t . Here the partonic soft function S part (ℓ ± , µ) can be computed as a perturbative series in α s . Power corrections at O(Λ QCD Q/s t,t ) are determined from S np in a manner discussed below, while power corrections at O(s t,t /Q 2 ) involve new factorization theorems containing subleading soft functions (which have been worked out for inclusive B-decays [20,21,22,23]). The soft function carries information on how soft radiation is associated to the definition of the invariant mass variables M t,t . To be definite we consider hemisphere mass definitions where the soft function for both Eqs. (3) and (4) is [11] S(ℓ ± , µ) Here k +a s is the total plus-momentum of soft hadrons in X s that are in hemisphere-a, k −b s is the total minus momentum for soft hadrons in the other hemisphere. The soft function for thrust is related to the hemisphere soft function by with τ ≡ 1 − T , and we emphasize that S(ℓ + , ℓ − ) is independent of the top-mass. In general soft functions are matrix elements of Wilson lines, which in our case are In order to predict the invariant mass distribution in the peak and the tail regions we would like to connect Eqs. (3) and (4). In this paper we consider the task of constructing an appropriate soft-function that contains both S np and S part and which can be applied in the peak and the tail region. In order to be useful the result must remain consistent for scales µ ∼ s t,t /Q, both in the tail region where s t,t ≫ QΛ and in the peak region where s t,t ∼ m t Γ t + QΛ. We will consider all large logs to have already been summed by renormalization group evolution from Q down to these µ's. So the task is to determine the soft function matrix element at these µ's, where it should contain no large logs.
To begin, consider modeling the soft function by where S part (ℓ ± , µ) is the partonic soft function computed in perturbation theory, and S mod (l ± ) is a nonperturbative model function that is µ-independent and contributes only for ℓ ± ∼ Λ QCD . In Ref. [17] an analog to Eq. (8) was used in the study of b → sℓ + ℓ − to alleviate the issues mentioned about Eq. (2). Taking S part to O(α s ) this formula provided a simple way of incorporating the cutoff OPE moment constraints of Ref. [19] in the model for the nonperturbative B-meson soft function. Here we will argue that, suitably refined, Eq. (8) can be used to design soft functions for jets that are consistent with the desired properties stated earlier. Defining moments we will demand that S mod is normalized, S [0,0] mod = 1. We will also demand that higher moments are finite where we have S Since S [0,0] mod = 1 we have the desired result that S(ℓ ± , µ op ) = S part (ℓ ± , µ op ) at leading power. Computing the renormalized soft function in Eq. (5) to order α s ( Fig. 1 with no n f -bubbles) it factors as 1 1 We note that the factorized form of the soft function with respect to the two hemisphere light-cone variables ℓ ± in Eq. (11) allows for the possibility to choose two different µ's at which to stop running the two jet functions B ± in the factorization theorems (3) and (4). While we do not expect that relation (11) is maintained for non-logarithmic corrections beyond the one-loop level, one can prove that the factorized form is maintained to all orders as far the scale-dependence is concerned, as in Eq. (15) [24]. Thus it is possible to treat the situation where s t and st are widely separated and to account for the resulting non-global logarithms [25] by choosing both renormalization scales differently.
We see explicitly that large logs in S part (ℓ −l, µ) are minimized for µ ∼ ℓ −l. Hence when ℓ andl are parametrically different it is the larger of the two that is important for the proper setting of the renormalization scale in the soft function. This is compatible with the expansion in Eq. (10). In the convolution with the jet functions in the tail region in Eq. (4), the logs in S part are minimized for µ = µ op , and S part (ℓ ± , µ op ) can be determined by a truncated series in α s (µ op ). Thus for Eq. (4) the result in Eq. (8) works at any order in perturbation theory. We would also like S(ℓ ± , µ) to give a viable model for the peak region S np (ℓ ± , µ) in Eq. (3) when it is applied at a low scale µ = µ low > ∼ Λ QCD . Here ℓ ± ∼l ± in S part (ℓ ± −l ± , µ) for the convolution in Eq. (8). This convolution builds the proper µ-dependence into S(ℓ ± , µ), since the µ-dependence is determined by perturbation theory exactly as in S part (ℓ ± , µ). Thus it avoids the issue of having a µ-dependence related to the soft function anomalous dimension in the model parameters in S mod . The convolution with S part also generates a perturbative tail, implying that S(ℓ ± , µ) is not normalizable. To see this define the cutoff moments Using Eq. (8) with Eq. (12) and S up to terms of O(α 2 s ) or O(Λ QCD /L). Rather than a deficiency, this behavior of S L[0,0] is a necessary feature, as it is consistent with the renormalization equations for S(ℓ ± , µ). Only S mod needs to be normalized.
For the peak region, perturbative improvements to S part in Eq. (8) that cause a large change to S, could in principle be compensated by changes to the model parameters in S mod . However, it is quite desirable to make S part and S mod as independent as possible, so that the interpretation of the model parameters remains unchanged as we perturbatively improve S part . A measure for this independence is the convergence of the perturbative expansion for S part at µ low . In general the convolution in Eq. (8) generates double logarithmic terms, The choice of µ low should be small enough to avoid these potentially large logarithms, but large enough to ensure the validity of the perturbative expansion in α s (µ low ). Thus a satisfactory choice of µ low might be difficult to find, and requires careful examination. To test this issue we can determine the logarithmic series for S part (ℓ ± , µ), by finding the partonic soft function in renormalization group improved perturbation theory at LL order, NLL order, etc. The renormalization group improved S part satisfies the exact relation where U s is the LL, NLL, etc. evolution kernel. As indicated this kernel factors in the variables ℓ + and ℓ − to any order in perturbation theory [24]. Using this RG-improved S part (ℓ ± , µ) the full S(ℓ ± , µ) in Eq. (8) also satisfies the evolution equation (15) exactly, with a µ-independent S mod (l ± ). When the logs are small we can expand the RG-improved result to a fixed order in α s (µ), and the resulting S part and S satisfy the RG to this order. We will use this truncated version of the NLL series for S part (ℓ ± , µ) to test for a choice of µ which minimizes large logs in the soft function. This will also provide a test for the stability of model parameters to the addition of perturbative corrections.
Lets construct the NLL partonic soft function using a Fourier transform as in [26]. At NLL order the partonic soft functions factorize S N LL where the position space kernel is The LL results for ω and K involve Γ cusp 0 and β 0 and the NLL results involve Γ cusp 1 and β 1 , which also agrees with Ref. [27]. Here r = α s (µ)/α s (µ 0 ), C F = 4/3, β 0 = 11 − 2/3n f and To obtain a suitable boundary condition to solve Eq. (16) exactly, we note that the series of [θ(ℓ) ln k (ℓ/µ)/ℓ] + plus-functions in S part (ℓ, µ) become a series of ln k [i y µe γ E ] inS part (y, µ). Thus in position space we can take a boundary condition where all the logs are absent. For example, using the LO boundary conditionS part (y, µ = −ie −γ E /y) = 1 in Eq. (16) we It is straightforward to verify that this partonic soft function satisfies the evolution equation in Eq. (16). Specifying higher order boundary conditions forS part will then properly specify the subleading non-log terms in the series forS part . For instance,S part (y, µ = −ie −γ E /y) = 1 − πC F α s (−ie −γ E /yµ)/8 fixes the NLO boundary condition of Eq. (12). Thus the general solution to Eq. (15) is This result allows us to determine the LL and NLL series. Order by order in perturbation theory the Fourier transform (FT) can be carried out analytically since In addition to the leading logs, this inverse Fourier transform gives contributions to non-log terms from the expansion of e ǫγ E /Γ(1 − ǫ), which are subleading to the momentum space NLL series. As long as such subleading terms are unambiguously defined order by order and obey the RGE, one is free to include them in the NLL result. For our purposes we define the LL, NLL, etc. results as the resummed series obtained in position space, since it is in this space that the evolution equations are the simplest. With the NLO boundary condition and NLL evolution we find where L j = 1/µ θ(ℓ) ln j (ℓ/µ)/(ℓ/µ) + . Note that the coefficients for the terms beyond NLL order are incomplete, namely α 2 s L 0 , α 3 s L 2,1,0 , and α 2,3 s δ(ℓ). We show coefficients for these terms because of our convention of specifying the series in position space and using the full transform to momentum space. To obtain the complete α 2 s L 0 and α 3 s L 2 terms we would need to include the non-cusp part of the two-loop anomalous dimension.
Having determined the desired form of S part in Eq. (8) and a means to test for large logs, we now turn to the nonperturbative information in S mod and the overlap with perturbation theory. To satisfy the moment constraints on S [n,m] mod one can consider a two parameter model with exponential tails [18] where N (a, b) ensures f exp is normalized to one, and b = 0 controls the noninclusive correlation betweenl + andl − . Physically the range −1 < b < 0 is favored [18]. In the past this and other models used for soft functions in jet physics are taken to be nonzero for ℓ ± ≥ 0. This is a natural constraint given that it is satisfied to any order in perturbation theory for S part (l ± ). Withl ± ≥ 0, Eq. (8) enforces ℓ ± ≥ 0 in S(ℓ ± , µ). However, a better approximation is to take a soft-function with a gap so that the soft-function model vanishes forl ± < ∆,  Here ∆ encodes the minimum hadronic energy deposit in each hemisphere. 2 Since the model parameter ∆ ∼ Λ QCD it has an O(1) effect in the tail region where the soft function is nonperturbative. Among the model parameters ∆ plays a special role because it enables a hadronic interpretation for the variablesl ± ≥ ∆ in S mod (l ± ). Through the convolution in Eq. (8) this gap is transferred to give ℓ ± ≥ ∆ in S(ℓ ± , µ). This transfer relies on the fact that we have a partonic threshold at zero-momentum, i.e. that S part (ℓ ± −l ± ) has support only for ℓ ± ≥l ± . However, this transfer is not entirely straightforward because in perturbation theory the partonic threshold has a renormalon which yields an O(Λ QCD ) ambiguity in ∆. In the Borel transform of the hemisphere soft function considered here, this renormalon corresponds to a pole at u = 1/2. Since the soft function is universal for massless jets and top quark jets this renormalon is also behind the u = 1/2 Borel pole identified by Gardi [29] in an analysis of event-shape distributions in full QCD for massless partons. The nature of this soft function renormalon is similar to the well known O(Λ QCD ) renormalon of the heavy quark pole mass definition, but is not equivalent to it; rather it is specific to the soft function for jets. For example, the u = 1/2 renormalon pole of the soft function that occurs in inclusive B decays is solely related to the heavy quark pole mass, and is eliminated by switching to a short-distance threshold mass, see for example [19]. For the case of the top jet event shape distribution considered here the pole mass renormalon is contained in the jet functions [11], and is of no concern for the construction of the soft function. Only gluon fields appear in the matrix element defining our S in Eq. (5).
Using standard renormalon calculus either based on gluon propagators dressed with massless fermion bubble chains or on the modified gluon propagator the one-gluon exchange graphs in Fig. 1 give the Borel transform This parameterizes the leading O(Λ QCD ) renormalon ambiguity of the tree-level soft function, and has the same form as a shift in the zero point of δ(ℓ + )δ(ℓ − ) expanded to first order. It is also consistent with the result found by Gardi [29] for thrust, accounting for Eq. (6). Equation (25) can be generalized to soft function diagrams with an arbitrary number of gluons with one gluon modified by Eq. (24). Since one can use the soft limit for the modified gluon momentum (compared to the momenta of the unmodified gluons) only diagrams where the dressed gluon is external need to be considered. The computation of the contributions from the dressed gluon then factorizes from the remaining gluons yielding This result parameterizes the leading O(Λ QCD ) renormalon ambiguity of the soft function at any order. The Borel pole at u = 1/2 leads to instabilities in the perturbative predictions as we systematically include perturbative corrections to S part . As we will see below, such instabilities are for example reflected in S becoming negative in certain ranges of ℓ ± , or in an instability of the ℓ ± values where S is maximal. Physically, this ambiguity ties together the perturbative physics that we aimed to associate with S part and the hadronic information in S mod , and it must be resolved by experimental information.
In order to remove the ambiguity and allow for a stable determination from experimental data we would like to use a renormalon free scheme for the gap. Thus we take ∆ =∆ + δ where∆ is a renormalon-free model parameter for the hadronic threshold, and δ = δ 1 +δ 2 +. . . has a perturbative expansion which cancels the renormalon ambiguity in S part . Shifting variables tol ± =l ± − δ we have To cancel the renormalon ambiguity we must expand Eq. (27) in δ simultaneously with our expansion for S part = S 0 part + S 1 part + . . ., so that Here δ i ∼ O(α i s ) can be defined with any prescription that removes the O(Λ QCD ) renormalon ambiguity, and simultaneously this prescription will define a scheme for the hadronic param-eter∆. Note that ∆ is renormalization group invariant, thus∆ inherits a scale-dependence if δ is not renormalization group invariant. Moreover, we note that quadratic and higher powers of δ i that appear in Eq. (28) are required to ensure the consistency of the perturbative scheme. The terms linear in δ i are the ones relevant for removing the leading O(Λ QCD ) ambiguity, having the same form as Eq. (26).
In order to motivate a definition for a subtraction scheme associated to δ consider the first moment S L[1,0] from Eq. (13). For now the upper cutoff L is arbitrary. Starting from Eq. (27) we use the OPE as in Eq. (10) and expand to linear order in δ to obtain where in the second line we dropped α s corrections to the power correction, and here part which should be canceled by the δ-term in Eq. (29). A suitable form for δ to render the leading order factorization theorem and the first power correction renormalon free is .
(31) Note that different choices of L correspond to different schemes for renormalon-free gap parameters∆. Other ways to define δ are also feasible. From the expression for S part given in Eq. (21) we obtain Note that the one-loop δ 1 term is exact, while the two-loop term δ 2 relies on our NLL approximation of Eq. (21). Lets examine the impact of renormalon subtractions on the soft function. In  (14), so its shape and normalization change when varying µ.
In Fig. 2 the subtracted curves also show a somewhat smaller correction to the ℓ value where their maximum is located than the unsubtracted curves, but this effect is more dependent on the choice of parameters, such as the L value, see Fig. 3. At O(α s ) the perturbative series for the peak position has not yet approached its asymptotic behavior, but we expect the improvement in convergence for the peak position of the soft function to become more pronounced when higher order perturbative results for the soft function are considered.
To test whether S part suffers from large logs for particular values of µ, the O(α 2 s ) NLL predictions for the soft function from Eq. (21) are shown as the blue dot-dashed and dashed lines in Figs. 2 and 3. The dot-dashed curves do not have renormalon subtractions, and again exhibit negative dips. The dashed curve use our renormalon free∆, with subtractions given by the terms in the last set of square brackets in Eq. (28) and δ 1 and δ 2 from Eq. (32). We see that at this order the renormalon subtractions continue to eliminate the negative dip at small ℓ values. The behavior of the peak location for the two-loop NLL result is in general not dramatically improved, but this is simply because the O(α 2 s ) soft function given in Eq. (21) is based on a logarithmic approximation in a region where the logs are not large, and hence does not contain the large renormalon terms of the full two-loop soft function. Finally, for the lower right panel of Fig. 3, we see an indication for an instability The impact of the renormalon subtraction is also significant for the differential cross section. Let us first consider the peak region based on the factorization theorem (3). Since we only wish to illustrate the impact of the soft function, we use tree-level jet functions B ± (ŝ) = 1/(ŝ 2 + Γ 2 t ) for Q/m t = 5, Γ t = 1.43 GeV, m t = 172 GeV. We also ignore common normalization factors, and evolution factors that sum large logarithms down to the low renormalization scale µ low of the soft function, since they affect all predictions in the same way. (For the case of top-quark jets, a complete analysis including all these terms is carried out in Ref. [24].) Fig. 4 displays this differential cross section for equal invariant masses M = M t = Mt over M − m t for the three parameters sets of Fig. 2. Again we find that using a renormalon free gap parameter improves the convergence of the predictions and avoids the problem of negative dips in the cross-section. Interestingly, the curves show even better convergence compared to the soft function alone, and show nice convergence for the peak location. We find that this is true in general and related to the additional smearing that is provided by the width of the jet function. These results illustrate that the removal of the O(Λ QCD ) renormalon contributions in the soft function is essential to obtain a renormalon-free mass measurement from the peak position of the invariant mass distribution. We emphasize again that the renormalon issue in the soft function treated here is entirely independent of the pole mass renormalon problem, which appears in the massive jet function and the top quark pole mass.
Finally, let us examine the tail region of the differential cross section, using again tree-level   Fig. 5a are displayed, but now with the renormalon subtraction. Since the perturbative contributions in S part are at the scale µ op it is mandatory to choose L of order µ op to avoid large logarithmic terms, as can be also seen from Eq. (31), and we adopt the specific scheme choice L = µ op . Comparing the curves in Figs. 5a,b we see that the renormalon subtraction substantially improves the perturbative convergence. Figure 5d shows the difference between the solid/dashed and the dot-dashed curves from Fig. 5b.
Comparing it to Fig. 5c we see that the renormalon subtractions lead, as anticipated, to a significantly better perturbative behavior for values one would extract from the data for the power correction. 3 This illustrates that the renormalon subtracted predictions are essential for extracting stable and renormalon-free model parameters from experimental data. A scheme such as the one used here, where L = µ, works well for both the tail and peak regions, avoiding large logs. If a result for the gap model parameter is determined from data in a scheme where L = µ, then Eq. (32) can be used to relate the result to other schemes, such as for L = µ/2.
To conclude, we have provided a prescription for designing soft function models in jet production, that can be applied both in the peak region where the soft function is nonperturbative and in the tail region where the soft function can be expanded with an OPE. The method entails the convolution of the partonic soft function with a normalized model function that encodes the nonperturabive information, Eq. (8). It automatically implements consistent renormalization scaling behavior in the MS scheme, making the design particularly useful when dimensional regularization is employed for perturbative calculations. As a novel feature we argue that the soft function models need to exhibit a gap which accounts for the fact that for real hadrons there is a minimal hadronic energy. This gap is also required to devise a systematic scheme to remove the leading O(Λ QCD ) renormalon that is contained in the partonic soft function. In Eqs. (27,28,31) we have provided a simple definition for such a scheme and demonstrated that the removal of the renormalon avoids large uncertainties in predictions of the soft function and hence the cross-section in the peak region. In the tail region it also reduces the size of fluctuations in the power corrections, since they are otherwise affected by the O(Λ QCD ) renormalon. It is possible to generalize our method to treat also subleading O(Λ n QCD ) renormalons with n > 1, which are expected to have smaller effect on the soft function stability. Subtraction of these subleading renormalons might improve the numerical stability at higher order in perturbation theory of model parameters in S mod not related to the gap.