Renormalization-group improvement in hadronic $\tau$ decays in 2018

One of the main sources of theoretical uncertainty in the extraction of the strong coupling from hadronic tau decays stems from the renormalization-group improvement of the series. Perturbative series in QCD are divergent but are (most likely) asymptotic expansions. One needs knowledge about higher orders to be able to choose the optimal renormalization-scale setting procedure. Here, we discuss the use of Pad\'e approximants as a model-independent and robust method to extract information about the higher-order terms. We show that in hadronic $\tau$ decays the fixed-order expansion, known as fixed-order perturbation theory (FOPT), is the most reliable mainstream method to set the scale. This fully corroborates previous conclusions based on the available knowledge about the leading renormalon singularities of the perturbative series.


Introduction
Since the 1990s, inclusive hadronic decays of the τ lepton have been acknowledged as a reliable source of information about QCD. In particular, the strong coupling, α s , can be extracted with competitive precision from these decays. Since the works by Braaten, Narison and Pich [1] and later by Le Diberder and Pich [2], which finally shaped the standard strategy to extract α s from this process, several important developments have occured. On the experimental side, the precision has improved a lot thanks to the LEP experiments; the latest (re)analysis of ALEPH data was published in 2014 [3]. On the theory side, our understanding of the theoretical input from QCD necessary to achieve an accurate α s determination has improved as well. In parallel, there was similar progress in the global knowledge about α s from other processes in the past 25 years. The uncertainty in the PDG recommendation for α s (m Z ) went down from about 5% in 1994 [4] to a mere 0.9% in the latest edition [5,6], while individual extractions from the lattice are achieving uncertainties below 1% (see, for example, [7]). Although the extraction of α s from τ decays remains appealing -it is performed at rather low-energies and provides, therefore, a nontrivial test of asymptotic freedom -it must be carefully scrutinized given the state of affairs.
In the last few years, a reassessment of the α s extraction from τ decays was motivated by the publication of the result for the α 4 s correction in the relevant perturbative QCD series, which is the next-to-next-to-next-to-leading order (N3LO) correction [8,9]. This tour de force calculation, a five-loop QCD result involving about 20,000 Feynman diagrams, was completed in 2008 -more than 15 years after the publication of the α 3 s result [10,11]. Since then, many aspects of the extraction of α s from τ decays have been reexamined. In this note, we will focus on the perturbative series, in particular on its renormalization group improvement.
The QCD description of hadronic tau decays must rely on finite-energy sum rules, which exploit analyticity in order to circumvent the breakdown of perturbative QCD at low energies. The theoretical predictions are then obtained from a contour integral in the complex plane of the variable s -the invariant mass of the final-state hadrons. When performing this integration, one must set the renormalization scale. The two most common procedures are known as fixed-order perturbation theory (FOPT), in which the scale is kept fixed, and contour-improved perturbation theory (CIPT) [2,12], in which the scale runs along the contour of integration. The two lead to different series and this difference, which is larger than the error ascribed to each series individually, is one of the main sources of theoretical uncertainty in the extraction of α s from hadronic τ decay data.
Before entering the specifics of τ decays let us remind some basic facts about perturbative expansions in QCD. As discovered by Dyson in 1952, the perturbative series in powers of the coupling in realistic quantum field theories are divergent expansions [13], no matter how small the coupling is. The fact that the first few terms of these series do provide meaningful results, i.e. they seem to agree reasonably well with experiment, led Dyson to conjecture that these series must be asymptotic expansions: a special type of divergent series that are useful in practice. Asymptotic expansions approach the true value of the function being expanded up to a finite order, after which the series starts to diverge. Their usefulness is illustrated by the famous Carrier's rule which states that "Divergent series converge faster than convergent series because they don't have to converge" [14].
The idea is that the series may approach the true value much faster than a convergent expansion, which is actually a fortunate feature in QCD, since the computation of higherorder corrections becomes quickly impractical. It also implies that a good asymptotic expansion, in comparison with a convergent one, can have consecutive terms that decrease less in magnitude when compared to their predecessor, precisely because "it does not have to converge".
The divergence of the series is due to the factorial growth of its coefficients -this behaviour, in turn, can be mapped to regions of specific loop diagrams. The "optimal truncation" of an asymptotic series of this type is often achieved by truncating it at its smallest term. This goes under the name superasymptotic approximation [14]. The error that is made in such an approximation is typically of the order of e −p/α , where p > 0 is a constant and α the expansion parameter. The quantity e −p/α is non-perturbative and does not admit a power series in α, but vanishes when the expansion parameter goes to zero, as one would intuitively expect. 1 This suggests that the issues related to the fact that the series is asymptotic become more prominent when the coupling is larger. In QCD this means lower energies, such as in τ decays where the relevant scale, of the order of the τ mass is ∼ 2 GeV. One should say that these rules do not have the status of theorems, mathematical proofs are rare here. As a matter of fact, there is no proof that the series in QCD is asymptotic to start with, but everything indicates that this is indeed the case.
The discussion about FOPT and CIPT in τ decays and a decision about which one is the most reliable procedure cannot be taken out of this context. They are both asymptotic series (at best), therefore divergent. In particular, some of the arguments put forward in the literature in favour of CIPT [18] 2 based on the relative size of the first few terms of the series are insufficient, since they tend to downplay, or simply ignore, these basic features of perturbative expansions in QCD. Also, in analysing the series it does not make sense to talk about a radius of convergence, because we are dealing with divergent expansions. Inevitably, a final conclusion about the reliability of the two different procedures requires knowledge about higher-order coefficients of the series. In the absence of higher-order loop computations, one must resort to other methods to estimate those terms.
A possibility is to use our partial knowledge about the renormalons of the series. It is well known that the behaviour of a series of this type at intermediate and high orders is dominated by the renormalons close to the origin. Under reasonable assumptions, one can then construct an approximation to the Borel transform of the series 3 using the leading renormalons and match this description to the exactly know coefficients. This allows for an extrapolation to higher orders and one is able to obtain an estimate for the higher-order coefficients. These type of construction has been studied in detail in Refs. [19,20]. The main conclusions that can be drawn from this strategy are twofold.
1. Under reasonable assumptions, i.e., without any artificial suppression of leading renormalon singularities, FOPT is the most reliable method to set the renormalization scale in hadronic τ decays. Because in CIPT a subset of terms, associated with the running of the coupling, are resummed to all orders important cancellations are missed and the series does not provide a good approximation to the "true" valueunderstood as the value obtained from the Borel sum of the reconstructed series.
2. The fact that FOPT is to be preferred is linked to the renormalon singularity as- 1 Given the logarithmic running of αs in QCD the error of the truncated perturbative QCD expansion becomes e −p/αs(Q) ∼ Λ 2 Q 2 p . These non-perturbative power corrections in 1/Q 2 are, of course, related to the higher-order terms in the Wilson's OPE. In the case of τ decays, the leading one is given by the gluon condensate and scales as 1/Q 4 . There is an infinite series of such terms, one for each gauge-invariant operator that contributes to the OPE and, therefore, the QCD expansion becomes a double expansion in αs and in 1/Q 2 [15]. Effects related to asymptotic nature of the latter are related to the so-called duality violations [16,17]. 2 This argument [18] is essentially unaltered since the publication of Ref. [2] in 1992. 3 Essentially its inverse Laplace transform.
sociated with the gluon condensate. Should this singularity be, for some unknown reason, much suppressed then CIPT would be best.
These conclusions, albeit providing strong support to FOPT, are somewhat model dependent since they do rely on the partial knowledge about the renormalons and could be affected by the inclusion or the removal of a specific singularity from the model. It is, therefore, desirable to study this issue from a model-independent point of view in order to corroborate, or to discredit, the results obtained from the renormalon models.
Here we will discuss recent results presented in Ref. [21] where we used the mathematical method of Padé approximants [22] to extract information about the higher-order coefficients of the series. Padé, or rational, approximants are a reliable model-independent tool that has regained importance in recent years and has found applications in many aspects of particle physics [23][24][25][26][27]. In Ref. [21] we applied the method systematically to the problem of estimating higher orders in the perturbative QCD description of hadronic τ decays. We first used the large-β 0 limit of QCD where the series is exactly known to all orders in α s to test the method. This was done having in mind the concrete situation of QCD, i.e., reconstructing the series solely from its first four coefficients. The method has proven to be robust and sufficiently precise to allow for a conclusion about the reliability of FOPT and CIPT, correctly reproducing the fact that FOPT is to be preferred in the large-β 0 limit. We then turned to QCD and applying the same methods reconstructed the higher orders of the series. Our main conclusions were • The results from Padé approximants and its variants are robust. This conclusion is supported both by the tests in large-β 0 and by the fact that we are able to obtain the N3LO coefficient in QCD from the lower order ones with good precision.
• The reconstruction based on the model-independent Padé approximants favours FOPT and lends support to the renormalon models of Refs. [19,20].
• The six-loop coefficient of the Adler function is found to be c 5,1 = 277 ± 51. This result is in line with some other estimates [19,28], but has a smaller uncertainty.
In the remainder of this note, we will review the main results of Ref. [21] to which we refer for further details.
2 Overview of the theory

