Determination of perturbative QCD coupling from ALEPH $\tau$ decay data using pinched Borel-Laplace and Finite Energy Sum Rules

We present a determination of the perturbative QCD (pQCD) coupling using the V+A channel ALEPH $\tau$-decay data. The determination involves the double-pinched Borel-Laplace Sum Rules and Finite Energy Sum Rules. The theoretical basis is the Operator Product Expansion (OPE) of the V+A channel Adler function in which the higher order terms of the leading-twist part originate from a model based on the known structure of the leading renormalons of this quantity. The applied evaluation methods are contour-improved perturbation theory (CIPT), fixed-order perturbation theory (FOPT), and Principal Value of the Borel resummation (PV). All the methods involve truncations in the order of the coupling. In contrast to the truncated CIPT method, the truncated FOPT and PV methods account correctly for the suppression of various renormalon contributions of the Adler function in the mentioned sum rules. The extracted value of the ${\overline {\rm MS}}$ coupling is $\alpha_s(m_{\tau}^2) = 0.3116 \pm 0.0073$ [$\alpha_s(M_Z^2)=0.1176 \pm 0.0010$] for the average of the FOPT and PV methods, which we regard as our main result. On the other hand, if we include in the average also the CIPT method, the resulting values are significantly higher, $\alpha_s(m_{\tau}^2) = 0.3194 \pm 0.0167$ [$\alpha_s(M_Z^2)=0.1186 \pm 0.0021$].