QCD in hadronic τ decays
Here, we briefly recall the main theoretical ingredients needed for the QCD analysis of hadronic τ decays. We refer to Refs. [19][20][21] for further details.
The main observable in hadronic τ decays is the ratio R τ which represents the total decay width normalized to the width of τ → eν e ν τ . Here, we restrict the analysis to nonstrange channels which allows us to safely neglect effects due to quark masses. There are then two observables R τ,V and R τ,A where the decay is mediated by vector and axial-vector ud currents, respectively. They can be parametrized as where S EW and δ EW are small electroweak corrections and V ud the CKM matrix element, δ NP encloses all non-perturbative corrections both from OPE condensates and from duality-violations. The unity in between square brackets is the partonic result while δ (0) , which is the main object of this work, represents the perturbative QCD corrections.
The relevant quark-current correlators are where |Ω represents the physical vacuum and the currents are parts. Because the correlators depend on conventions related to the renormalization procedure, it is advantageous to work with the Adler function, defined as the logrithmic Exploiting the analyticity of the correlators involved, the perturbative corrections are written as integral in the complex plane with fixed |s| = m 2 τ as [19] where x = s/m 2 τ and W (x) is the weight function determined by kinematics. The perturbative expansion of D starts at O(α s ) and can be cast as where L = log(−s/µ 2 ) and a µ = α s (µ)/π. In this expansion, the only independent coefficients are the c n,1 ; the others can be obtained imposing renormalization group (RG) invariance, and are expressed in terms of the c n,1 and β-function coefficients [19,29]. The logarithms can be summed with the scale choice µ 2 = −s ≡ Q 2 giving where r n = c n+1,1 /π n+1 and the numerical coefficients correspond to the choice µ 2 = Q 2 , N f = 3, in the MS scheme.
To obtain the perturbative corrections to R τ,V /A one needs to perform the integral in Eq. (3). In the process, one needs to adopt a procedure in order to set the scale µ, which enters, implicitly, through Eq. (4). A running scale, µ 2 = Q 2 , as done in Eq. (5), gives rise to the aforementioned Contour-Improved Perturbation Theory (CIPT), where the running of α s along the contour is resummed to all orders. In this case, δ (0) can be written as (6) Another option is to employ a fixed scale µ 2 = m 2 τ , which gives rise to Fixed Order Perturbation Theory 4 . Then, because α s is evaluated at a fixed scale, it can be taken outside the contour integrals, which are performed over the logarithms that appear in Eq. (4) as Therefore, δ FO can also be written as an expansion in the coupling where the coefficients depend then on c n,1 , on the β-function coefficients, and on the integrals J FO n . In QCD, this expansion reads, for N f = 3 and in the MS scheme, where we give the numerical result of the known contributions to the first unknown coefficient.
Because the perturbative series is divergent, it is convenient to work with the Borel transformed series, which can have a finite radius of convergence, defined, in terms of the expansion in α s , as The original expansion can be understood as an asymptotic series to the inverse Borel transform provided the integral exists. In our context, the series R can represent either the reduced Adler function, D of Eq. (5), or δ (except for u = 1), where β 1 is the leading coefficient of the QCD β-function. 5 The UV renormalon at u = −1, being the closest to the origin, dominates the large order behaviour of the series, which must, therefore, be sign alternating at higher orders. As seen in Eq. (5), this sign alternation is still not apparent in the first four coefficients of the QCD expansion in the MS scheme, which are known exactly.