I. INTRODUCTION
The physics of semihadronic τ lepton decays is an important area of QCD, because it describes QCD at relatively low momenta Q m τ ∼ 1 GeV and, at the same time, has high precision experimental results. The latter are principally from ALEPH Collaboration [1][2][3][4], where the spectral function ω(σ) was measured with high precision. 1 The extraction of the value of the running QCD coupling α s (Q 2 ) at such low momenta Q 2 ≈ m 2 τ represents a test of QCD, especially when comparing, via renormalisation group equation (RGE) evolution, with the extraction of the running coupling from experiments at higher energies Q 2 m 2 τ where the perturbative methods of evaluation work very well [5][6][7].
The theoretical framework for the calculation of the QCD corrections r τ to the τ decay ratio R τ ∝ Γ(τ → ν τ hadrons), 2 and of other related sum rules, is well-established [10][11][12]. The perturbative part of the related quark current correlator is known up to O(α 4 s ) [14]. The nonperturbative corrections to r τ are also well understood and were shown to be small [12,15].
This problem of FOPT vs CIPT was addressed in the works [23,24]. There it was argued, on the basis of the large-β 0 (LB) approximation and on numerical evidence, that the truncated FOPT method accounts for certain renormalon cancellations in r (D=0) τ ≡ r τ (m 2 τ ) (D=0) and in related Finite Energy Sum Rules (FESRs), and that the truncated CIPT method does not account for such cancellations. 3 . Such arguments necessarily involve an extension of the perturbative part of the Adler function, d(Q 2 ) (D=0) , beyond ∼ α 4 s so as to account for the theoretically expected renormalon structure. The resulting FOPT and CIPT (truncated) evaluations of r (D=0) τ and of other related FESRs were then compared with the evaluation of these quantities when the Adler function is calculated as the inverse Borel transformation (Borel sum) and the renormalon ambiguity in the Borel sum is fixed by the Principal Value (PV) prescription.
In this work, we perform a QCD analysis of various sum rules related with the strangeless semihadronic τ -decays, following in part the work of Ref. [20]. In order to discern the role of the renormalon singularities, we use an extension of the Adler function d(Q 2 ) (D=0) beyond the order α 4 s , based on the renormalon-motivated construction of Ref. [26]. We apply, in the theoretical Operator Product Expansion (OPE) of various FESRs and of Borel-Laplace sum rules, the (truncated) FOPT and CIPT methods and the Borel sum (PV) method, and then we extract the corresponding values of α s from the ALEPH experimental data. All the sum rules are double-pinched, and the V+A channel of the ALEPH data was used; we believe that these two aspects suppress sufficiently the duality violating effects, cf. [20] (cf. also [12,15,[27][28][29][30][31][32][33]). We further argue (beyond the LB approximation) that in the considered sum rules important renormalon contributions of the Adler function get cancelled in the truncated FOPT approach, and that such a cancellation is in general not expected in the truncated CIPT approach. The Borel sum approach, on the other hand, is expected to sum correctly in the sum rules the main renormalon contributions of the Adler functions. The extracted values of α s appear to be consistent with these considerations; namely, they turn out to be similar in the truncated FOPT approach and in the Borel sum approach, and they are consistently higher in the truncated CIPT approach. For these reasons, our main numerical results for α s are obtained from the combination of the FOPT and Borel sum results. 4 The paper is organised in the following way. In Sec. II we recapitulate the main elements of the QCD sum rules in the context of the semihadronic τ decays, and their relation to the Adler function. In Sec. III we resume the main aspects of the renormalon-motivated extension of the Adler function in the MS scheme, as constructed in Ref. [26]. In Sec. IV we present the specific sum rules (Sec. IV A) to be considered in the numerical analysis, and the methods of evaluation of the Adler function extension: FOPT, CIPT and Borel sum (PV) (Sec. IV B). In the related Appendix we show how certain renormalon contributions of the Adler function get cancelled in the various considered sum rules (FESRs and Borel-Laplace sum rule), at any loop level (i.e., beyond the LB approximation) and in any renormalisation scheme. In Sec. V we then present the numerical results for the extracted values of the coupling α s and of the lowdimension condesates. In Sec. VI we make conclusions, summarise our results and make a brief comparison with the results of other works.

II. SUM RULES AND ADLER FUNCTION
The Adler function D(Q 2 ) is a logarithmic derivative of the quark current polarisation function Π(Q 2 ) where Π(Q 2 ) stands for the total (V+A)-channel polarisation function These functions appear in the quark current-current correlator Π J,µν (q) = i d 4 x e iq·x T J µ (x)J ν (0) † = (q µ q ν − g µν q 2 )Π (1) where q 2 ≡ −Q 2 is the square of the momentum transfer. Further, J µ are up-down quark currents, J µ =ūγ µ d and uγ µ γ 5 d for J = V, A, respectively. In the V +A sum (2), the contribution Π  The polarisation function has a theoretical expression in the form of OPE [34] Π (th) (Q 2 ; µ 2 ) = − 1 2π 2 ln where µ 2 is the squared renormalisation scale, and O 2k V +A are vacuum expectation values (condensates) of dimension D = 2k (≥ 4). The O(a) terms in the Wilson coefficients turn out to be negligible [35]. 5 Using the OPE (4), the corresponding Adler function is obtained using the relation (1) According to the general principles of Quantum Field Theory, the polarisation function Π(Q 2 ; µ 2 ), and thus the Adler function D(Q 2 ), are holomorphic (analytic) functions of Q 2 in the complex Q 2 -plane with the exception of the negative semiaxis, Q 2 ∈ C\(−∞, −m 2 π ). The associated QCD sum rules are obtained then in the following way. If g(Q 2 ) is any holomorphic function in the complex Q 2 -plane, then the integration of the integrand g(Q 2 )Π(Q 2 ) along the closed path C 1 + C 2 presented in Fig. 1 gives zero by Cauchy theorem The closed contour C 1 + C 2 for integration of g(Q 2 )Π(Q 2 ), where the contour radius is R = σmax (≡ σm) (≤ m 2 τ ). C1+C2 dQ 2 g(Q 2 )Π(Q 2 ) = 0, (6) which then leads to the g-function associated QCD sum rule Here, the integration on the right-hand side is in the counter-clockwise direction in the complex Q 2 -plane, and we denoted with ω(σ) the spectral function of Π(Q 2 ) (along the cut) which was measured by OPAL [36,37] and ALEPH Collaborations [1][2][3][4] in strangeless semihadronic τ decays. We will use the ALEPH data as they have less experimental uncertainty; these data are presented in Fig. 2. In the sum rule The spectral function ω(σ) for the (V+A)-channel, measured by ALEPH Collaboration [1][2][3][4], without the pion peak. The pion peak contribution 2π 2 f 2 π δ(σ − m 2 π ) has to be added (where we took fπ = 0.1305 GeV). In the sum rules we will take σm = 2.80 GeV 2 in order to exclude the last two bins with large uncertainties.
(7) the theoretical polarisation function (4) can be replaced by the Adler function (5) by application of integration by parts where, as in Eq. (7), the integration on the right-hand side is in the counter-clockwise direction in the complex Q 2 -plane, D (th) (Q 2 ) is given by the OPE expansion (5), and G(Q 2 ) is The Adler function D(Q 2 ) as a (quasi)observable is a spacelike quantity, i.e., it is holomorphic in the complex Q 2plane with the exception of the negative semiaxis. On the other hand, the quantities (9) are timelike observables, they are functions of the squared energy σ m > 0 (= −Q 2 ). In the case of the sum rules (9), the timelike squared energy σ m is in an intermediate range σ m ∼ m 2 τ ∼ 1 GeV 2 (we have here σ m = 2.8 GeV 2 ). There exist several other timelike quantities in form of integrals of D(Q 2 ) that are phenomenologically important [38], among them: (a) the production ratio for e + e − → hadrons, R(s) [39,40], where the squared energy |Q 2 | = s > 0 is in principle not constrained; (b) the leading order hadronic vacuum polarisation contribution to the anomalous magnetic moment of µ lepton, (g µ /2 − 1) had(1) [41,42], where the dominant momenta Q 2 of D(Q 2 ) are in the deep IR regime Q 2 ∼ m 2 µ (∼ 0.01 GeV 2 ) [43,44].

III. ADLER FUNCTION: RENORMALON-MOTIVATED EXTENSION
In the sum rule (9), the theoretical expression for the Adler function is the OPE Eq. (5), where the leading-twist (D = 0) QCD part is given by the perturbation expansion (pt) where a(µ 2 ) ≡ α s (µ 2 )/π is the pQCD coupling, d 0 = 1 in our normalisation, and we will work in the MS renormalisation scheme. Here, µ 2 = κQ 2 is the renormalisation scale (0 < κ ∼ 1 is the renormalisation scale parameter), and the κ-dependence of the coupling is determined by the (five-loop) MS RGE [45] da(κQ 2 ) where β 0 = (11−2N f /3)/4 (= 9/4 for N f = 3) and β 1 = (1/16)(102−38N f /3) are universal (i.e., scheme independent) in mass independent schemes, and β j (j ≥ 2) depend on the renormalisation scheme. We will always use the five-loop MS RGE when varying the renormalisation scale, independent of the truncation index. The first four terms in the expansion (11) (i.e., the coefficients d j , j = 0, 1, 2, 3) are exactly known [14,46,47]. In Ref. [26], a renormalonmotivated extension of this expansion to all orders was constructed. It was based on the following considerations. The perturbation expansion (11) in powers of a can be reorganised in another expansion where a n+1 (Q 2 ) are logarithmic derivatives which can be expressed in powers of a (by using the RGE) These relations can be inverted and have the form The relations (15) imply linear relations between the expansion coefficients d n and d k where the coefficients k s (n+1−s) are (κ-independent) combinations of the beta-function related coefficients c j ≡ β j /β 0 [26]. We note that d 0 = d 0 (= 1 in our normalisation). If we formally replace in the expansion (13) the logarithmic derivatives by the corresponding powers, a n+1 (Q 2 ) → a(Q 2 ) n+1 , we obtain another associated quantity d( which agrees with d(Q 2 ) (D=0),pt Eq. (11) only at the one-loop level, and is κ-independent only at the one-loop level ]. It turns out that the exact κ-dependence of the coefficients d n (κ) has the one-loop-type form As a consequence, the Borel transform B[ d] of the power expansion (18) has the simple one-loop-type (or: large-β 0 -type) κ-dependence This would suggest that this Borel transformation has the renormalon structure of the form of the the large-β 0 Borel transform of the Adler function , i.e., in terms of single or multiple poles (and not noninteger-multiplicity poles) The ansatz as made in Ref. [26] includes these leading renormalon singularities, as well as the "zero"-multiplicity u = 2 infrared renormalon singularity ∼ ln(1 − u/2): where κ = 1 and the values of the parameters, for the MS scheme case, are We refer to Ref. [26] for details on how these parameter values were obtained. We point out that the model, by construction, reproduces the first four exactly known expansion coefficients ( d 0 = 1, d j , j = 1, 2, 3). Further, the next unknown expansion coefficient (at κ = 1, in MS) is predicted to be d 4 = 338.19 ( d 4 = 37.77); this prediction comes from consideration of this approach in the lattice-related MiniMOM scheme 6 where the number of adjustable parameters was one less (i.e., without the term ∝ d IR 3,1 ). 7 The value of the parameter α was determined on the basis of the knowledge of the subleading Wilson coefficient c D=4 1 of the D = 4 condensate. On the other hand, the values of the other five parameters in the ansatz (23) were fixed by the knowledge of the first five perturbation coefficients d j (j = 0, 1, . . . , 4; where d 4 = 338. 19).
The expression (23) generates the coefficients d k [cf. Eq. (20)]. Then, the coefficients d n are obtained via the relations (17). As shown in Ref. [26], the coefficients d n obtained in this way then lead to the following Borel transform of the quantity d(Q 2 ) D=0 Eq. (11) where the main beyond-one-loop effects are contained in the coefficients γ p and γ p the coefficients b j,k and the residue ratios d X p,k / d X p,k in the MS scheme (cj = 0 for j ≥ 5), with N f = 3. Note that D = 2p for IR renormalons, and D = −2p for UV renormalons.
We recall that c j ≡ β j /β 0 . The numerical values of the coefficients C (D) j,k and the residue ratios d X p,k / d X p,k are given here in Table I (cf. Ref. [26]). Further, the coefficient α is whereĉ (4) . The result (25) was obtained from the expression (23) to a large precision, by generating first the coefficients d n from the coefficients d k via the relations (17) and going up to high n (n max = 70). It is interesting that the form of the expression (25) is also expected by the arguments of the theory of renormalons.
We can interpret the transition from the coefficients d k to the coefficients d n [the relation (17)], or equivalently, the transition from the Borel transform Eq. (23) of the quantity d(Q 2 ; κ) D=0 to the Borel transform Eq. (25) of the quantity d(Q 2 ) D=0 , as a procedure of "dressing" with the beyond-one-loop effects. Nonetheless, we point out that the coefficients d k contain all the information about the quantity d(Q 2 ) D=0 (to all loop levels), because they are in one-to-one correspondence with the coefficients d n , cf. Eq. (17). This in spite of the fact that the Borel transform Eq. (23) of the quantity d(Q 2 ; κ) D=0 behaves under the variation of the renormalisation scale parameter κ as if it were the Borel transform of the quantity d(Q 2 ) D=0 in the one-loop approximation, cf. Eq. (21).
In Table II we present the values of some of the coefficients d n and d n (for κ = 1). In the Table we include the ratios d n /((n + 1)!(−β 0 ) n ) and d n /J(n) (X) , where J(n) (X) (X=0, 1) describe the leading (∼ 1) or next-to-leading (up to ∼ 1/n) asymptotic behaviour factor of d n as follows from the expression containing the UV (u = −1) renormalon contribution [i.e., the contribution to d n from the term containing d UV 1,2 in Eq. (25)] We see from the Table that the two ratios d n /((n + 1)!(−β 0 ) n ) and d n /J(n) (1) converge to specific values at large n (approximately to −0.0221 and −0.0257, respectively), which confirms that the p = 1 UV renormalon contribution is really the dominant contribution to these coefficients at large n. The ratio involving d n in the last column converges even faster if we included the terms O(1/n 2 ) in the asymptotic form (i.e., d n /J (2) (n)). We point out that the starting point for the construction of the higher order coefficients d n (n ≥ 4) of an extended Adler function in [26] (and here) was not the Borel transform of the (D = 0) Adler function, B[d](u), but a renormalonmotivated ansatz for the Borel transform of the auxiliary quantity d, B[ d](u), which has a particularly simple strucure of poles with integer multiplicity, Eq. (23). On the other hand, the works [23,24] construct and use a renormalonmotivated ansatz for the Borel transform B[d](u) (which has poles of noninteger multiplicity) in order to generate the higher order coefficients d n . The authors of [55,[57][58][59][60] generated the higher order coefficients d n by a combination of a renormalon-motivated ansatz for B[d](u) and application of an optimal conformal mapping in the Borel plane. For a classical review on renormalons, we refer to [61], and for some recent developments on the subject of renormalons we refer to [62][63][64][65]. The MS coefficients dn and dn (with κ = 1) of the considered renormalon-motivated Adler function extension: the coefficients dn are generated by the Borel transform Eq. (23) of the extended auxiliary quantity d(Q 2 ; κ) Eq. (18), and the coefficients dn are then generated by the relations (17). The values of the first four coefficients (n = 0, 1, 2, 3) coincide with the exactly known values. See the text for details. n dn dn dn/((n + 1)!(−β0) n ) dn/J(n) (0) dn/J(n) (1)  In the previous Section we described the renormalon-motivated extension of the known truncated perturbation series for the Adler function d(Q 2 ) (D=0) . We will use various methods of evaluation of this function in the sum rule approach described in Sec. II, and will use various weight functions g(Q 2 ) in the sum rules Eq. (6). In the analysis, we will use the ALEPH experimental data, and will extract the corresponding values of the (MS) QCD coupling α s (m 2 τ ).
We will consider the FESRs with moments a (2,n) associated with the following weight functions g (2,n) (n = 0, 1, 2, . . .): We recall that we assume that O D are Q 2 -independent. The weight function G (2,n) (Q 2 ) is related with g (2,n) (Q 2 ) via the relation (10), and the theoretical expression (30d) represents the right-hand side of the sum rule (9) (minus unity) where for the entire Adler function D (th) (Q 2 ) the OPE Eq. (5) was used, and Q 2 ≡ σ m exp(iφ) (−π ≤ φ < +π) on the countour. The coefficient (n + 3)/(n + 1) appearing in the weight functions was used so that the unity in the OPE expansion (5) of the Adler function gives exactly the unity in the contour integration on the right-hand side of the sum rule (9). We will use the sum rules with the above moments for n = 5 ± 1, by assuming that the contributions of the high dimension condensates in Eq. (30d) are negligible. This will allow us to extract the values of α s (m 2 τ ) without consideration of these condensates.
In addition to these sum rules, we will consider the sum rules with the (double-pinched) Borel-Laplace transforms B(M 2 ) (where M is a complex scale parameter) where in Eq. (31d) the dimension D = 2k condensates contribute to the Borel-Laplace In our application of these (double-pinched) Borel-Laplace sum rules, we will use only the real part of the Borel-Laplace, ReB(M 2 ; σ m ), because in this way the D = 4 condensate contribution dominates over the D = 6 condensate contribution when M 2 varies along the ray M 2 = |M 2 | exp(iπ/6), and D = 6 dominates over D = 4 when M 2 = |M 2 | exp(iπ/4). 8 Further, while including the (small) D = 8 contribution, we will neglect higher dimension contributions In general, while smaller scales |M 2 | tend to minimize the duality violations they make the (higher) condensate contributions larger [20]. Further, larger values of |M 2 | lead to large experimental uncertainties [ω (exp) (σ) has larger uncertainties at large σ]. We consider as a reasonable range and we will use this range in our fits. Further, for Ψ ≡ Arg(M 2 ) we will use the rays with Ψ = 0, π/6 and π/4 [we recall that, along the last two rays, specific condensate contributions are suppressed in ReB th (M 2 ; σ m )].
B. Evaluation methods for the Adler function d(Q 2 ) (D=0)

Fixed-order (FOPT)
The basis of the known fixed-order (FOPT) approach is the application of the Taylor expansion to the D = 0 (leading-twist) part d(Q 2 ) (D=0) of the Adler function D (th) (Q 2 ) Eq. (5) on the contour Q 2 = σ m exp(iφ) on the right-hand side of the sum rule (9) When appying this Taylor expansion to the expansion (13) of d(Q 2 ) (D=0) in logarithmic detivatives a n+1 (κQ 2 ), and using the identity which is a direct consequence of the definition (14), we obtain where the φ-dependent expansion coefficients d n (φ; κ) are the following combination of the coefficients d n (κ) of the expansion (13) We note that d n (φ = 0; κ) = d n (κ). When inserting the expansion (37) in the right-hand side of the D = 0 part of the sum rule (9), where where the expansion coefficients r where d 0 (κ) = 1 and the coefficients K n− (σ m ) are the following contour integrals: The expansion (39) represents the FOPT expansion in logarithmic derivatives a n+1 (κQ 2 ) at Q 2 = σ m , which we denoted as ( FOPT). This expansion involves in practice a truncation, say at a N (κσ m ), which we will denote with the superscript ( FO, [N ]). This expansion can be reorganised in terms of powers a(κσ m ) k using the relations (15), and then truncating at the power a(κσ m ) N ; this represent the usual FOPT approach, and we will denote it with the superscript (FO), and its truncated version with (FO; where in complete analogy with the relations (17). We point out that, while the approaches ( FOPT) and FOPT give in principle equal results, in practice it is not so due to the truncation. Namely, both types of series are divergent due to the renormalon-dominated growth of r n and r n coefficients when n increases; a truncation is needed (at a N and a N , respectively), which then gives somewhat different results.

Contour-improved (CIPT)
The contour-improved method is represented by the direct integration along the contour of the integrand d(Q 2 ) (D=0) G(Q 2 ) on the right-hand side of the sum rule (9) where d(Q 2 ) (D=0) has the form of the perturbation series (11) truncated at a specific power a(κQ 2 where d [N ] (D=0);pt is the truncated series Here, the renormalisation parameter κ-dependence appears because of truncation.

Principal-Value (PV)
In this approach, the Adler function is evaluated as the Principal Value (PV) of the inverse Borel transformation where the paths C ± go from u = 0 to u = +∞ within the upper and lower half of the complex u-plane; application of the Cauchy theorem shows that the details of these two paths are irrelevant, because the Borel transform B[d](u; κ) has singularities only along the real axis. For example, we can choose for C + the path going as straight line from u = 0 to u = iε (for any ε > 0) and then parallel to the real axis from u = +iε to +∞ + iε. Another example is the path from u = 0 along a ray u = |u| exp(iφ 0 ) (|u| from zero to +∞) where φ 0 is a fixed angle, 0 < φ 0 < π/2. We note that in the integration (46) in the sum rules, the value of Q 2 is in general complex nonreal [Q 2 = σ m exp(iφ)].
In applying the integration (46), we use for the Borel transform B[d](u; κ) the expression (25), but now at a given general κ-parameter value, and truncated. In practice, this truncation requires to include a polynomial correction form so that the full expansion coefficients d n (κ) of the Adler function are restored. The leading part of the renormalon growth of the coefficients d n (κ) is contained in the PV of the inverse Borel integral of the truncated singular transform B[d](u; κ). The latter transform, at a general κ, is obtained in the following way. The starting point is the Borel transform of the auxiliary quantity d(Q 2 ; κ) (D=0) , cf. Eqs. (23) and (21), which can be written for the general case of the renormalisation scale parameter κ as where These relations are obtained by using the κ-dependence Eq. (21) and the expression Eq. (23), performing the corresponding expansions of exp(u(ln κ + K)) around u = 2, 3, −1, and ignoring the terms This truncation means that we include in the expression (48) only the singular contributions. We note that for κ = 1 [κ = exp(− K)] the values of the residues d X p,k ( κ) reduce to the values d X p,k given in Eq. (24). Then, according to the data of Table I, the corresponding truncated Borel transform of the Adler function is [in analogy with Eq. (25)] The truncation here consists of not including the terms of higher powers of ( It is this (singular) Borel transform contribution that we use in the evaluation of the PV of the inverse Borel integration, Eq. (46). We recall that the power indices γ p and γ p (p = 1, 2, . . .) are given in Eqs. (26).
The values of the parameters C (D) j,k and the ratios d X p,k / d X p,k are given in Table I (cf. also Table II in Ref. [26]); the latter ratios are independent of the renormalisation scale parameter κ (thus independent of κ). We will use the central values given in Table I. For example, d IR 2,1 ( κ) = 1.7995 d IR 2,1 ( κ), etc. The expression for the evaluation of the D = 0 part of the Adler function in the described PV-approach thus acquires the form where the B[d](u; κ) sing is given in Eq. (50) and the correction polynomial δd(Q 2 ; κ) (D=0) is given in Eq. (47). In practice, it turns out that the correction polynomial values are large as are also the PV integral values in the expression (51). However, the two terms for the (D = 0) Adler function give to the sum rules in general contributions with opposite sign, and the sum of the two terms there is smaller (by about two orders of magnitude) then each term. We point out that the correction polynomial expression (47) appears in the expression (51) because in the sums in the singular part of the Borel transform, Eq. (50), truncations were made. To understand this more clearly, the Borel transform B[d](u; κ) Eq. (50), but in its nontruncated form, implies that the coefficients d n (κ) have the form The truncation consists of neglecting the indicated relative corrections O(1/n k ) at the end of the brackets in Eq. (52), which then gives us the "singular" parts d n (κ) → (d n (κ)) sing , and the correction coefficients δd(κ) n = d n (κ) − (d n (κ)) sing appearing in the correction polynomial Eq. (47). It turns out that the series in the brackets of Eq. (52), in inverse powers of n, are relatively slowly converging for n < 10-20, because of the relativey large values of the numerators there. Therefore, the truncation effect and the correction coefficients δd(κ) n are relatively large for such n. Specifically, when κ = 1, we have |δd n (κ)/d n (κ)| ∼ 10 1 for 0 ≤ n ≤ 3, and |δd n (κ)/d n (κ)| 1 when 4 ≤ n ≤ 10. However, for very large values of n (n > 15), the coefficients δd n (κ) in the correction polynomial (47) become relatively negligible: δd n (κ)/d n (κ) → 0 when n → ∞. 9 Specifically, when κ = 1, we have |δd n (κ)/d n (κ)| < 0.05 when n > 15, and |δd n (κ)/d n (κ)| < 0.025 for n > 25. Despite the fact that the first few terms in the polynomial (47) give large contributions, the sum (47) [and thus the sum Eq. (51)] is better behaved when the truncation index N there is N ≥ 5, as it does not have the leading parts of the renormalon contributions in its coefficients δd n .

A. Double-pinched Finite Energy Sum Rules with high index
When we apply FESR with the moments a (2,n) (σ m ), cf. Eqs. (30), we can see from Eq. (30d) that the D > 0 part of the moment a (2,n) (σ m ) depends only on two condensate values, O 2n+4 and O 2n+6 , which are of high dimension when n increases. We will assume that these high-dimension condensates give a negligible contribution (cf. also Ref. [20]) when n is large. When equating the theoretical and the ALEPH experimental values of these moments, we can extract the QCD coupling value. It turns out that, when n increases, this value appears to stabilise reasonably at n ≈ 5 (however, see the discussion at the end of this Subsection). This is true for each of the previously described evaluation methods (FOPT, CIPT, PV). As a result, we extract the following values: These values were obtained with the truncation orders N t = 9, 10, 10 for FOPT, CIPT, PV. The truncation order in the PV approach refers to the maximum power a Nt in the correction polynomial (47). The truncation order in CIPT and PV approaches was chosen in the following way: it is such order that, when we increase the truncation index from N t − 1 to N t , the variation of the theoretical moment a (2,5) th is minimal. 10 On the other hand, the choice of the truncation order N t = 9 in the FOPT approach was chosen by looking at the stability of separate renormalon contributions (see below and Figs. 3). In the results (53) we separated the various uncertainties according to their sources. The symbol (κ) indicates the variation when the renormalisation scale parameter κ (= µ 2 /Q 2 ) is varied around κ = 1, up to κ max = 2 and down to κ min = 0.5. The symbol (N t ) indicates the variation when the truncation number is varied around its central value N t to N t ± 2; (n) indicates the variation when we extract α s from a (2,n) with n = 5 ± 3 (variation of α s when n varies is not weak, cf. Table III, therefore we take δn = ±3). On the other hand, the symbol (d 4 ) in the results (53) indicates the uncertainty due to the d 4 coefficient, where we vary the coefficient d 4 around its central value (as predicted by the considered renormalon-motivated extension) d 4 = 338.19, where we chose this variation to be d 4 = 338.19 ± 338.19; the resulting variation of this type is obtained by using N t = 5, i.e., when the last term in the truncation includes d 4 and keeping all the other parameters of the model unchanged. 11 Finally, the symbol (amb) in Eq. (53e) represents an estimate of uncertainty due to the Borel integration ambiguity for the Adler function. 12 We mentioned earlier, in the text after Eq. (24), how the value d 4 = 338.19 was obtained. We recall that the six parameters (24) Table I would be extracted. Therefore, we decided to estimate such uncertainties of α s (m 2 τ ) from the d 4 -variation in the simpler way as described in the previous paragraph. This means that we kept the (exactly) known coefficients d j (j = 0, . . . , 3) unchanged and took (artificially) N t = 5; in PV approach the values of the parameters of the renormalon model, Eqs. (24) and Table I, were kept (artificially) unchanged, and in the correction polynomial (47)  show that around this index (N t = 6 ± 1) we have partial cancellation of the strong instabilities from the u = 2 and u = 3 IR renormalons. On the other hand, N t = 9 is approximately the index where both of these contributions 10 We evaluate the differences |a (2,5) th (Nt) − a (2,5) th (Nt − 1)| with increasing Nt = 4, 5, . . . and look for such Nt where this difference is minimal. 11 In the PV approach, we varied δd 4 of the correction polynomial (47) (with Nt = 5) around its central value by ±338. 19. 12 This variation was obtained from the variation of the Adler function (51) by ±δd( where the integrand is the same as in Eq. (51) [cf. also Eq. (50)]. 13 We separate these parts in the coefficients dn = d IR2 (20) and (23). give stationary points (minimum and maximum, respectively). This is the reason why in FOPT approach we chose N t = 9. The u = −1 UV contribution is suppressed in a (2,5) (σ m ), as seen in Fig. 3(c), and this behaviour is expected on theoretical grounds as explained in the Appendix A 1.
Instead of the FOPT method Eq. (42), we could apply the tilde-variant ( FOPT) method Eq. (39). The results are similar to those of the usual FOPT method, Eqs. (53a)-(53b): α s (m 2 τ ) ( FO) = 0.314 +0.011 −0.031 , now with N t = 6 the optimal truncation index (the u = 2 and u = 3 IR contributions have local extremes at N t = 6, minimum and maximum, respectively; the u = −1 UV renormalon contribution is negligible). The uncertainties are significantly larger, though, due to the larger uncertainties of the type (κ) and (d 4 ). 14 This has to do with the fact that for Q 2 ∼ σ m (∼ 1 GeV 2 ) the ratio a n (Q 2 )/a(Q 2 ) n [= 1 + O(a)] becomes large 15 for n ≥ 4, and thus the series of terms d n a n+1 behaves in pQCD at such energies considerably worse than the power series of terms d n a n+1 . For this reason, we will use the FOPT method in this work only for illustrative and comparative purposes. In Figs. 4(a) 19) may indicate that the taken uncertainty δ d 4 (≡ δd 4 ) = ±338.19 of the coefficient d 4 is possibly too large. 15 We believe that this is related with the vicinity of the Landau singularities of the pQCD coupling a(Q 2 ) at such low |Q 2 | ∼ 1 GeV 2 ; such numerical problems do not occur in holomorphic variants of QCD in which the coupling a(Q 2 ) has no Landau singularities, cf. Ref. [26].   of α s from specific moments a (2,n) (σ m ) with high n (such as n = 5 ± 1), under the assumption that the corresponding high-dimension condensate contributions are negligible [cf. Eq. (30d)], is a realistic approach. For this reason, in Table III we present the values of α s (m 2 τ ) extracted from moments a (2,n) (σ m ) (with the condensate contributions neglected) for a wide range of n, in the FOPT, CIPT and PV approaches, with the corresponding truncation indices N t = 9, 10, 10. Only the experimental uncertainties are included in the Table. We can see in the Table that the extracted values of α s continue to grow when n increases beyond n = 5. In the last line of the Table we present the values of α s (m 2 τ ) extracted in these approaches from the moment a (0,0) (σ m ), i.e., the moment which has the weight function a simple constant g (0,0) ( This moment has no condensate contributions. In the Table we can see that the values of α s (m 2 τ ) extracted from a (2,n) (σ m ) at very high n (> 20) appear to increase toward the values extracted from the moment a (0,0) (σ m ). This can be understood in the following way. If we denote x ≡ q 2 /σ m = −Q 2 /σ m , the weight function g (2,n) of the moment a (2,n) , Eq. (30a), in the limit n → ∞ and for |x| < 1 converges to the constant 1/σ m = g (0,0) As we see, the double-pinch factor (1 − x) 2 in the weight functions g (2,n) , which was there in order to suppress the duality violation effects, tends to disappear in the large n limit. This, together with the results in Table III, where in the theoretical part we included, in addition to the D = 0 contribution, also the D = 4, 6, 8 contributions, cf. Eq. (33). One of the advantages of using these sum rules, in comparison to the (high-index) FESRs of the previous Section V A, is that we have now an additional continuous complex parameter M 2 . As argued in Sec. IV A, the argument Arg(M 2 ) = Ψ of these parameters [M 2 = |M 2 | exp(iΨ)] is preferrably in the range 0 ≤ Ψ < π/2, and we use specifically Ψ = 0, π/6, π/4, and |M 2 | in the range [0.9, 1.5] GeV 2 [cf. Eq. (34)]. In practice, we cannot ensure the equality (56) for a continuous set of values of M 2 . Therefore, we decided to minimise the following sum of the squares of deviations between the theoretical and experimental values: where M 2 α is a specific sufficiently dense set of points along the rays with Ψ = 0, π/6, π/4 and 0.9 GeV 2 ≤ |M | 2 ≤ 1.5 GeV 2 . In practice, we chose 11 equidistant points along each of the three rays, i.e., the sum (57) contains 33 terms. 16 The various uncertainties are of the same type as those explained in the previous Sec. V A, Eqs. (53). The truncation numbers N t were chosen in a somewhat similar way as in Sec. V A. Namely, we consider the moments a (2,0) (σ m ) and a (2,1) (σ m ) as functions of N t , using in these moments, at each N t , the corresponding α s and O D values obtained from the mentioned global fits of Borel-Laplace sum rules with the same N t . The optimal truncation index N t is then determined to be such at which the best stability of these moments is achieved (cf. also the discussion of Figs. 6-9 later on). This time the truncation numbers turn out to be N t = 8, 5, 6 for FOPT, CIPT, PV, respectively; and the variation of N t around these values we take in general as N t → N t ± 2. 18 The various uncertainties are obtained in the same way as in the previous Section, with the exception of the experimental uncertainty which can be regarded here only as an estimate. Namely, for various values of M 2 α , the quantities ReB exp (M 2 α ; σ m ) are correlated with each other in a complicated manner, i.e., their covariance matrix is complicated and its inversion becomes numerically unstable when the set of the M 2 α values increases. 19 This is also the reason why for the minimisation we used the simple sum of squares Eq. (57) and not a sum involving the inverse of the covariance matrix of the Borel-Laplace sum rules. Therefore, for the estimate of the experimental uncertainties, we proceeded in the following way. We evaluated the sum of squares of the type of Eq. (57), for a small number of points M 2 α (two points along each ray, i.e., the initial and final; the sum has thus six terms), and varied each of the quantities α s and O D separately around the value of the minimum of χ 2 , until this value of χ 2 increased by unity. For example, χ 2 ( O 4 ± δ O 4 exp ) = χ 2 min + 1 (≈ 1). The number of terms (six) in these sums of squares was taken so low in order to not subestimate the experimental uncertainties. Nonetheless, as we can see, the estimates of the experimental uncertainties obtained in this way are still significantly lower than the various theoretical uncertainties.
In Table IV we present the results for α s and the condensates. 20 The final uncertainties in the condensate values are obtained in the same way as for α s (m 2 τ ), i.e., by combining various theoretical uncertainties and the experimental uncertainty.
In Figs 18 For the FO method, we take Nt = 8 ± 2. For the CI method, we take Nt = 5 +2 −1 ; the case Nt = 3 is not included as it does not use all the exactly known coefficients d j , and the corresponding extracted value of αs is significantly higher than in the cases of higher Nt > 3. For the PV method, we take Nt = 6 +2 −1 ; the case Nt = 4 is not included because the fit quality is much worse there (χ 2 = 2.2 ∼ 10 0 ). 19 Cf. the discussion in App. C of [66], where unpinched Borel-Laplace was used, in the context of a QCD with holomorphic coupling. 20 Instead of O 4 V +A we present the corresponding values for the gluon condensate, aGG = 6 O 4 V +A + 6f 2 π m 2 π , where 6f 2 π m 2 π ≈ 0.00199 GeV 4 . these moments when the contributions of the condensates are set equal to zero (but α s values are those used in the full moments). Further, the experimental band (based on ALEPH data) is included. We can see that the full moments (i.e., those with the condensates included) are rather stable under the variation of N t (especially at N t = 4-7) and are consistent with the experimental values. In Figs. 6(b) and 7(b) we can see that the relatively best stability of these results under the variation of N t is at N t ≈ 5 for CIPT and N t ≈ 6 for PV. On the other hand, the results without the condensate contributions are unstable under the variation of N t , and in general deviate significantly from the experimental band. We point out that the values of α s and of the condensate values O D V +A (D = 4, 6, 8) were obtained from a global anaylsis involving fits of the theoretical Borel-Laplace quantities ReB(M 2 ; σ m ) to the corresponding experimental bands, i.e., quantities with a significantly different structure than those of the FESR moments a (2,n) (σ m ).
The behaviour of the moments a (2,0) (σ m ) and a (2,1) (σ m ) in the case of the FOPT methods shows qualitatively similar behaviour as in the case of the CIPT and PV methods presented in Figs. 6 and 7. Again, as in the previous Sec. V A, we can apply here, instead of the FOPT approach Eq. (42), the tilde-variant ( FOPT) Eq. (39). The results turn out to be very similar to those of the FOPT approach Eqs. (58a)-(58b): α s (m 2 τ ) ( FO) = 0.307 +0.024 −0.021 , and now N t = 6 is the optimal truncation index (also N t = 5, 7 appear to be acceptable). Again, as in Sec. V A, the uncertainties are significantly larger than in the FOPT method, principally because of the larger uncertainty of the type (d 4 ), cf. footnote 14. There, we will use the ( FOPT) results only for illustrative and comparative purposes.
The results for the moments a (2,0) and a (2,1) as a function of the truncation index N t , in the FOPT and FOPT approaches, are presented in Figs 9(b) indicate that the stability is achieved at N t ≈ 8 for FOPT and N t ≈ 6 for FOPT.
Concerning the (local) stability of the results for the momenta a (2,0) (σ m ) and a (2,1) (σ m ) under the variation of the truncation index N t , one question that appears is whether we get such a stability also when the values of the fit parameters (α s and O D V +A ) are not fitted at each N t but are kept fixed. For this, it is sufficient to consider the D = 0 contributions a (2,0) (σ m ) (D=0) and a (2,1) (σ m ) (D=0) as a function of N t at a fixed value of α s . In Figs. 10-11 we present these results, for the corresponding fixed central value of α s (m 2 τ ) which is chosen as the central value of each corresponding method -cf. Table IV; and α s (m 2 τ ) = 0.3074 for FOPT. We can see that these contributions in general show no local stability under the variation of N t , although Fig. 10(b) indicates that N t ≈ 6, 7 might be a reasonable value for the CI and PV method, respectively. We also notice that Figs. 6(a)-9(a) show that the condensate contributions (corrections) are numerically significant in a (2,0) (σ m ) and a (2,1) (σ m ), which may cast doubt on the (illustrative) analysis in Sec. V A where the central values of α s (m 2 τ ) were extracted from a (2,5) (σ m ) in Eqs. (53) when neglecting the condensate corrections. Nonetheless, closer inspection of Figs. 6(a)-9(a) reveals that the consensate corrections are in general significantly smaller for a (2,1) (σ m ) than for a (2,0) (σ m ) (the only partial exception being the PV approach). This indicates that high dimension condensates probably do not contribute significantly [cf. Eq. (30d)], and that consequently a (2,5) (σ m ) momentum sum rule has only small corrections from condensates (of dimension D = 14, 16).
We point out that the Borel-Laplace QCD sum rules were first introduced in [34], and later applied in the literature, e.g. in Refs. [66,67]; these Borel-Laplace sum rules had no pinch factor (1 + Q 2 /σ m ) n . Part of the analysis in the work of Ref. [20] uses single-pinched Borel-Laplace sum rules, for M 2 > 0; the condensate contributions are not included (but they are included in FESRs), and the extraction of α s with the Borel-Laplace there is always for a specific chosen value of M 2 > 0 at a time.
Concerning the described global fit with Borel-Laplace sum rules, the following question may be raised. The IR renormalon structure of the used extended Adler function includes only the IR renormalons at u = 2 and u = 3, but not u = 4, cf. Eqs. (23) and (25). This would at first suggest that only the first two condensates, O D V+A with D = 4 and D = 6, should be used in the Adler function to counter the corresponding renormalon ambiguities. But we used  Table IV, but now O10 V +A is included in the global fit. The condensates OD V +A are presented in units of the first three condensates (D = 4, 6,8) in the Adler function for our global fit, i.e., one more. The main reason is that the use of only D = 4 and 6 condensates is not enough because of the simplifying assumptions that we made for the OPE structure Eq. (5). Namely, we assumed in Eq. (5) that the condensates are Q 2 -independent. 22 However, for the condensate term with D = 6 (and those with D ≥ 8) this is not correct, as indicated by the structure of the Borel transform of the Adler function in the LB approximation [61,68,69],  On the other hand, if we use in the global fits with Borel-Laplace sum rules only two condensate terms, D = 4 and D = 6 (without D = 8) with Q 2 -independent condensates, then it turns out that the (two) condensates do not stabilise the resulting moment a (2,1) (σ m ) as a function of the truncation index N t . In fact, they make the variation of a (2,1) (σ m ) with N t even worse than for the pure D = 0 parts a (2,1) (σ m ) (D=0) , in stark contrast with the results in Figs. (7) and (9). 23 We recall that a (2,1) (σ m ) depends on D = 6 and D = 8 condensates (when these condensates are considered Q 2 -independent), and our Adler extension formally does not require D = 8 condensate (which would counter the u = 4 IR renormalon pole ambiguity effects). Therefore, the numerical results of our global fits suggest that the D = 8 condensate (Q 2 -independent) in our analysis simulates the role of the effects of the running of the 1/a(Q 2 ) factor in the first term on the right-hand side of Eq. (59), in the Borel-Laplace sum rules and in a (2,1) (σ m ). A global fit analysis using the more explicit form (59) remains outstanding, but we expect it to give results similar to those presented here.
These questions notwithstanding, our global fit analysis with OPE with condensates assumed to be Q 2 -independent, Eqs. (4)-(5), can be repeated by including one more condensate term, of dimension D(≡ 2k) = 10. This then gives us the results presented in Table V. The values of index N t were kept unchanged in comparison to Table IV. If we determined N t , as in the previous case of O 10 V +A = 0, as the value at which the resulting momenta a (2,0 )(σ m ) and a (2,1) (σ m ) are least N t -sensitive, then we would obtain in the present case (when O 10 V +A is varied): for FOPT N t = 8 (unchanged); for CIPT N t = 4-5, and for PV N t = 5. Nonetheless, here we kept N t unchanged (N t = 8, 5, 6 for FOPT, CIPT, PV, respectively), so that the comparison of the results of Tables IV and V gives us the effects of the OPE-truncation change only (D max = 8 → 10), without interference of the effects of the N t -change. The resulting uncertainties of the extracted values of α s (m 2 τ ) under the variation of N t are anyway similar in the two cases D max = 8 and D max = 10. 24 Comparison of Tables V and IV shows that the OPE-truncation effects are moderate: the values of α s (m 2 τ ) change by less than the uncertainties given in Table IV, and even the values of the condensates in most cases change by less than 50 percent. The values of condensates of dimension D = 8 and D = 10 in Table V indicate that their contributions to sum rules are small and tend to cancel each other. 25 We can estimate the OPE-truncation uncertainty in the extracted values of α s (m 2 τ ) as the difference between the corresponding central values in Tables V and IV, and We can see that, due to the OPE truncation effects, the uncertainty of α s (m 2 τ ) increases moderately in the FOPT case, and remains practically unchanged in the CIPT and PV cases.
As in the case of D max = 8, we can evaluate in the case D max = 10 the FESR momenta a (2,0) (σ m ) and a (2,1) (σ m ) as a function of N t , and obtain results analogous to those in Figs. 6-9. We will not present such Figures, but it turns out that now the stability of these momenta under the variation of N t is even stronger, and the agreement with the experimental values is even better.
The behaviour of the extracted values of condensates with increasing dimension D is qualitatively similar to that in the work of [20] where various pinched FESRs were applied to the ALEPH τ -decay spectral functions. The dependence on the OPE-truncation variation (variation of D max ) is somewhat milder in our analysis, though. The fact that the extracted values of condensates with D ≥ 8 in our analysis are small, or tend to cancel each other, is possibly related with the fact that our D = 0 Adler function extension contains the first two IR renormalons (u = 2 and u = 3), cf. Eq. (23) [cf. also Eqs. (25) and (50)]. However, also the FESR fit results of [20], in the V + A channel, show a similar trend when D max = 8 or 10.
The works of [21,71] (cf. also [72] where R ee (s) is used), on the other hand, give for (sufficiently pinched) FESRs the solutions of OPE with considerably larger absolute values of the condensates O D V +A , for many terms (up to D = 16). 26 This shows that there are at least two very different sets of OPE solutions to the τ -decay data, which correspond to the same or approximately same spectral function ω(σ) V +A (ALEPH) for (sufficiently pinched) FESRs. The results of the works [20] and [21,71] suggest that the duality violations are well suppressed in the V + A channel for FESRs which are at least doubly-pinched. 27 In our global fit analysis, we used doubly-pinched Borel-Laplace sum rules, whose weight functions g M 2 (Q 2 ), Eq. (31a), are additionally suppressed in the timelike limit by the exponential factor exp(cos(Ψ)Q 2 /|M 2 |) → exp(− cos(Ψ)σ m /|M 2 |) (where we took: Ψ ≡ arg(M 2 ) = 0, π/6, π/4).
where the characteristic functions G  23); they involve simple positive or negative powers of t and ln t (cf. [26] for details). The main difference between the FOPT, CIPT and PV methods, on one hand, and this evaluation method, on the other hand, is that this method does not involve truncation. However, when this resummed version, d(σ m e iφ ) (D=0);res , is used in the sum rules, e.g. in the contour integrals (30d) and (31d), the integrations over t at 0 < t < 1 involve the (pQCD) coupling a(te − K σ m e iφ ) at low momenta. For small φ ≈ 0 and small t this means that the integrations are performed close to the Landau cuts of a(Q 2 ) in the complex Q 2 -plane, i.e., at 0 < Q 2 < Λ 2 Lan. (∼ 0.1 GeV 2 ), and this makes the evaluation numerically unreliable. The extracted values of the parameters also indicate this problem. Namely, the central values of α s in the a (2,5) -approach and in the global fit approach with this resummation method are disparate, α s (m 2 τ ) = 0.377 and 0.246, respectively. We will not use these results, as they are significantly affected by the mentioned problem of Landau singularities. 28 In this context, we mention that this resummation approach works well when the QCD coupling has no Landau singularities [26].

VI. SUMMARY OF THE RESULTS AND COMPARISON WITH LITERATURE
The main results of the paper are in Eqs. (60). For the purpose of additional comparison of different methods (FOPT, PV, CIPT), the results in Eqs. (53) are also important.
We can argue that the FOPT and PV methods have the following feature in common: (a) the FOPT [or ( FOPT)] perturbation series for the sum rules, as argued in the Appendix, explicitly have the leading renormalon contribution of the Adler function d(Q 2 ) (D=0) suppressed in them; 29 ; (b) the PV approach in the sum rules isolates the dominant parts of the contributions from the renormalon singularities (UV and IR) of the Adler function, and resums them with the PV convention, while the perturbation series of the correction part in this approach is largely free of the renormalon contributions.
On the other hand, the CIPT approach to the sum rules keeps unchanged the entire coefficients d n of the perturbation series of the Adler function in the resummation of the sum rules, thus importing the strong renormalon-dominated divergence of d n 's (when n increases) in the sum rule evaluation. It is true that the CIPT approach also transforms the powers a(Q 2 ) n [or the log derivatives a n (Q 2 )] of the Adler function into different functions via contour integration with specific weight functions, but this change in general does not account for the renormalon cancellations which are neither reflected in the (unchanged) expansion coefficients d n of the CIPT series. We believe that these aspects are the main reason why the extracted values of α s from the (truncated) CIPT approach differ significantly from the (truncated) FOPT and PV methods (while the latter two methods give mutually similar results). These conclusions are valid not just in the analysis of Sec. V A of the moments a (2,5) (σ m ) (D=0) where the condensate contributions were neglected, but also in the analysis of Sec. V B where the first three condensates were included.
In the Appendix we argued that the FOPT (and FOPT) expansion of the moments a (2,n) (D=0) (with n large, such as n = 5) has the UV renormalons (at u = −1, −2, . . .) as well as some of the IR renormalons (at u ≥ 2) suppressed by one power, in comparison to the Adler function d(Q 2 ) (D=0) . Specifically for the first IR renormalon (at u = 2) of the Adler function this implies that it is almost entirely supressed in the FOPT evaluation of the moments a (2,n) (D=0) for n ≥ 1. 30 On the other hand, the FOPT (and the FOPT) expansion of the Borel-Laplace sum rules ReB(M 2 ) has only the UV renormalons suppressed by one power, in comparison to the Adler function d(Q 2 ) (D=0) ; but the IR renormalons are not suppressed. Therefore, one might expect that the FOPT (and the FOPT) global fit analysis with Borel-Laplace sum rules would give us more unstable results and less reliable value of the extracted α s than such an analysis with the moments a (2,n) (D=0) (n = 5). However, the inclusion of the condensate contributions (with D = 4, 6,8) in such an analysis takes care of the fact that the IR-renormalon contributions of the Adler function are not suppressed in the 28 A general discussion of the Landau singularity problems in pQCD couplings is given, e.g., in [74]. 29 We recall that the leading renormalon contribution is the double-pole u = 1 UV renormalon in the perturbation series of the auxiliary quantity d(Q 2 ; κ) (D=0) , and its analog in the perturbation series of d(Q 2 ) (D=0) . We point out that this suppression of the leading renormalon contribution in the sum rules is true not just in the large-β 0 approximation, but in the exact approach, as shown in the Appendix. 30 This means, the renormalon contribution ∼ 1/(2 − u) γ 2 in the Borel transform of the Adler function d(Q 2 ) (D=0) is suppressed to ∼ 1/(2 − u) γ 2 −1 in the Borel transform of a (2,n) (D=0) for n ≥ 1, cf. also Eq. (25) and (50). We recall that in our considered renormalonmotivated Adler function extension there are only renormalons UV1 (at u = −1), IR2 (at u = 2) and IR3 (at u = 3).
FOPT (and the FOPT) Borel-Laplace sum rules [cf. also the more detailed discussion around Eq. (59)]. Our analysis with Borel-Laplace was indeed performed with the (low-dimension) condensate contributions included. Therefore, we can argue that the FOPT (and the FOPT) global fit with Borel-Laplace sum rules gives us reliable extraction of the values of α s (and of the condensates).
The PV global fit with Borel-Laplace sum rules and condensates is also expected to give reliable results, because the renormalon structure of the Adler function is taken into account correctly (in an isolated, resummed form) in such sum rules.
However, the CIPT global fit with Borel-Laplace sum rules and condensates is again expected to present problems, because the truncated CIPT approach does not suppress the leading UV renormalon contributions and deteriorates in a significant way the interplay between the IR renormalon effects (in the D = 0 Borel-Laplace part) and the condensate contributions. In this context, we recall that the truncated CIPT is neither a perturbation power series nor does it represent a resummation of the evaluated quantity (because it is truncated). 31 The numerical results presented in this work [principally Eqs. (53), (58), (60), Table IV] appear to confirm the arguments given above. Namely, the extracted values of α s are closer to each other when the FOPT (and FOPT) and PV evaluation methods are used, while the extracted value of α s becomes significantly larger when the (truncated) CIPT evaluation method is used. This is true in the analysis of Sec. V A where the high order moments a (2,n) were considered, and in the analysis of Sec. V B where Borel-Laplace sum rules with condensates were considered. Therefore, in our main predictions for α s we will include the (truncated) FOPT and PV evaluation methods, but not the (truncated) CIPT method. Further, as argued at the end of Sec. V A, the results (53) are not reliable because of the unaccounted nonperturbative effects (from condensates and the quark-hadron duality violation effects). However, the results (53) serve principally as an additional comparison of the three methods (FOPT, PV and CIPT) as mentioned above.
Furthermore, we believe that the fact that we used a renormalon-motivated extension of the coefficients d n (n ≥ 4) of the Adler function does not introduce large model ambiguities. One reason is that the extension is motivated on the known renormalon structure of the Adler function, and simultaneously reproduces correctly the first four coefficients d n (n = 0, 1, 2, 3). The other reason is that our methods used truncation indices N t (i.e., truncation at a Nt ) which were often low (N t = 5, 6). 32 On the gounds mentioned above, our main results are represented by the global fit results (with the double-pinched Borel Laplace sum rules) of the truncated FOPT and PV approaches, Eqs. (60a), (60c). We obtain our central result by averaging between these two results. This then gives the following averaged results of the global fits: It turned out that the central result of the FOPT approach is practically equal to that of the FOPT approach, but the uncertainties are higher. We did not include the FOPT approach result in the average (62). The uncertainty ±0.0073 in Eq. (62a) was obtained by adding in quadrature the deviation between the average value 0.3116 and the the central value 0.3075 of Eq. (60a), and the minimal uncertainty ±0.0061 of the two results Eqs. (60a) and (60c) (cf. a similar reasoning in Ref. [20]). The result Eq. (62b), at the canonical scale Q 2 = M 2 Z (and N f = 5) was then obtained by RGE-evolution using the five-loop MS β-function [45] and the corresponding four-loop quark threshold matching [76]. 33 If, however, we included in the average also the CIPT result (60b), the average central value and the uncertainties would significantly increase In this context, we note that the relations and differences between the FOPT and CIPT approach in FESRs of the semihadronic τ decays were investigated in the works [75] from the point of view of Borel transforms and Borel 31 Somewhat related arguments for the preference of the FOPT methods over the CIPT methods, in FESRs of semihadronic τ -decays, were presented in Refs. [23][24][25]75]. 32 We recall that the truncation index Nt was determined in each method in such a way that a relative stability of the full moments a (2,0) (σm) and a (2,1) (σm) is achieved under the variation of Nt. 33 The threshold matching was performed at the scales Q 2 thr = κm 2 q with κ = 2, andmq ≡mq(m 2 q ) equal to 4.2 GeV (q = b) and 1.27 GeV (q = c). 0.312 ± 0.007 (FOPT+PV) Boito et al. [21] DV in a (i,j) 0.296 ± 0.010 0.310 ± 0.014 -0.303 ± 0.012 Pich & R.-S. [20] DV in a (i,j) 0.298 ± 0.031 0.312 ± 0.047 -0.302 ± 0.032 sums. A Borel transform was constructed for the CIPT of FESRs, by first rewriting the CIPT series of such FESRs formally as a (FOPT-type) series in powers of a(σ m ), r (CI) n a(σ m ) n+1 , using for the Adler function either the largeβ 0 approximation [68,69] or the renormalon-motivated model of Ref. [23]. 34 It was shown that the resulting Borel transform has a significantly different structure of nonanalyticity than the Borel transform of the FOPT FESRs (for the latter, cf. Appendix A). For example, in the large-β 0 approximation the u = 2 IR renormalon of the Adler function is completely suppressed in the Borel transform of the FOPT FESRs momenta a (2,n) (σ m ) (with n ≥ 1), 35 while this is not the case for the u = 2 IR renormalon in the Borel transform of the CIPT FESRs momenta a (2,n) (σ m ). The Borel transforms of the CIPT FESRs do not reflect the D > 0 structure of the OPE of the Adler function Eq. (5), or equivalently, the corresponding OPE of the FESRs Eq. (30d). This is in contrast with the Borel transform of the FOPT FESRs which do respect this D > 0 OPE structure as explained in Appendix A. The authors of [75] suggest that the CIPT FESRs would require different, nonstandard OPE corrections, i.e., corrections which would not correspond to the (D > 0) OPE corrections (5) and (59) of the Adler function. In the present work, we did not try to implement such nonstandard OPE in the CIPT evaluations. For these reasons, and for the reasons explained earlier in this Section, we consider it correct to include only the FOPT and PV results, leading to Eqs. For comparison, we present in Table VI the values of α s (m 2 τ ) extracted from ALEPH τ -decay data by various groups, using various sum rules and various methods of evaluation. In the Table the results are presented to three digits. The result of Ref. [60] is an update of the result of Ref. [58], and uses a (PV) summation of a renormalonmotivated Borel transform with a conformal mapping. The results from Ref. [19] in the Table are given for their V+A channel analysis. We can see in the Table that the results of the exhaustive analysis of Ref. [20] gave a result for CIPT approach very similar to ours, while their FOPT analysis gave a result significantly higher than ours. The latter occurred principally because in our case of FOPT evaluation the optimal truncation index turned out to be relatively high (N t = 8) which indicates that the (renormalon-motivated) extension of the Adler function beyond the order a 5 (N t = 5) plays a role in our case of FOPT evaluations. Nonetheless, we recall that the FOPT method in our global fits gave a lower index value N t = 6 and a similarly low central value α s (m 2 τ ) = 0.307 (but higher uncertainties). We can see in Table VI that the duality violation analysis of Ref. [21] (cf. also [37,71,79]) gives even significantly lower values of α s . On the other hand, it was argued in Ref. [20] that the uncertainties in this DV-model should be larger.
In a recent work [80], the mentioned DV-model strategy was used in FOPT-evaluated FESRs of semihadronic τ - 34 We note that this sum r (CI) n a(σm) n+1 , strictly speaking, is not a (FOPT-type) perturbation series, because each coefficient r (CI) n in this sum is itself a series in powers of a(σm). 35 In the large-β 0 approximation, the same type of relations are valid for the Borel transform of a (2,n) (σm) as in Eq. (A15) for the Borel transform of a (2,n) (σm). Further, the u = 2 IR renormalon residue is in our considered model numerically significant in B[ d](u), Eqs. (23)- (24), as well as in B[d](u), cf. Eqs. (25), (50) and Table I (first line, last column). 36 This truncation was used in [20] (cf. also [77,78]) where Nt = 5 FOPT and CIPT methods were used and d 4 = 275 ± 400 (at κ = 1).
We use d 4 ≈ 338 ± 338, i.e., the central value d 4 ≈ 338 as suggested by the described renormalon-motivated extension of the D = 0 Adler function.
decays, using an experimental spectral V-channel function based on data from various experiments (ALEPH, OPAL, BABAR, and supplemented by e + e − → hadrons data), and they obtained the result α s (m 2 τ ) = 0.3077 ± 0.0075, which is very close to our result α s (m 2 τ ) = 0.3075 +0.0066 −0.0069 obtained from Borel-Laplace sum rule global fit with FOPT method, cf. Eq. (60a) and Tables IV and VI. In this context, we point out that our analysis used the combined V+A channel of ALEPH data, and involved double-pinched sum rules (a (2,n) and double-pinched Borel-Laplace). We believe that both of these aspects suppress significantly the possible duality violation effects in our analysis.
The Mathematica programs on which our calculations were based are available from the www page Ref. [81].

Acknowledgments
This work was supported in part by FONDECYT (Chile) Grants No. 1200189 and No. 1180344, and Project PIIC (UTFSM).
Appendix A: Renormalon structure of Adler function-related sum rules In this Appendix we will present relations between the renormalon structure of the Adler function d(Q 2 ) (D=0) and the considered FESRs a (2,n) (σ m ) and double-pinched Borel-Laplace sum rules.
The (D = 0) parts of the theoretical side of the FESRs a (2,n) (σ m ), cf. Eqs. (30b) and (30d), consist of the following elements: where x ≡ q 2 /σ m = −Q 2 /σ m = −e iφ , and the perturbation expansion of d(Q 2 ) (D=0) in powers a n+1 is given in Eq. (11), and in logarithmic derivatives a n+1 in Eq. (13) [cf. Eq. (14)]. The auxiliary quantity d(Q 2 ; κ) (D=0) is defined then via the expansion Eq. (18) as expansion in powers a n+1 . This auxiliary quantity is independent of the renormalisation scale parameter κ only when a(κQ 2 ) runs according to the one-loop RGE, due to the (exact) identities (19). On the basis of the sum rule quantity δ (A3) Further, if we apply in d(σ m e iφ ; κ) (D=0) in the sum rule (A2) the one-loop (1 ) RGE running of a(κσ m e iφ ) around then the quantity δ ( d) x n Eq. (A2), in this (1 )-approximation, turns out to be which means that the Borel transform of δ The expression (A5) was obtained by using in the definition (A2) the (1 ) version of the identity (A3) 38 , exchange the order of integration over du and dφ, and use the (1 )-identity (A4b); the integration over dφ is then trivial leading to the identity (A5). The quantity δ Eq. (A5) is κ-independent because the Adler function auxiliary quantity d(σ m e iφ ; κ) (D=0) in the integrand in Eq. (A2) is κ-independent when a(κQ 2 ) runs according to the one-loop RGE We point out that the coefficients d n (κ) remain unaffected by this replacement a(κQ 2 ) → a (1 ) (κQ 2 ), leading from Eq. (A3) to (A8). We can see also explicitly that the expression Eq. (21)] and therefore by Eq. (A6) and When combining Eqs. (A9)-(A10), we see that the integrand on the right-hand side of Eq. (A5) [cf. also Eq. (A6)] is κ-independent and thus δ ( d;1 ) x n is κ-independent.
where m = 0, 1, 2, and A = σ m /M 2 is a complex constant. The resulting expressions for these integrals are The last identity (A19b) is obtained as the limit of the expression (A19a) when u = N + and → 0 (with N ≥ m integer).
In Eq. (A21b) we explicitly wrote down the auxiliary quantity r(σ m ; κ) at any loop level, in order to point out that it has identical coefficients as the one-loop version r (1 ) (σ m ) in Eq. (A21a). These generated coefficients r n (κ) contain the full (i.e., beyond one-loop) information about the original sum rule quantity r(σ m ) whose two variants of the perturbation expansion ['lpt' and 'pt', in analogy with Eqs. (13) and (11) for the Adler function] are r(σ m ) lpt = r 0 a(κσ m ) + r 1 (κ) a 2 (κσ m ) + . . . + r n (κ) a n+1 (κσ m ) + . . . , r(σ m ) pt = r 0 a(κσ m ) + r 1 (κ)a(κσ m ) 2 + . . . + r n (κ)a(κσ m ) n+1 + . . . , where we recall the definition (14) of the logarithmic derivatives a n+1 . Here, r n (κ) and r n (κ) are interpreted as coefficients of the FOPT expansion and of the FOPT expansion of the sum rule quantity r(σ m ), respectively, with the renormalisation scale parameter κ. The two sets are related in the same way Eq. (43) as the corresponding coefficients of the Adler function expansion coefficients Eq. (17). It can be checked that the described construction of expansions (A22) represents the FOPT expansion and of the FOPT expansion of the sum rules, by comparing the obtained coefficients with those obtained in the direct application (via Taylor expansions) of the FOPT and the FOPT expansion as described in Sec. IV B 1, Eqs. (39) and (42).
Furthermore, in order to understand better why the construction, leading via Eq. (A20) to expansions Eqs. (A22), gives the usual FOPT and FOPT expansions, we note that the logarithmic derivatives a n+1 as defined in Eq. (14), 39 This is in contrast to the usual arguments in the literature which refer to the large-β 0 (LB) approximations of physical quantities. The latter approximations give us information only on the LB-parts d  VII: The (MS) coefficients rn and rn (with κ = 1) of the FOPT expansion in powers of a(σm) for the FESR r(σm) = a (2,1) (σm) (D=0) [= rτ (σm) (D=0) ], where the considered renormalon-motivated Adler function extension was used. The values of the first four coefficients (n = 0, 1, 2, 3) coincide with the exactly known values. See the text for details. n rn rn rn/(n!(−β0) n ) rn/(Γ(γ 1 + n)(−β0) n ) rn/(Γ(γ 1 + 1 + n)(−β0) n ) We can regard the powers (a (1 ) ) n+1 in the construction described in this Appendix as an instrument of provisional replacement: (a) in the full physical quantities we replace the (full) couplings a n+1 : a n+1 → (a (1 ) ) n+1 ; (b) thus the Borel transforms of the auxiliary power series can be constructed and the simple one-loop RGE running can be used in the relations involving such Borel transforms; (c) at the end, the inverse replacements (a (1 ) ) n+1 → a n+1 are made [cf. Eqs. (A21a) and (A22a)]. The above arguments also show that the sum rules (which are timelike quantities) have FOPT perturbation expansions (A22b) for which the same renormalon-related arguments [26] can be applied as for the perturbation expansions of spacelike quantities such as the Adler function Eq. (11), except that now, instead of in general complex (and nonnegative) Q 2 , we have Q 2 = σ m > 0. 40 Specifically, if r(σ m ) has a double-pole (DP) or a single-pole (SP) u = −1 UV renormalon, the 'lpt'-expansion coefficients for large n behave as: r n ∼ (n + 1)!(−β 0 ) n (DP) and r n ∼ n!(−β 0 ) n (SP). And the usual perturbation expansion ('pt') cefficients r n behave as r n ∼ Γ(γ 1 + 1 + n)(−β 0 ) n [1 + O(1/n)] (DP) and r n ∼ Γ(γ 1 + n)(−β 0 ) n [1 + O(1/n)] (SP) where γ 1 = 1 − c 1 /β 0 [cf. Eqs. (26) and (29)]. To illustrate this, we present in Table VII, which is analogous to Table II made for the coefficients of the Adler function, the expansion coefficients r n for the FESR moment r(σ m ) = a (2,1) (σ m ) at increasing n, and we can see that at large n they are dominated by single-pole (SP) u = −1 UV renormalon contribution (columns 4 and 5), not double-pole (DP, column 6) where no convergence of the corresponding ratio is seen when n increases. For better visualisation of the behaviour of the various contributions (X=UV1, IR2, IR3) to the expansion coefficients, we present in Table VIII, for κ = 1 the ratios of the lpt-coefficients d X n / d n and the pt-coefficients d X n /d n of the renormalon-motivated Adler function extension d(Q 2 ) (D=0) , and in Table IX the corresponding ratios of the ( FOPT) lpt-coefficients r X n / r n and the (FOPT)  pt-coefficients r X n /r n of the moment a (2,1) (σ m ) (D=0) . The three types of the coefficients d (1 + u) 2 .