Padé approximants
A Padé approximant (PA) to a function f (z) [22], denoted P M N (z), is defined as the ratio of two polynomials in the variable z of order M and N , Q M (z) and R N (z), respectively, with the definition R N (0) = 1. Let us consider a function f (z) which assumes a series expansion around z = 0 as 5 We define the QCD β-function as The Padé P M N (z) is said to have a "contact" of order M + N with the expansion of the function f (z) around the origin of the complex plane: the expansion of P M N (z) around the origin is the same as that of f (z) for the first M + N + 1 coefficients From the reexpansion of the approximant P M N (z) one can read off an estimate for the coefficient f M +N +1 , the first that is not used as input [23]. Estimates of this type will be of special interest in this work.
The successful use of Padé approximants to obtain quantitative results about the function f (z) requires only a qualitative knowledge about the analytic properties of the function. The PAs can also be used to perform a reconstruction of the singularity structure of f (z) from its Taylor expansion. Convergence theorems exist for the cases of analytic and single-valued functions with multipoles or essential singularities [22]. Even for functions that have branch points the PAs can be used, in many cases, successfully. In these cases, for increasing order of approximation, the poles of the PAs tend to accumulate along the branch cut, effectively mimicking the analytic structure of the function [22].
In this work, most of the times, the role of the function f (z) is played by the Borel transform of the Adler function, defined in Eq. (9). A key feature of the Borel transform, as already discussed, is its singularities along the real axis, the renormalons. It will be of interest to us to study how this singularity structure is mimicked by the PAs. It is important to note that when f (z) is a general meromorphic function some of the poles (and residues) of the approximant P M N (z) may become complex, even though the original function has no complex poles. 6 Such poles cannot be identified with any of the renormalon singularities, but they do not prevent the use of P M N (z) to study the function away from these poles. In fact, in the process of approximating a function with an infinite number of poles by an approximant that contains only a handful of them, the appearance of these extraneous poles is expected to happen [24].
The approximation of functions with branch points and cuts -as is the case for the Borel transform of the Adler function in QCD -is more subtle. In this case, a possible strategy is the manipulation of the series to a form which is more amenable to the approximation by Padés. Let us consider the particular case of a function f (z) = A(z) (µ−z) γ + B(z) with a cut from µ to ∞ with exponent γ and a reminder B(z) with little structure (both A(z) and B(z) are to be analytic at z = µ). Following the method of Baker called D-log Padé approximant [22], we can form PAs not to f (z) but to which turns out to be a meromorphic function to which the convergence theorem applies. The use of appropriate Padé approximants to F (z) determines in an unbiased way both the pole position, z = µ, and the residue, −γ, which corresponds to the exponent of the cut of f (z). No assumption about neither µ nor γ is made; they are determined directly from the series coefficients. The approximation of F (z) by a PA yields an approximant for f (z) that is not necessarily a rational function. To be more specific, the Dlog-PA approximant to f (z) obtained from using P M N to approximate F (z), that we denote Dlog M N (z), is 6 When the meromorphic function is of the Stieltjes type the poles will always be along the real axis. The functions we approximate in this work are not of this type.
where P M N (z) = Q M (z) R N (z) is the aforementioned PA to F (z). Due to the derivative in Eq. (13), the constant f (0) is lost and must be reintroduced in order to properly normalize the Dlog M N (z). In practice, the non-rational approximant Dlog M N (z) can yield a rich analytical structure, in particular the presence of branch cuts -not necessarily present in the function f (z) -is to be expected.
3 Results in large-β 0 Before discussing our results in QCD, we will present results in the so-called large-β 0 limit, which is a good laboratory for the strategy we present here. Results in this limit are obtained by first considering a large number of fermion flavours, N f , keeping α s N f ∼ 1. In this framework, the qq bubble corrections to the gluon propagator must be resummed to all orders. Using this dressed gluon propagator one can then compute all the corrections with highest power of N f at every α s order to a given QCD observable [15]. The results in large-β 0 are obtained by replacing the N f dependence by the leading QCD β-function coefficient (β 1 in our notation) which incorporates a set of non-abelian gluon-loop diagrams. Accordingly, the QCD β-function is truncated at its first term. 7 In this limit, the Borel transform of the reduced Adler function, defined in Eq. (9) can be written in a closed form as [15,34,35] where the scheme parameter C measures the departure from the MS, which corresponds to the choice C = 0. The result clearly exhibits the renormalon poles, both the IR, that lie along the positive real axis, and the UV ones, that appear on the negative real axis. They are all double poles, with the sole exception of the leading IR pole at u = 2, related to the gluon condensate, which is a simple one. It will be important to consider the Borel transform of δ (0) as well which reads [19,21] B[δ (0) ](u) = 12 The analytic struture of this last Borel transform is much simpler than that of B[ D](u). Now all the UV poles are simple poles, because of the zeros of sin(πu). For the same reason, the leading IR pole of B[ D](u), at u = 2, which is simple in large-β 0 , is cancelled in B[δ (0) ](u) -a result first pointed out in Ref. [36] for the Borel transformed spectral function. Our analysis with PAs benefits greatly from these cancellations since the Borel transformed function is now much less singular. 8 A simpler analytic structure can be much more easily mimicked by the PAs. We also note that the leading UV pole has a residue about ten times smaller than in the Adler function counterpart. This, together with an enhancement of the residue of the double pole at u = 3, postpones the sign alternation of the series and enlarges the range of convergence of the Taylor series. PAs constructed to the expansion of Eq. (16) benefit from these features of B[δ (0) ](u) and lead to smaller errors by virtue of Pommerenke's theorem, granting better coefficient's determination [22]. 7 Strictly speaking, the large-β0 limit would be the "large-β1" limit, in our notation. 8 The fact that the only poles that remain double in Eq. (16) are the ones at u = 3 and u = 4 is not a coincidence. This reflects the fact that δ (0) is maximally sensitive to the dimension-six and dimension-eight OPE condensates. Consequences of this general result for the choice of weight functions in αs analyses from τ decays will be investigated elsewhere [37].
The coefficients c n,1 of the reduced Adler function can be reconstructed from the Borel transform by performing the expansion around u = 0 and using Eqs. (5) and (9). The first six coefficients of the Adler function in the large-β 0 limit, denoted D Lβ , read (N f = 3, MS) to be compared with their QCD counterparts given in Eq. (5). We observe that the sign alternation due to leading UV renormalon sets in at the sixth order (in the MS). These coefficients lead to the following large-β 0 FOPT expansion of δ (0) : FO,Lβ (a Q ) = a Q + 5.119 a 2 Q + 28.78 a 3 Q + 156.7 a 4 Q + 900.8 a 5 Q + 4867 a 6 Q + · · · , (18) to be compared with Eq. (8). Now the sign alternation of the coefficients is postponed and sets in only at the 9th order because of the suppression of the leading UV pole in Eq. (16).
In comparison with the results in full QCD, the large-β 0 limit is a good approximation, in the case of the Adler function, only up to α 2 s . However, for δ (0) FO,Lβ this approximation is still good up to the last known term, i.e. α 4 s . The reason for this better agreement lies in the fact that these coefficients depend also on the β-function coefficients -which are largely dominated by β 1 in QCD -as well as on the integrals of Eq. (7).
In Ref. [21] we have performed a careful and systematic study of the use of Padé approximants to obtain the higher-order coefficients of the series of Eqs. (17) and (18). We have verified that the procedure displays convergence and that the leading renormalon poles can be correctly reproduced. We have also discussed how renormalization scheme variations, partial Padé approximants [22], as well as D-log Padé approximants can be used in order to improve the quality of the approximation. Finally, we were able to design an optimal strategy to predict the higher orders based only on the first four coefficients of the series, which are the only ones available in QCD. Here we will focus on the results from this strategy.
The optimal strategy of Ref. [21] exploits the fact that the Borel transform of δ (0) FO displays a much simpler singularity structure. As shown in Eq. (16), this Borel transform does not have the pole at u = 2 and all other poles are simple poles (with the exception of the ones at u = 3 and u = 4). The leading UV renormalon is therefore more isolated from the IR ones. It can be expected that the use of Padé approximants directly to this Borel transform should yield better and more stable results than in the case of the Adler function. We should note that a rational approximant to δ (0) contains enough information to allow for a full reconstruction of the Adler function since the coefficients c n,1 can easily be read off from the FOPT expansion of δ (0) as We start here by applying Padé approximants directly to the series in α s /π, given by Eq. (18). As we have observed, the FOPT series in large-β 0 is rather well behaved and, at intermediate orders, its asymptotic nature is not visible yet. This is mapped into a simpler analytic structure in the Borel plane. It is therefore likely that in this case the approximation of the series by Padé approximants in a Q will lead to a good description. In Fig. 1 (lower left panel) we display an example of the results obtained (detailed numerical coefficients can be found in [21]). The agreement with the exact results is quite impressive, as seen when comparing with the upper panel of Fig. 1.
Another elegant and efficient way to obtain the higher-order coefficients is to resort to D-Log Padés constructed to the Borel transform of δ optimal way to improve the convergence while remaining completely model independent. The success of this strategy can be understood from the study of the function F (u) = d du log B[δ (0) ](u) , introduced in Eq. (13). The leading analytic structure of F (u) is much simpler than that of the Adler function. The poles at u = 0, u = 1, and u = 2 are exactly cancelled out by the presence of π cot(πu) leaving only a leading UV pole at u = −1, an IR pole at u = 3 and a subleading IR pole at u = 4. It is therefore expected that the D-log Padés should perform well in the present case, since the isolated simple poles can be reproduce by the rational approximant without the need of "spending" too many coefficients.
In Fig. 1 (lower-right panel) we present results for a D-Log Padé applied to B[δ (0) ]. The predictions for c 5,1 have a rather small relative error and the sign alternation is well reproduced by the Padés using only the first four coefficients as input. Their Borel integral provide excellent estimates for the true value of the series. However, one must note that the results from the D-Log Padés applied to B[δ (0) ] are less good than those of the Padés applied to series in α s /π. For example, for Dlog Padés the coefficient c 4,1 is typically wrong by a factor of about two while before it was only a few percent off. Nevertheless, the description of the Borel transformed δ (0) by D-Log Padés has the advantage that the factorial growth of the coefficients is automatically reproduced and an asymptotic series is obtained, in line with the exact result. Furthermore, Fig. 1 shows that these small imperfections do not prevent an excellent reproduction of the exact series.
In Fig. 1 we can see that both methods allow for an excellent reproduction of the exact series up to around the 10th order. It is important to stress that the superiority of FOPT, which is well established in large-β 0 , is very well reproduced in both cases even though we use as input only the first four coefficients of each series.

Results in QCD
In large-β 0 approximants constructed to δ FO and B[δ (0) ] resulted optimal. Furthermore, for the known terms of the perturbative series for δ (0) FO in large-β 0 and in QCD the coefficients rather similar. This suggests that the regularity of the series is preserved in QCD, which indicates that it can be well approximated by Padé approximants constructed directly to the series in α s /π as well. Moreover, although Eq. (16) is strictly valid only in large-β 0 , because it relies on the one-loop running of the coupling, modifications to this result would be solely due to higher-order beta function coefficients. We can therefore expect that a suppression of the leading IR singularity at u = 2, as well as a suppression of all the other renormalons except for the ones at u = 3 and u = 4, would survive in full QCD and render this Borel transform more amenable to approximation by rational functions.
We start with Padé approximants applied to the α s /π expansion of δ (0) FO . We begin with a post-diction of c 4,1 using P 1 2 (a Q ) and P 2 1 (a Q ). The results for six higher-order coefficients obtained from these approximants are shown in Tab. 1. The relative error from the central values of c 4,1 is now only ∼ 13%. This is quite remarkable when put into perspective since, before the true value of c 4,1 was computed, a forecast of this coefficient using other methods and including additional information (taking into account known terms of order α 4 s N 3 f and α 4 s N 2 f ) yielded c 4,1 = 27 ± 16 [28,38,39], a central value which was 45% off. This gives an idea of how powerful optimal PAs can be.
We turn then to the results obtained for P 3 1 and P 1 3 which are also shown in Tab. 1. Now, the forecasts of c 5,1 are 304.7 and 301.3, respectively. We note a striking stability of the results for c 5,1 and c 6,1 ; even c 7,1 and c 8,1 are remarkably similar in all of the four approximants considered. The use of the PAs to sum the asymptotic series also leads to consistent result in all cases, as can be seen in the last column of Tab. 1. Our experience from large-β 0 indicates that this stability and the good prediction of c 4,1 strongly corroborate the robustness of the results. We have checked that the use of D-log Padé approximants is also very successful. We are, therefore, in a position to conclude that using PAs to δ (0) FO in QCD is at least as stable as in large-β 0 . We should then investigate the approximants constructed to its Borel transformed.
As in the previous section, the quality of the forecast of c 4,1 as well as stability arguments lead us to conclude that the D-log Padés are the optimal approximants to B[δ   not shown in Tab. 2, leads to slightly lower values for the coefficients (e.g., c 5,1 = 237), but even these apparent instability can be well understood in terms of a partial cancelation between a pole and a zero present in the P 1 1 used for its construction [22]. We, therefore, consistently discard this approximant. It is also interesting to observe that all the D-log Padés of Tab. 2 predict that the sign alternation of the series starts at order 11. This suggests that the UV singularity in QCD is less prominent than in large-β 0 which should postpone the sign alternation, a fact that can be corroborate by scheme variations of the type of [40] as discussed in detail in [21]. Finally, the Borel sum of the series obtained from these D-log Padés is also very consistent (last column of Tab. 2).
The picture that emerges from the results of this section is that the use of δ FO and its Borel transform lead to the best model-independent approximants in QCD -as is the case in large-β 0 . The quality of the predictions of c 4,1 as well as the stability of the results among different approximants signal that we have managed to obtain a robust description of δ (0) and of the Adler function at higher orders.
We extract our final estimates for the higher-order coefficients from the eight approximants of Tabs. 1 and 2 including, thus, those that have only three coefficients as input parameters. By doing so, we take advantage of Padés that belong to different sequences and can obtain a more reliable error estimate for our final coefficients. Since one of the most striking features of these results is their stability, we will not try to favour one approximant over another, even though one could try to inspect their analytic structure in detail with this goal in mind. Our final estimate of the coefficients and of the true value of δ (0) is obtained as the average of the eight results of Tabs. 1 and 2. To these averages we add an error equal to the maximum spread found between the coefficients obtained from two different approximants. This error should certainly not be interpreted in a statistical sense; it gives an interval where the value of the coefficient is expected to lie.
This procedure applied to the six-loop coefficient, c 5,1 , leads to which largely covers all the results obtained from our optimal approximants. Therefore, in a sense, our error estimate could even be considered as too conservative -even if much smaller than other estimates in the literature. For example, in Ref. [19] the estimate c 5,1 = 283 ± 142 is used, while in Ref. [28] one finds c 5,1 = 145 ± 100 (using only partial information about the five-loop coefficient). The value obtained from the principle of Fastest Apparent Convergence (FAC) in Ref. [8] is c 5,1 = 275, remarkably close to our final central value, given in Eq. (20). On the basis of what we know about the series coefficients, it seems extremely unlikely that the six-loop coefficient would not be within these bounds. Results for coefficients c 6,1 and higher are given in Tab. 3. The final values for the Adler function coefficients are extracted with reasonable errors up to c 10,1 . One should remark that due to the α s suppression at these higher orders, an error that seems large c 9,1 c 10,1 c 11,1 c 12,1 (1.6 ± 1.4) × 10 6 (6.6 ± 3.2) × 10 7 (−5 ± 57) × 10 7 (2.1 ± 1.5) × 10 10 in the coefficient does not translate into a very large uncertainty in the sum of the series.
The situation changes only for c 11,1 . For this coefficient, six of the PAs of Tabs. 1 and 2 predict that the sign alternation sets in. However, two of the approximants do not, which leads to the huge error. Therefore, we find some indication that the sign alternation of the Adler function coefficients sets in at the eleventh order (in agreement with [19]). This instability signals that our results cease to be fully reliable at the 11th order. We apply the same procedure described above to obtain an estimate for the true value of the δ (0) using the results in the last columns of Tabs. 1 and 2. Using α s (m 2 τ ) = 0.316 ± 0.010 [5], this leads to δ (0) = 0.2050 ± 0.0067 ± 0.0130, (21) where the first error is the estimate from the spread of the PAs and the second error is due to the uncertainty in α s . This result agrees with other estimates found in the literature using other methods [19,[40][41][42]) With the coefficients of Tab. 3 we are finally in a position to plot, in Fig. 2, the perturbative expansions of δ (0) and compare them with the true value of the series obtained from Eq. (21). The bands in the perturbative expansions of Fig. 2 represent the uncertainty from the series coefficients, given in Tab. 3, while the band in the Borel sum of the series is the first error Eq. (21). The uncertainties we are able to obtain from the optimal Padé approximants allow us to conclude that FOPT is the favored renormalization-scale setting procedure in the case of QCD. The CIPT series, even though it looks more stable around the fourth order, does not approach well the central value of the sum of the series. The recommendation that FOPT is the best procedure in QCD was advocated in Ref. [19] in the renormalon-model context. Here it is reobtained in a model-independent way.

Conclusion
In this work we have used the mathematical method of Padé approximants to obtain a description of the perturbative QCD series for hadronic τ decays beyond five loops. We have discussed strategies to optimize the use of the available knowledge -namely the first four coefficients. The Borel transform of the series can be used to explain why these strategies are so efficient, as can be cross-checked from the exact results in the large-β 0 limit. The method is shown to provide accurate and reliable predictions for the higher orders and for the sum of the series. This can then be used to study the problem of renormalization-scale setting in hadronic τ decays and the result of this analysis is that fixed-order perturbation theory, FOPT, is favoured within our model-independent reconstruction of the series.
Perturbative expansions in QCD, such as the one for hadronic τ decays, are divergent series that are assumed to be asymptotic. Any conclusion about the renormalization group improvement of the series must be drawn in this context, which automatically invalidades arguments based on the "convergence" of the different series. Since a few years, there is solid renormalon-based evidence that FOPT is the best method to set the scale in τ decays [19,20]. In this work, we have used a completely model-independent method, namely the Padé approximants, to reconstruct the higher orders [21]. Our final results are rather similar to the ones obtained from the renormalon-based methods and fully corroborate the conclusions of Refs. [19,20]. Therefore, as of 2018, the evidence in favour of FOPT is significant and makes this procedure, most likely, the best one to be used in phenomenological studies of hadronic τ decays.