Mathematical Aspects of the Asymptotic Expansion in Contour Improved Perturbation Theory for Hadronic Tau Decays

Recently, it was demonstrated that the discrepancy between the fixed-order (FOPT) and contour-improved (CIPT) perturbative expansions for $\tau$-lepton decay hadronic spectral function moments, which had been affecting the precision of $\alpha_s$ determinations for many years, is related to the CIPT expansion being inconsistent with the standard formulation of the operator product expansion (OPE). Even though the problem can be alleviated phenomenologically for the most part by employing a renormalon-free scheme for the gluon-condensate matrix element, the principal inconsistency of CIPT remains. The CIPT expansion is special because it is not a power expansion, but represents an asymptotic expansion in a sequence of functions of the strong coupling. In this article we provide a closer look at the mathematical aspects of the asymptotic sequence of the functions the CIPT method is based on, and we expose the origin of the CIPT inconsistency as well as the reasons for its apparent good convergence at low orders. Our results are of general interest, and may in particular provide a useful tool to check for the consistency of expansion methods that are similar to CIPT.


I. INTRODUCTION
The comparison of moments of the inclusive hadronic τ -decay invariant-mass spectral distribution obtained at LEP [1] with the corresponding theoretical predictions, based on finite-energy sum rules involving the Adler function, is one of the most precise methods to determine the QCD strong coupling α s at the scale of the τ lepton mass [2]. The theoretical predictions for these finite-energy sum rules involve perturbative series for weighted contour integrals over the invariant mass of the vacuum polarization function that encode the main dependence of the moments on α s . For many years a systematic theoretical discrepancy has persisted for these perturbative series [3,4] related to two renormalization scale-setting prescriptions, called fixed-order perturbation theory (FOPT) and contour-improved perturbation theory (CIPT) [5,6]. While FOPT represents a simple expansion in powers of the strong coupling at a certain renormalization scale, CIPT is based on integrations over the strong coupling renormalization scale and is therefore an expansion in nontrivial functions of the strong coupling. In moments for which the leading nonperturbative correction, coming from the dimension-four gluon-condensate (GC) matrix element in the operator production expansion (OPE), becomes suppressed due to the integration with the corresponding weight function (which are the most relevant for α s determinations) the CIPT scale setting leads to systematically smaller values for the truncated perturbative series, resulting in larger values for the extracted strong coupling.
It was claimed already some time ago in Refs. [7,8] based on the study of renormalon models of the Adler function, that the discrepancy between CIPT and FOPT is due to an inconsistency of the CIPT expansion. Recently, in Refs. [9,10], it was shown that for these moments the CIPT expansion, in contrast to FOPT, does not lead to a corresponding suppression of the moment's quartic sensitivity to infrared (IR) momenta, and the associated suppression of the infrared renormalon does not take place. This renders CIPT, as a matter of principle, inconsistent with the standard formulation of the OPE, even though the CIPT series typically exhibits an apparently excellent behavior at low orders. As far as the hadronic τ decay spectral-function moments are concerned, where this issue is numerically dominated by the dimension-four GC renormalon [9,10], the CIPT problem can be alleviated phenomenologically to a large extent by employing a renormalon-free scheme for the GC matrix element [11,12]. Even though the CIPT problem for hadronic τ -decay spectral function moments can therefore be considered as resolved for phenomenological purposes and a quantitative description of the CIPT inconsistency (called the "asymptotic separation") has been provided in Refs. [9,10], the deeper mathematical background of why CIPT turns out to be inconsistent with the OPE has not yet been studied.
It is the purpose of this article to provide such a mathematical examination. An important motivation for our study is that, once the mathematical origin of the problem is specified in a way independent of the application to the τ hadronic spectral function moments, the consistency of other perturbative expansion methods where integrations over the strongcoupling renormalization scale are employed can be checked. For the most part of our study we use the leading logarithmic (one-loop) approximation for the strong coupling evolution and the large-β 0 approximation, where essentially at all stages of our analysis we can rely on fully analytical results. We also check by numerical analyses that all our conclusions are (beyond any reasonable doubt) valid as well in full QCD. Our main finding can be summarized in the following simple statement: "A convergent series in FOPT in general leads to a divergent series in CIPT". We provide the general mathematical criteria on which this statement is based, and also clarify why CIPT frequently appears to provide a more convergent series expansion than FOPT at low orders.
The article is organized as follows: In Sec. II we set up our notation and prove that the set of expansion functions of CIPT indeed form well-defined asymptotic sequences, which ensures that CIPT provides well-defined asymptotic expansions. However, these asymptotic sequences lack an important type of uniformity property due to zeros for coupling values that become ever smaller at large orders. Here we also prove analytically that a generic factorially divergent series contained in the Adler function (related to an infrared renormalon), which yields a FOPT series with a finite radius of convergence when the corresponding OPE correction vanishes, leads to a divergent CIPT series for any value of the strong coupling. This fact was already discussed in Refs. [9][10][11] based on phenomenological finite-order studies, but not proved rigorously based on the all-order behavior of the actual series. In Sec. III we discuss the reexpansion of CIPT series in terms of FOPT and vice versa. We demonstrate that, while the functions defining the asymptotic expansion of CIPT can be represented by FOPT series with a finite radius of convergence, the reverse is not true. In fact, the asymptotic expansion of a single power of the strong coupling in terms of the CIPT expansion functions is divergent for any value of α s . We provide arguments supporting that this property is related to the nonuniformity and the zeros of the CIPT expansion functions. This finding is the basis of the boldfaced statement above and the central result of this work. All these results are based on explicit analytic expressions derived using the leading logarithmic strong coupling evolution and the large-β 0 approximation. Finally, in Sec. IV we provide numerical evidence that these findings also apply in full QCD, where explicit all-order analytical expressions are not available. In Sec. V we conclude. Lastly, we also add Appendix A, where we state a number of mathematical definitions, theorems, and corollaries that we use in the main body of the article concerning the convergence properties of series, series of functions and the concept of asymptotic expansions. Even though we assume that they are familiar to many researchers, we quote them in the article explicitly for completeness, since some of them, such as the Weierstrass double series theorem, are rarely used in high-energy physics applications. In particular we quote the definition for socalled asymptotic sequences which generalize the concept of asymptotic power expansions. This generalization is relevant because CIPT is not an expansion in powers of the strong coupling but in more complicated functions of the strong coupling with a nontrivial analytic structure. We recommend the reader not familiar with this general definition of asymptotic expansions or the Weierstrass theorem to start the article with the Appendix A.

A. Notation and Conventions
We write the perturbative series for the reduced Adler function for strange plus nonstrange light-quark production in the form 1 D(s) = ∞ n=1c n,1 a n (−s) = ∞ n=1 a n n k=1 kc n,k ln k−1 (−x) , for the rescaled strong coupling. In Eq. (2.1) and in what follows we use the shorthand notation a ≡ a(s 0 ) and x = s/s 0 , where frequently in phenomenological applications one has s 0 = m 2 τ . The series for the τ hadronic spectral function moments which are derived from the Adler function are defined by the contour integrals The coefficientsc n,k are defined in the perturbative expansion of the vacuum polarization function, which reads Π(s, µ 2 ) = − 1 4π 2 [L + n=1 a i (µ 2 ) n k=0c n,k L k ], with L = ln(−s/µ 2 ). The reduced Adler function is defined by the relationD(s) = −4π 2 s dΠ(s) ds − 1 and is renormalization scale invariant.
where the contour starts/ends at x = 1 ± i0 (or s = s 0 ± i0) and is a counterclockwise path in the complex x-plane (or s-plane) around the origin [6,13]. The weight function W (x) is a polynomial with the property W (1) = 0 for physical applications. Note that we suppress the dependence of the moments on the physical scale s 0 since it comes solely from the argument of a = a(s 0 ). In the following we also frequently use the shorthand notation δ The CIPT expansion for δ (0) W (x) starts from the perturbative series forD(s) given in powers of a(−s), and the contour integral is done over powers of the complex-valued strong coupling a(−s) = a(−xs 0 ). This yields where The resulting series is an asymptotic expansion in terms of the functions H n,ℓ (a) with fixed ℓ depending on the reference strong coupling a, which is the quantity determined from the comparison to experimental data in strong coupling determinations. The FOPT expansion for δ (0) W (x) starts from the perturbative series forD(s) given in powers of a = a(s 0 ). Complex phases then appear from the powers of ln(−x) in the coefficients of this series, and the contour integral is computed over the polynomials in ln(−x). This yields where The integrals I k,ℓ can be determined analytically and read: where Γ(n, a, b) ≡ Γ(n, a) − Γ(n, b) and I 0,ℓ = δ ℓ,0 . Note that the second equality of Eq. (2.8) stems from Eq. (A2) and the third from the identity in Eq. (A5). The notation (b) k stands for the Pochhammer symbol (b) k = Γ(b + k)/Γ(b). Note that for ℓ = 0 only the leading term of the sum over n in the second line survives, yielding the expression for I k,0 given in the last line. The series δ FOPT,ℓ is a common power expansion in a. Once the CIPT and FOPT moment series are determined, it is not possible any more to change between the FOPT and CIPT expansions through a change of renormalization scale, since the CIPT series represents a coherent analytic weighted combination of complex-valued strong coupling series with different renormalization scales. It is this difference which is at the core of our investigations.
In the large-β 0 approximation the QCD β-function for the evolution of a(µ 2 ) has the simple leading logarithmic form da(µ 2 )/d ln as the relation of the coupling at different scales. In the large-β 0 approximation the Adler function is determined from single-gluon exchange diagrams dressed with infinitely many insertions of one-loop massless-quark vacuum polarization bubbles with n f flavors. Subsequently the replacement n f → −3β 0 /2 is imposed to yield an approximation for the full QCD results including gluonic corrections. This approximation provides the correct O(α s ) next-to-leading order (NLO) QCD corrections and is known to have many qualitative features of full QCD concerning corrections beyond NLO. It has the advantage that essentially all calculations can be carried out analytically to all orders. It is this property which we rely on in most of the following sections of this article. Still, it is essential to reconfirm that the qualitative insights gained in the large-β 0 approximation also apply in full QCD. This is what we address in Sec. IV.
In the large-β 0 approximation the all-order perturbative series for the Adler function can be conveniently written down in closed form [14] where the Borel functions in the brackets need to be Taylor expanded in powers of u and the relation´∞ 0 du u n−1 e − u x = Γ(n) x n is used. The two equalities yield the two expansions in Eq. (2.1) in powers of a(−s) or a. The function B(u) reads [14] is the first-order polygamma function and d 2 (p) = Each (single or double) pole at u = p along the positive real u axis corresponds to an equal-sign factorially diverging asymptotic series contribution in the coefficients of the Adler function, indicating a sensitivity of the Adler function to IR momenta Λ of O(Λ 2p ). These poles (as well as the corresponding factorially diverging series contributions) are called "IR renormalons". Each such IR renormalon contribution has a one-to-one association with a higher dimensional OPE correction term ∼ ⟨O p ⟩/s p 0 , which effectively compensates for the resulting ambiguity of the perturbation series. Here, ⟨O p ⟩ stands for a dimension 2p nonperturbative low-energy QCD matrix element that cannot be determined with perturbative methods [15][16][17].
For a moment with weight function W (x) = (−x) ℓ all OPE matrix-element corrections of the form const. × ⟨O p ⟩/s p 0 for p ̸ = ℓ are eliminated by the contour integration as can be seen from the residue theorem. As a consequence, perturbative expansion methods consistent with the OPE must yield a convergent series, at least within some region of a, for those terms in the Borel function B(u) with a single pole renormalon of the form 1/(p − u) [15][16][17]. In the following two subsections we prove analytically that FOPT satisfies this requirement, while CIPT does not. Therefore, we consider the moment δ (0) ℓ in the FOPT and the CIPT expansions arising from the term B p (u) = 1/(p − u) with p = 2, 3, . . ., which yields the coefficientsc n,k,p = , (2.13) in the Adler function series of Eq. (2.1).

B. FOPT Series
The contribution of the single pole 1/(p − u) renormalon to the FOPT moment series coefficients d FOPT n,ℓ,p , defined using the coefficients of Eq. (2.13) in Eq. (2.6), can be written down immediately from the expressions given above. We found closed analytic expressions for these coefficients, simply carrying out the sums in the rightmost expression of Eq. (2.6), for any ℓ ∈ N 0 and p ∈ N: Note that we have d FOPT 1,ℓ>0,p = 0. Using the expansion of the incomplete gamma function Γ(n, z) for |n| → ∞ quoted in Sec. A 3 we can obtain the large-n asymptotic expression for the coefficients d FOPT n,ℓ,p which we can use to apply the root test according to Theorem A.2. From Eqs. (A4) and (A5) we find where the sum over k on the RHS refers to the k-th leading term as n → ∞ when p ̸ = ℓ. This sum is actually absolutely convergent. It is conspicuous that there is a universal formula for the (asymptotic) expansion in the cases p = ℓ and p ̸ = ℓ as n → ∞. Furthermore, for p ̸ = ℓ the leading large-n expression (for k = 1) does not depend on p, and its only dependence on ℓ is through the factor (−1) ℓ . The latter property implies that for p ̸ = ℓ the leading term as n → ∞ in Eq. (2.15) always cancels for physical weight functions which have the property W (1) = 0. Interestingly, also the first j subleading terms cancel if also the first j derivatives of W (x) vanish at x = 1 (W ′ (1) = W ′′ (1) = · · · = W (j) (1) = 0), i.e., if W is j-fold pinched.

So for the kinematic weight function
From the π n global factor displayed in Eq. (2.15) we see that the limit superior of Eq. (2.16) is obtained also for any linear combination of d FOPT n,ℓ,p̸ =ℓ coefficients in ℓ including those arising from physical weight functions where the leading k = 1 term or potentially subleading terms for larger k cancel. Here we used that the limit superior of | cos(nπ/2)| 1/n or | sin(nπ/2)| 1/n can be obtained from the infinite subseries with even (or odd) n for which the modulus is unity. This proves that for ℓ ̸ = p the FOPT moment series is absolutely convergent for any complex a within the circle of convergence |a| < 1/π. 2 The FOPT expansion of Eq. (2.6) is therefore consistent with the OPE. We also mention that for any region within this circle of convergence the FOPT series is also uniformly convergent according to Corollary A.5.1. We can analytically sum the FOPT series, and for the case of monomial weight functions and positive real a we find: 3 For illustration, we have displayed the truncated moment series δ FOPT,ℓ,p̸ =ℓ for the values (ℓ, p) = (0, 2), (1,2) in Fig. 1 as red dots using a = 0.2256 which corresponds to α (n f =3) s (s 0 ) = 0.315. The dashed horizontal lines indicate the series sums. The visible oscillations of the FOPT series are substantially damped for physical weight functions due to the cancellation of the leading asymptotic contributions mentioned before.
We note that since the size of the strong coupling α s is only known with an uncertainty and furthermore depends on the observable (through the value of s 0 ), we do not discuss the 2 We note that in Ref. [18] the radius of convergence of anomalous dimensions in Soft-Collinear Effective Theory and boosted Heavy-Quark Effective Theory as well as for the MS-mass were studied in the large-β 0 approximation. It was found that the series converge for |a| < 2.5, which represents a much larger region of convergence than for FOPT moments with ℓ ̸ = p. 3 At first sight it may appear that the small a expansion of Eqs. (2.17) and (2.18) is only a divergent asymptotic series as it involves the asymptotic expansion in Eq. (A6) for the incomplete gamma function for a large second argument. However, the coefficients combine to the convergent coefficients d FOPT n,ℓ,p̸ =ℓ when all terms contributing to a single power a n are combined. The same comment also applies to the small a expansions of the H n,ℓ (a) functions in Eqs. (2.22) and (4.7). particular case |a| = 1/π. This case does not have any specific meaning in practice. So when we talk about convergent series, we always refer to absolutely convergent series in the rest of this article, and it is the fact that a finite interval of convergence for a exists that matters for the consistency with the OPE. We also acknowledge that an indirect proof of the convergence using the Borel representation of the spectral function moment series based on Eq. (2.11) with the unexpanded Borel function B p (u) = 1/(p − u) for p = 2, 3, . . . and using the renormalon calculus has been given in the Appendix of Ref. [11]. The proof presented above directly deals with the actual series. 4

C. CIPT Series
In the large-β 0 approximation the CIPT expansion functions H n,ℓ (a) can be readily computed by rewriting the integral of Eq. (2.5) in terms of the phase angle x = e iϕ , which gives where the integral in the second equality is obtained by a change of variable to t = − 1 functions have poles at a = ±i/π since the phase integral becomes singular at one of the boundaries ±π (or one of the t ± is zero and starts at the pole singularity at t = 0). We define the phase integral in the interval [−π, +π] strictly along the real axis. Since

Im[a]
|a| 2 ∓ π this corresponds to a complex path in t that is in the straight 4 The proof based on the Borel transform in Ref. [11] also applies in full QCD.
vertical downward direction. The definition with a fixed ϕ integration path, where the distance to the origin at ϕ = 0 is bounded, ensures that H n,ℓ (a) is analytic at least in some finite neighborhood of the origin a = 0 and that lim a→0 H n,ℓ (a) = 0 from all directions in the complex a plane. In fact, at the origin the H n,0 (a) vanish like a n , while the H n,ℓ≥1 (a) vanish like a n+1 because of the complex phase e iℓϕ in Eq. (2.19). In other words, H n,ℓ (a) = O(a n+1−δ ℓ,0 ) as a → 0. 5 Our path prescription provides an unambiguous definition of the integrals for all complex a besides where there are cuts. Such cuts arise when the strong coupling Landau pole at ϕ = i/a traverses through the path of ϕ that goes along the real axis from −π to +π (or the straight path of t from t − to t + traverses through t = 0). With our definition these cuts are located along the straight lines (−i∞, −i/π) and (+i/π, +i∞) on the imaginary axis of the complex a plane. Any other path definition would lead to the same poles and branch points at a = ±i/π, but to a different curve of the cut connecting the branch points. A path for ϕ deformed into the positive (negative) imaginary plane corresponds to a path for x with |x| ≤ 1 (|x| ≥ 1) in Eq. (2.5). Thus, if the ϕ integration contour would be deformed far away from the real axis, it is possible that the minimal distance of the cut curve to the origin at a = 0 becomes smaller than 1/π. Our definition is the most obvious one, as it is the standard path used for real a, see Eq. (2.5), and leads to branch cuts along the imaginary axis away from the origin. However, these cuts only truly arise when the Landau pole leads to a simple pole in the integrands in Eq. (2.19), i.e., when there is a finite residue that can make a contribution when the Landau pole at ϕ = i/a crosses the integration path. Since the residue has the general form 1 2πi ℓ n−1 Γ(n) e −ℓ/a , cuts only arise for (n, ℓ) = (1, 0) or when ℓ ≥ 1. There are no cuts for ℓ = 0 and n ≥ 2. Except for the poles and the cuts (if they arise), the H n,ℓ (a) are analytic in the entire complex a plane. In any case, for any sensible contour path choice (with ϕ being close to the real axis or |x| being equal or close to 1) in Eq. (2.5)] the H n,ℓ (a) functions are all analytic within the circle with radius 1/π around the origin, where they can also be expanded in a convergent a n Taylor series.
The analytic results for complex a read where in addition we have H n,ℓ (0) = 0 since one cannot simply insert a = 0 in Eq. (2.22). The term proportional to Θ(a) arises from the definition of the incomplete gamma function in Eq. (A1) which otherwise leads to a mismatch related to the residue at t = 0 already quoted above. We remind the reader that the residue becomes relevant due to the use the incomplete gamma functions, which implies a path deformation in the t variable through real +∞. The step function θ(x) appearing in Θ(a) is defined as θ(x) = 1 for x ≥ 0 and θ(x) = 0 for x < 0. We note that the particular dependence of Θ(a) on the step function θ(x)-definition at x = 0 arises due to the convention used for the incomplete gamma function when the second argument is negative real. The result as shown above arises due to the definition Γ(a, x) ≡ Γ(a, x + i0) for negative real x, as it is used e.g. in Wolfram Mathematica [19]. We note that the analytic evaluation of the integrals in Eq. (2.19) in terms of the incomplete gamma function (where the t path is always deformed through positive real infinity) implies that for values of a on the cut curves, i.e. when i/a is in the interval (−π, +π), the ϕ integration is deformed around the pole into the positive imaginary half-plane. Also note that for the first two equalities in Eq. (2.22) we have used Eq. (A3). Interestingly, they provide the asymptotic expansion of H n,ℓ≥1 (a) as n → ∞, where the k = 1 term represents the leading contribution. The analytic properties of the functions H n,ℓ (a) imply that they can be expanded in Taylor series around a = 0 which are absolutely convergent for |a| < 1/π: The expansion coefficients s ℓ n,k are directly related to the expressions for I k,ℓ defined in Eq. (2.7) and read where we note that s ℓ n,n = δ ℓ,0 . Therefore, the s ℓ n,k form upper-triangular matrices if we take n and k as the row and column indices, respectively. Using the third equality of Eq. (2.8), we also find the infinite series representation which is absolutely convergent, but also provides an (asymptotic) expansion as k → ∞, where the j = 0 term is the leading contribution. For ℓ = 0 only the first term in the j-sum survives. Since there is an infinite subsequence in k such that the modulus of the cosine functions is in the interval [δ, 1] for a 0 < δ < 1 and using Stirling's formula for the gamma functions (contained in the Pochhammer symbol) we find lim sup k→∞ |s ℓ n,k | 1/k = π . (2.27) This reconfirms our statement above concerning the convergence radius of the Taylor expansion around the origin of the H n,ℓ functions in the complex a plane. For illustration we have displayed the results for |s ℓ n,k | 1/k as a function of k in Fig. 2(a) for various combinations of (n, ℓ). We have explicitly checked that summing the series in Eq. (2.24) indeed converges to the analytic expression in Eqs. (2.20), (2.21), and (2.22) for all complex a with |a| < 1/π. It is an interesting fact that, similar to the FOPT coefficients d FOPT n,ℓ,p̸ =ℓ , as discussed following Eq. (2.15), the leading and potentially subleading asymptotic terms cancel in linear combinations of ℓ that arise for physical weight functions W (x), but the limit superior remains π even for such linear combinations due to the global factor of π k shown in Eq. (2.26).
From the expressions in Eqs. (2.21) and (2.22) it is straightforward to determine the limit superior of |H n,ℓ (a)| 1/n as n → ∞. The result reads Interestingly, we yet again find that the leading (and potentially subleading) asymptotic contributions as n → ∞ cancel in linear combinations of the H n,ℓ functions that arise for physical weight functions W (x) in analogy to the discussion following Eqs. (2.15) and (2.26). The limit superior given in Eq. (2.28) also applies for these linear combinations, as we can see from the form of the subleading asymptotic terms given in Eq. (2.22).
From Eqs. (2.21) and (2.28) we see that the n scaling of the H n,ℓ functions is in powers of a/ √ 1 + a 2 π 2 . The fact that |a|/ √ 1 + a 2 π 2 < min(|a|, 1/π) for real a, is one of the reasons why the CIPT expansion typically exhibits a better behavior of the series at low orders than the power a n expansion of FOPT, see also Fig. 3, were H n,0 (a)/(a/ √ 1 + a 2 π 2 ) n−1 is shown for 1 ≤ n ≤ 10 in the interval [0, 1]. It is clearly visible that this scaling property does not only arise for large n, but is active for all n values.
It is now straightforward to examine the convergence of the single pole 1/(p − u) renormalon series contribution in the Adler function to the CIPT moment series δ (2.29) We see that the limit superior diverges for any value of a. The outcome is the same for any linear combination of H n,ℓ functions in ℓ including those arising from the physical weight functions as can be seen from the form of the subleading asymptotic terms in Eq. (2.22). We have displayed the truncated moment series δ CIPT,ℓ,p̸ =ℓ , in Fig. 1 for (ℓ, p) = (0, 2), (1, 2) and a = 0.2256 as the blue dots. We can clearly see that CIPT shows a much better apparent behavior at low orders than FOPT (red dots). The CIPT series does in particular not show the oscillatory behavior of FOPT. This is due to the scaling of the H n,ℓ functions just mentioned above, as well as the zeros of the H n,ℓ functions visible in Fig. 3, which further suppress the size of H n,ℓ (a). These zeros, which appear to have a considerable benefit at this point, turn out to play an important additional role in the further considerations of this article. However, the fact that the CIPT series is divergent for any value of a renders δ the CIPT series appears to approach at intermediate orders 10 < n < 15 and the (correct physical) value of the FOPT series (red dashed horizontal line) renders the CIPT method phenomenologically inconsistent because the Adler function series contains IR renormalon contributions. We want to remind the reader, however, that CIPT can still be used as a viable phenomenological method once the dominant IR renormalons (with the smallest values of p) contained in the Adler function are removed by subtractions [11,12] as the problematic aspects of CIPT just discussed can then be suppressed to a negligible level in the presence of the experimental uncertainties in τ -decay spectral data and truncations of the perturbative series. In the absence of any factorially divergent behavior in the coefficientsc n,1 , such that lim sup n→∞ |c n | would be finite (which is not true for the Adler function), the CIPT series in Eq. (2.4) would have a finite interval of convergence and also yield not only an absolutely but also a uniformly convergent series in any subinterval of the interval of convergence. The important property of uniform convergence in such a situation 7 can be seen from Eqs. (2.21) and (2.22). The expressions show that for real a the functions H n,ℓ (a) [or any physical linear combination of H n,ℓ (a)] are bounded by |a|/ √ 1 + a 2 π 2 n−1 for n larger than some integer n. This is illustrated in Fig. 3 for the functions H n,0 (a), which are bounded by the function |a|/ √ 1 + a 2 π 2 n−1 for all n ∈ N. Since the latter expression is bounded in any subinterval, that CIPT series would also be uniformly convergent in any such subinterval according to the Weierstrass criterion for uniform convergence in Theorem A.5.
We conclude this section with the confirmation that the CIPT expansion in the functions H n,ℓ (a) in Eq. (2.4) is indeed an asymptotic expansion from the mathematical perspective. Even though the CIPT expansion has now been used in phenomenological analyses for a 10]. 7 We note that the uniformity of the convergence of a series in functions ϕ 1 (x), ϕ 2 (x), . . . and the property of uniformity of the asymptotic sequences {ϕ n (x)} discussed below are two different mathematical aspects that must not be confused.
number of decades, it is useful to proceed with the discussion on the mathematical aspects of CIPT at this basic level. Since the following considerations are somewhat mathematical, we refer the reader not familiar with the quoted definitions and theorems to Appendix A 2.
According to the general definition of an asymptotic expansion as a → 0, quoted in Definition A.11, it is mandatory that the H n,ℓ functions form an asymptotic sequence {H n,ℓ (a)} = {H 1,ℓ (a), H 2,ℓ (a), . . .} as a → 0 for any complex domain R containing the origin as a limit point. This means that they obey the "little o" order relation H n+1,ℓ = o(H n,ℓ ) as a → 0 for all n ∈ N and for any ℓ ∈ N 0 , see Definition A.9. The definition of the little-order relation is quoted in Definition A.8. The property of an asymptotic sequence then ensures that the coefficients of the asymptotic expansion can be determined in an unambiguous way by the recursion relation formulated in Corollary A.11.1. Since the H n,ℓ functions are nonzero in at least some neighborhood of a = 0 (excluding a = 0, where the H n,ℓ functions vanish), the little-order relation already follows from the property for any n ∈ N. This can be seen from the form of their Taylor expansions in Eq. (2.24) and the fact that the H n,0 (a) vanish like a n and the H n,ℓ≥1 (a) vanish like a n+1 as a → 0. However, the {H n,ℓ (a)} form nonuniform asymptotic sequences, see Definition A.10, as one has to consider smaller and smaller domains around the origin a = 0 such that the ratio shown in Eq. (2.30) is small. This is due to the zeros already mentioned above and also visible in Fig. 3. To see this at the analytic level let us first consider the case ℓ = 0. The ratio of two consecutive H n,0 functions reads . (2.31) From the oscillatory properties of the sine function it is easy to see that, if R contains intervals of the real axis (containing the origin or having the origin as a limit point), there is no neighborhood U ϵ of a = 0 such that the ratio in Eq. (2.30) is smaller than a given small ϵ for all n. In other words, there always exists such a neighborhood U (n) ϵ for any n, but the maximal distance of the points in U (n) ϵ to a = 0 shrinks to zero as n increases if a can be real. This is because the {H n,0 (a)} have zeros besides a = 0 along the real a axis, which happen to approach zero as n increases. These zeros arise from the phase oscillations due to the iaϕ term in Eq. (2.19) and follow the condition (n − 1) arctan(aπ) = kπ with k ∈ Z which holds only if arctan(aπ) ∈ R. Concretely, they are located atã (0) (n, k) = 1 π tan kπ n−1 for − n−1 2 < k < n−1 2 , where k = 0 corresponds to H n,ℓ (0) = 0. While the number of zeros increases with n, the zerosã (0) (n, k ̸ = 0) all vanish like 1/n as n → ∞ for any k. The zeros are displayed graphically in Fig. 4 exemplarily for a number of n values. For ℓ > 0 zeros with analogous qualitative characteristicsã (ℓ) (n, k) arise as well since the factor (−x) ℓ in Eq. (2.19) only provides an additional modulation of the phase cancellations. However, finding an analytic expression for the zeros is harder. We have displayed them in Fig. 4 as well for a number of ℓ values obtained from numerical evaluations. An approximate analytic insight on the zeros for ℓ > 0 can, however, be gained with the large-n asymptotic behavior of the H n,ℓ≥0 (a) functions as n → ∞ already given in Eq. (2.22).

III. HAVING A DEEPER LOOK
After having discussed the analytic properties of the CIPT expansion functions H n,ℓ (a) in the previous section we now proceed to have a deeper look into the mathematical circumstances why the expansion in terms of the asymptotic sequences {H n,ℓ (a)} is divergent in cases where physics requires a convergent perturbative series. The aim is to identify the mathematical property of the functions H n,ℓ (a) that leads to this behavior, independently of the application in the context of the τ hadronic spectral function moments. We remind the reader that the motivation of this analysis is to identify mathematical criteria to scrutinize also other expansion methods where integrations over the renormalization scale are employed.
To this end we consider transformations of the CIPT or FOPT series by reexpanding one of the series in terms of the other expansion functions. We thus also need the expansion of any power a n in terms of the H k,ℓ (a) functions for each ℓ = 0, 1, 2, . . . in addition to the reverse expansion already given in Eq. (2.24). It has the form Note that for ℓ ≥ 1 the relation only matters for n = 2, 3, . . ., and the Kronecker delta δ ℓ,0 shown here and below arises for the reason explained in Footnote 5 in order to make terms that do not contribute explicit. So we consider the series transformations which represent the transformation of a single series into a double series. Note that at this point our notation does not imply any particular ordering prescription for the summation of the double series. So one should think of the emerging double series as two-dimensional arrays. For example, we have ∞ n=1 ∞ k=n e n,k ∼     e 1,1 e 1,2 e 1,2 · · · 0 e 2,2 e 2,3 · · · 0 0 e 3,3 · · · . . . . . . . . . . . .
where the horizontal lines arise from the series of the expansion function used for the series prior to the double series transformation. A particular summation prescription is indicated explicitly by additional parentheses. Thus, in the example ∞ n=1 ( ∞ k=n e n,k ) stands for summing the horizontal series first and the results of these horizontal sums afterwards. On the other hand ∞ k=1 ( k n=1 e n,k ) stands for summing the vertical series first and the results of these vertical sums afterwards. For the cases discussed below, the limits of the sums are a bit more involved due to the appearance of δ ℓ,0 , which is, however, not significant for the main purpose of the discussion.

A. Double Series Transformations
We start by discussing the transformation of a CIPT series into a double series of power terms, see Eq. (3.2). Consider a CIPT series ∞ n=1 c n H n,ℓ (a) that is convergent according to the root criterion for some a = a 0 within the region of convergence. It is then also absolutely convergent as well as uniformly convergent at least in the circle around zero with radius |a 0 |. The property of uniform convergence follows from the fact that the H n,ℓ (a) are bounded as we discussed already at the end or Sec. II C. In addition, the horizontal lines from the expansion of the H n,ℓ (a) functions in powers of a are absolutely convergent in the circle around zero with radius r < 1/π. According to the Weierstrass double series Theorem A.6 these conditions are sufficient so that doing either horizontal or vertical sums first leads to the same result, ∞ n=1 ( ∞ k=n+δ ℓ,0 c n s ℓ n,k a k ) = ∞ k=1+δ ℓ,0 ( k−δ ℓ,0 n=1 c n s ℓ n,k a k ) for |a| < max(a 0 , r). In other words, an absolutely convergent CIPT series can be reexpanded into a convergent FOPT series which sums to the same value. Furthermore, from the discussions in Sec. II we also know that it is possible that a divergent CIPT series can be reexpanded into an absolutely convergent FOPT series. So starting from a CIPT series and reexpanding it in FOPT does not make its convergence property worse, but can even improve it.
Let us now consider the transformation of a FOPT series ∞ n=1 d n a n into a double series of H n,ℓ (a) functions, see Eq. (3.3). It is a straightforward but rather cumbersome task to determine the coefficients of the expansion of a power a n in terms of the H k,ℓ (a) functions given in Eq. (3.1), but the expressions can also be determined in closed analytic form by inverting the triangular matrix of the s ℓ n,k coefficients. The result reads where B n are the Bernoulli numbers. It is now instructive to have a look at the asymptotic expression for these coefficients as k → ∞. Accounting for the factorial asymptotic growth of the Bernoulli numbers, see Eq. (A8), the result reads The rather interesting outcome is that the expansion of a simple power a n in terms of the H n,ℓ (a) functions does not have any region of convergence and is factorially diverging for any value of a. For illustration, we have displayed the expressions for |t ℓ n,k | 1/k as a function of k for different values of (n, ℓ) in Fig. 2(b). So if we start from an absolutely convergent FOPT series, the resulting CIPT series will, in general, be divergent for any a. The simplest example of an absolutely convergent FOPT series is that of a single term a n . We have shown the H k,ℓ expansions of a n truncated at order N for n = 1, 2, 3, 4, ℓ = 0 and a = 0.25 in Fig. 5. We see that all the series look quite good and convergent at low orders, but they eventually diverge. Furthermore, in the range of orders N , where the series stabilize (and closely approach a particular value), the truncated series show a finite systematic discrepancy to the actual value of a n (indicated by the horizontal dashed lines). Hence, for a general absolutely convergent FOPT series the horizontal sums ∞ k=n+(δ ℓ,0 −1) d n t ℓ n,k H k,ℓ (a) in Eq. (3.3) do not converge for any n or ℓ. The sum over k of the vertical series values s k = k+δ ℓ,0 n=2−δ ℓ,0 d n t ℓ n,k H k,ℓ (a) (which by themselves are finite due to the zeros of the array below the diagonal) is in general divergent as well, as we have demonstrated in the examples discussed in Sec. II C. It is possible that the sum of the vertical series values s k is convergent as we have seen just above for the opposite transformation. However, this case is not the general one. For the application to τ hadronic spectral function moments, the latter case would arise in the situation that the Adler function series of Eq. (2.1) were convergent. This situation does, however, not take place due to the unavoidable appearance of IR renormalons. It is also impossible that the divergent behavior of the coefficients t ℓ n,k as k → ∞ may somehow compensate for a divergent behavior in a FOPT series yielding a convergent CIPT series for any value of a as this would contradict our findings from the beginning of this subsection.
In the discussion above we have considered the CIPT expansion in terms of H n,ℓ functions for a single ℓ. In phenomenological applications, linear combinations of ℓ always arise due to the property W (1) = 0 already mentioned several times before. In this case, the same linear combinations of s ℓ n,k coefficients arise in the analogue of Eq. (2.24). The associated inverse t n,k coefficients, however, cannot be determined from the expressions in Eq. (3.5), and no closed analytic expressions in analogy to Eqs. (3.5) can be given due to the many possible forms for the weight functions W (x). Analytic expressions for definite values of n and k can, however, be obtained in a straightforward way. We have checked for many physical weight functions, including the important kinematic weight function W τ (x) = (1 − x) 3 (1 + x) that the limit superior for the resulting expressions |t n,k | 1/k as k → ∞ is divergent just as for a single ℓ. In particular, we do not find any cancellations of the kind discussed previously below Eqs. (2.16), (2.26) and (2.28) for the case of physical weight functions that change the outcome of the discussion for a single ℓ.

B. On the Origin of the Behavior of the CIPT expansion
The overall observation from the previous considerations is that using the CIPT expansion can turn a convergent FOPT series into a factorially divergent one, while the opposite can never ever happen. This means that CIPT is a priori inconsistent with the OPE and must be applied with great care, see Refs. [11,12]. The diagnostic instrument to probe this behavior is to analyze the convergence properties of the series in Eq. (3.1), where a simple power a n is expressed as a series in H k,ℓ (a) functions. This is the approach as to how the consistency of other expansion methods with the OPE, which are based on functions of the strong coupling, may be tested as well. If the series for a simple power a n is showing a factorial divergence or does not have a finite region of convergence, inconsistencies with the OPE may arise. Even though this diagnostic tool is quite efficient and straightforward to apply, it would be instructive to have a more direct qualitative insight into which particular property of the H n,ℓ (a) function is actually responsible for this behavior, as this knowledge may be highly useful in practice. It is the purpose of this section to explore this question.
An essential difference between the asymptotic power sequence {a n } and the asymptotic sequences {H n,ℓ (a)} is that the latter are nonuniform. The uniform convergence of series expansions in functions is an important aspect, as uniformity states that there is a certain level of global rapidity of the convergence valid for the whole interval for which the expansion is defined. The uniform convergence of the series expansion of a function F (x) = ∞ n=1 f n (x) in terms of the functions f n (x) is frequently a sufficient condition such that doing a linear operation on F (x) is equivalent to doing this operation on the individual f n (x) and a subsequent summation. Uniform convergence is frequently assumed to be valid without proof in many phenomenological particle-physics applications as it is frequently also true. Even though the H n,ℓ functions are bounded and series of the H n,ℓ functions can be uniformly convergent, as we mentioned in the discussion below Eq. (2.29), the nonuniformity of them as asymptotic sequences {H n,ℓ (a)} appears to disrupt their convergence properties in a serious way when a convergent power series is transformed into a CIPT series.
As we have seen in Sec. II C, this is because the ratio H n+1,ℓ (a)/H n,ℓ (a) to be small in a region around a = 0 requires the size of the region to shrink with 1/n when n increases. This is related to the zeros in the functions H n+1,ℓ (a) which approach a = 0 like 1/n when n increases. We now explore whether these zeros, which are actually one reason why the CIPT expansion exhibits a quite rapid convergent appearance at lower orders, could be the origin of the problem. In the following we consider a number of simple toy asymptotic sequences to gain some more concrete insight concerning the answer of this question.
Let us discuss first the series transformation arising from a usual change of renormalization scale. Starting from a series in the powers a n = [a(s 0 )] n we can consider the reexpansion in terms of the powers [a(µ 2 )] n as an expansion in functions h ∞ k=n s L n,k a k can be written down immediately and read s L n,k = (−L) k−n Γ(k) Γ(n)Γ(k−n+1) . The coefficients of the inverse relation a n = ∞ k=n t L n,k h k (a) are trivial to compute and read t L n,k = (−1) k+n s L n,k . We have lim sup k→∞ |s L n,k | 1/k = lim sup k→∞ |t L n,k | 1/k = |L| and both series are absolutely convergent for |a| < 1/|L|. We have displayed |s L n,k | 1/k = |t L n,k | 1/k for |L| = 1 as a function of k in Fig. 6. The {h L n (a)} form a uniform asymptotic sequence since |h n (a) at a = −1/L, but they are independent of n as well. So, as we can conclude from the Weierstrass double series Theorem A.6, using a fixed-order expansion at a different renormalization scale does not turn an absolute and uniformly convergent series into a nonconvergent one or may change the value of the series from the mathematical perspective, as long as L is not too large. The usual approach of renormalization group improved calculations is to adapt L such that the convergence happens in a rapid way.
Let us now consider the toy expansion functionŝ n,k a n and the series of a n in terms of thê h (m) n functions as a n = ∞ k=n t (m) n,kĥ n (a)} form asymptotic sequences, but they are nonuniform due to the zeros or pole singularities at a = 1/(ξn). It is obvious that a series inĥ (m) n functions for negative m has bad properties since the singularities make it impossible to approximate any function that is continuous in some finite neighborhood of a = 0. In Figs. 7(a) and 7(b) we have displayed the values for |s n,k are well-behaved. In this case, this does not come as a surprise. However, for positive m the situation appears not to be that bad due to the absence of the singularities and the emergence of zeros. So, let us have a closer look at the case m = 1. We have s (1) n,k = δ n,k − ξnδ n+1,k . The series only has two terms and therefore trivially converges for all a. However, the coefficients of the inverse series read t (1) n,k = Γ(k)ξ k−n /Γ(n) with t (1) n,k = 0 for k < n. We see immediately that this expansion is factorially divergent and does not have any finite region of convergence. For illustration we have displayed |t (1) n,k | 1/k for n = 1, 5 as a function of k in Fig. 7(b) as well. The interesting insight gained by the toy expansion functionsĥ (m) n (a) is that zeros (for positive m) as well as singularities (for negative m) yield the same badly diverging behavior for the t (m) n,k coefficients. We have tested a number of other toy expansion functions forming nonuniform asymptotic sequences due to zeros or singularities approaching a = 0 for large n, all leading to factorially diverging t n,k coefficients. Even though we do not claim that out findings provide the exact mathematical specification under which the t n,k coefficients factorially diverge as k → ∞, we believe that the presence (or absence) of the uniformity property of the asymptotic sequences {h n (a)} plays an essential role for the expansion functions to be consistent with the OPE. While it is obvious that expansion functions, where the nonuniformity is caused by pole-type singularities or other types of nonanalytic properties, are not suitable for phenomenological applications, the certainly surprising aspect is that even zeros, which may appear to have advantages for a rapid convergence at first sight, can have the same bad effect.

IV. GENERALIZATION TO FULL QCD
We conclude this article by considering the case of full QCD, where all terms of the QCD β-function are accounted for. Even though it appears quite unreasonable to argue that our observations concerning the asymptotic sequences {a n } and {H n,ℓ (a)} in the large-β 0 approximation may not apply in full QCD, it is worth having a closer look. However, in contrast to the large-β 0 approximation, we cannot obtain fully analytic results. We therefore rely on numerical studies and evidences rather than strict proofs. Still some closed formulas can be presented in the C-scheme for the strong coupling α s [20]. The C-scheme uses that only the one-and two-loop coefficients of the QCD β-function are renormalization-scheme invariant. It is therefore possible to adopt a scheme for α s where the QCD β-function takes the all-order form da(µ 2 ) d ln whereb 1 = β 1 /(2β 2 0 ) and β 1 = 102 − 38 n f /3. We furthermore demand that the C-scheme strong coupling differs from the common MS strong coupling by terms of order α 3 s and higher, which unambiguously fixes the C-scheme strong coupling to all orders. Furthermore, the QCD scale Λ QCD agrees with the one in the MS scheme. 8 Forb 1 = 0 the C-scheme agrees with the large-β 0 approximation. In the following we gather sufficient evidence that all qualitative insights obtained in the previous sections in the large-β 0 approximation are true as well in full QCD.

A. Spectral Function Moment Series in FOPT and CIPT
We start by considering the full QCD generalization of the spectral function moment analysis of Sec. II. The full QCD generalization (in the C-scheme) of the single pole Borel function in the large-β 0 approximation corresponding to an OPE matrix-element corrections of the form const. × ⟨O p ⟩/s p 0 in the Adler function reads B p,b 1 (u) = 1/(p − u) 1+2pb 1 (with p = 2, 3, . . .) and consistency with the OPE again demands that the resulting perturbative series of the spectral function moment δ 8 In Ref. [20] this strong coupling scheme was called the (C = 0)-scheme and it is numerically very close to the MS strong coupling, see also the Appendix of Ref. [11].
Let us first analyze the FOPT series coefficients d FOPT n,ℓ,p̸ =ℓ . We have derived them for many different choices for ℓ ̸ = p and find that the sequences |d FOPT n,ℓ,p̸ =ℓ | 1/n in n are indeed consistent with sequences that converge for all cases. We find that the limit superior approached by these series does not depend on the values of ℓ and p, but it does on the value ofb 1 . We also find that this limit superior agrees with the one obtained for the series of a ± ≡ a(−s 0 ± i0) in powers of a = a(s 0 ). In Fig. 8(a) we have displayed |d FOPT n,ℓ=0,p=2 | 1/n for orders up to n = 400 for different values ofb 1 . Carrying out a quadratic fit in 1/n we determined an accurate estimate for the limit superior which depends linearly onb 1 to very good approximation. We have displayed the outcome of this analysis for |d FOPT n,ℓ=0,p=2 | 1/n in Fig. 8(b) together with a fit function quantifying the coefficient ofb 1 . The outcome for a combined analysis using the coefficients d FOPT n,ℓ,p for ℓ = 0, 1, 2, 3, 4 and p = 1, 2, 3, 4 with ℓ ̸ = p and including also the series for a ± for orders 250 ≤ n ≤ 400 yields the result lim sup n→∞ d FOPT n,ℓ,p̸ =ℓ where the quoted uncertainty represents the values covered by the individual analyses. For b 1 = 0 we obtain π with very high accuracy, consistent with our analytic results in the large-β 0 approximation in Eq. (2.16). Our results reconfirm the convergence proof in Ref. [11] based on the renormalon calculus already mentioned at the end of Sec. II B and can be used to determine the radius of convergence of the FOPT series δ FOPT,ℓ [and the expansions of a ± in powers of a = a(s 0 )]. We have displayed the convergent FOPT moment series for the cases (ℓ, p) = (0, 2) and (1, 2) in Fig. 9 as the red dots. We note that we have also checked the FOPT series coefficients for physical weight functions with W (1) = 0. In analogy to our analytic findings for the large-β 0 approximation we found that there is a significant cancellation among the coefficients d FOPT n,ℓ,p̸ =ℓ , which does, however, not affect the limit superior. Also for physical weight functions the result given in Eq. (4.3) applies.
For n f = 3, which is relevant for hadronic τ decays, we haveb 1 = 32/81 ≃ 0.395 which yields lim sup n→∞ |d FOPT n,ℓ=0,p=2 | 1/n = 4.528 ± 0.020 giving α s (s 0 ) = 0.3083 ± 0.0013 as the convergence radius in the C-scheme. This corresponds to α s (s 0 ) = 0.3151 ± 0.0012 in the MS scheme which for s 0 = m 2 τ is right within the world average of the MS strong coupling at the τ lepton scale α s (m 2 τ ) = 0.312 ± 0.015. This interesting fact was already pointed out many years ago (albeit with much lower precision) by Pich and Le Diberder in Ref. [6]. So using FOPT may have issues as well. In practice, however, the behavior of the FOPT series at low orders is still quite good for all values close to the convergence radius due to the cancellation among the d FOPT n,ℓ,p coefficients for physical weight functions mentioned above. It is furthermore easy to check that the FOPT series behavior at intermediate orders is similar regardless whether α s (m 2 τ ) is chosen slightly above or below the convergence radius. At this point we remind the reader that here we discuss the behavior of FOPT for a particular contribution in the hadronic τ spectral function moment series. In full phenomenological applications, the FOPT series is asymptotic due to the unavoidable appearance of coefficients d FOPT n,ℓ,p for ℓ = p. Our finding, however, implies that high-precision determinations of the strong coupling from hadronic τ spectral function moments need to be interpreted with some care.
The CIPT expansion functions H n,ℓ (a) in full QCD have been determined analytically in Ref. [9] 9 as an infinite sum for positive real a. In the C-scheme the sum terminates due to the form of the β-function and the integral representation for the H n,ℓ (a) functions can be written as The integral reduces to Eq. (2.19) in the large-β 0 approximation, whenb 1 = 0. The results can be given in closed form in terms of a ± = a(−s 0 ± i0) and read H n≥2,0 (a) = 1 2πi 2b 1 (a n + − a n − ) n − a n−1 H n,ℓ>0 (a) = ℓ n+2ℓb 1 −1 a 2ℓb 1 e − ℓ a n Γ(n + 2ℓb 1 − 1) Since the first two terms in the expansion of the strong coupling at some scale in terms of powers of the strong coupling at another scale in full QCD and the large-β 0 approximation agree, the ratio of the H n,ℓ functions defined in Eq. (2.5) in full QCD to the ones in the large-β 0 approximation approach unity in the limit a → 0. It is therefore obvious that the {H n,ℓ (a)} also represent asymptotic sequences in full QCD. We also find that the H n,ℓ functions again exhibit zeros along the positive real a axis approaching zero as n → ∞. In Fig. 10 we display the positive zeros of the H n,ℓ functions in full QCD for a large number of ℓ and n values for illustration in analogy to the large-β 0 approximation shown in Fig. 4. This −i/π(−t 0 ) n ln(t + /t − ). So the result on the RHS given in Ref. [9] is correct, but it belongs to a different function than shown on the LHS. shows that the asymptotic sequences {H n,ℓ (a)} in full QCD are nonuniform as well. Using relation (A6) it is straightforward to derive the asymptotic expressions for the H n,ℓ>0 (a) as n → ∞ for positive real a yielding with θ n ≡ (2ℓb 1 + n) arctan[Im(a + )/Re(a + )] − ℓ Im(a + )/|a + | 2 . This allows us to also determine the limit superior of |H n,ℓ (a)| 1/n as n → ∞ in full QCD, which is surprisingly simple, lim sup n→∞ |H n,ℓ (a)| 1/n = |a + | . (4.10) As for the large-β 0 approximation, we can now see analytically that the CIPT series δ (4.11) We have displayed the divergent CIPT moment series for the cases (ℓ, p) = (0, 2) and (1,2) in Fig. 9 as the blue dots, showing, just like in the large-β 0 approximation, the apparent convergence of the CIPT series at intermediate orders and the discrepancy to the value of the convergent FOPT series. Finally, we also studied the large-order behavior of the s ℓ n,k coefficient for the series of the H n,ℓ functions in powers of a and the t ℓ n,k coefficients for series of a n in terms of the H n,ℓ functions, see Eqs. (2.24) and (3.1), respectively. Our findings are briefly summarized as follows: As for the large-β 0 approximation, the former series has a finite radius of convergence, while the latter does not converge for any value of a. Just like in the large-β 0 approximation, we find that the convergence radius arising from the s ℓ n,k coefficients for the expansion of the H n,ℓ functions in powers a n agrees with the convergence radius obtained from the FOPT series coefficients and the expansion of the a ± in a n powers, see Eq. of a yields a limit superior determined by the Taylor expansion of a ± . For the H n,ℓ>0 the same property follows from the fact that the combination of the asymptotic expansions of the two functions h n,ℓ in Eq. (4.7) yields a convergent series, see the comment in Footnote 3. Exemplarily, we have displayed the results for |s ℓ n,k | 1/k and |t ℓ n,k | 1/k as a function of k for different combinations of (n, ℓ,b 1 ) in Figs. 11(a) and 11(b), respectively. For the t ℓ n,k we again find that their sequence in k diverges factorially just as in the large-β 0 approximation for any ℓ. The series of a and a 2 in terms of the functions H n,ℓ (a) is shown for different values of (ℓ,b 1 ) in Figs. 12(a) and 12(b) exhibiting again the apparent convergence at intermediate orders, the discrepancy to the actual value and the eventual divergence. The results show a similar behavior for any other value of ℓ and the divergence also arises for physical weight functions.
Overall, we find beyond all reasonable doubt that all features of the CIPT series expansion using the H n,ℓ functions we found in the large-β 0 approximation are also present in full QCD.
In particular, any finite FOPT series is in general divergent when transformed into a CIPT expansion.

V. SUMMARY AND CONCLUSIONS
In this article we have provided a detailed discussion on the mathematical aspects of CIPT that has been used as a main method to make perturbative predictions for τ hadronic spectral function moments. In contrast to the FOPT expansion in powers of the strong coupling at a particular renormalization scale, [α s (s 0 )] n , CIPT yields expansions in functions where an integration over the renormalization scale is carried out. So CIPT provides series in nontrivial functions of the strong coupling that differ fundamentally from the solutions of the RGE equation for the strong coupling which are used in FOPT. While previous discussions on the same subject were predominantly based on explicit calculations of the series to a particular order [21][22][23] or the renormalon calculus, where different kinds of models for the Borel transform of the actual series were considered [7][8][9][10][24][25][26], we have carried out our analysis at a more basic level analysing directly the actual series using the elementary definitions and theorems on the (non)convergence of series and on asymptotic expansions.
For the most part of this article we have used the large-β 0 approximation, where all results can be obtained in an analytic way and all conclusions and considerations can be made in a mathematically rigorous way. We have reconfirmed that CIPT provides a well-defined asymptotic expansion, like the common power series expansion. However, the CIPT expansion functions represent nonuniform asymptotic sequences due to zeros along the positive axis that approach α s = 0 at large orders. We have proved that CIPT yields divergent series expansions without any region of convergence in cases where the OPE demands the existence of a region of convergence. This is in contrast to FOPT which indeed leads to series with a convergence region in these cases. We have found that the reason why the CIPT expansions frequently exhibits better series behavior than FOPT at low and intermediate orders is related to the zeros and because the CIPT expansion functions are bounded.
The most important finding of our analysis is that, while each CIPT expansion function has a finite region of convergence when written as a FOPT power series and represents an absolutely convergent series in this region, the inverse is not true. In other words, any FOPT power term [α s (s 0 )] n yields a divergent series when expressed as a sum of CIPT expansion functions regardless of the value of α s (s 0 ). There is a range of orders where the CIPT partial sum stabilizes before it diverges, but the truncated value of the CIPT series at these orders differs systematically from the value [α s (s 0 )] n . In other words, a convergent FOPT series will in general diverge in CIPT, and using CIPT will in general yield a degradation of the convergence properties of a series and an apparent convergence may yield an unphysical value inconsistent with the OPE. It is the latter property that makes phenomenological applications of CIPT dangerous. This does not exclude the use of CIPT as a phenomenological expansion method, but it must be used with great care. A phenomenologically consistent way to apply CIPT was suggested in Refs. [11,12].
Using models for expansion functions forming asymptotic sequences, we have shown that this property appears to be related to the zeros of the CIPT expansion functions which approach α s = 0 at high orders. It is kind of ironic that these zeros appear to be harmless at first sight and even beneficial as they suppress the numerical size of the CIPT expansion functions. We have provided numerical evidence that all our findings obtained in the large-β 0 approximation are also valid in full QCD. The primary use of our findings is that they can be used as a tool to also test the consistency of other expansion methods, where an integration over the renormalization scale is carried out.
Note Added: During submission of this work we received Ref. [27] where the concept of the "asymptotic separation" put forward in Refs. [9,10] was confirmed based on an alternative approach.
The root test specifies the property of absolute convergence which is sufficient for our discussions. Since we apply our considerations for series in the strong coupling which has uncertainties, the particular case L = 1 has no particular meaning for us, and will therefore not be discussed in any way. Since we do not only consider power series, but also series of functions we also introduce a generalization of the circle of convergence known for power series. For infinite series of functions ∞ i=1 f i (x) that converge to a function F (x) the property of uniform convergence is important as linear operations on F (x) can then typically be carried out by doing the operation on the individual functions f i (x) and summing afterwards.
Definition A.4 (Uniform Convergence). A series of functions ∞ i=1 f i (x) that converges to the function F (x) in the interval I is said to be uniformly convergent in an interval I ′ ⊆ I, if for any ϵ > 0, a number N (ϵ) independent of x exists such that | n i=1 f n (x) − F (x)| < ϵ for all n ≥ N (ϵ) and all x ∈ I ′ .
Weierstrass has formulated a simple bound criterion for uniform convergence. This theorem and the following corollary are important in our discussions as well.
Theorem A.5 (Weierstrass' Uniform Convergence Theorem). If each of the functions f i (x) is defined and bounded in the interval I, i.e., |f i (x)| ≤ γ i for all x ∈ I, and if the series of positive terms ∞ i=1 γ i converges, the series ∞ i=1 f i (x) converges uniformly in I. Corollary A.5.1. A power series ∞ i=1 c i x i with interval of convergence I converges uniformly in every subinterval I ′ ⊂ I.
We finally quote the powerful and important Weierstrass theorem on the reordering of double series without referring to the absolute convergence of the double series.
Theorem A.6 (Weierstrass' Double Series Theorem). Consider an infinite set of functions f i (x) that are analytic for |x| < r, so that the power expansions f i (x) = ∞ k=0 a (i) k x k exist and converge at least for |x| < r for all i. Furthermore, consider a convergent series of these functions F (x) = ∞ i=0 f i (x) that is uniformly convergent for |x| ≤ ρ for every ρ < r, so that the series converges in particular everywhere within the interval |x| < r and defines the function F (x) there. Then, the infinite sums A k = ∞ i=0 a (i) k are convergent and the infinite sum ∞ k=0 A k x k converges to F (x) for |x| < r, so that k )x k and is analytic for |x| < r.

Order Symbols, Asymptotic Sequences and Asymptotic Expansions
In gauge quantum field theories perturbative series are typically not convergent, but only asymptotic. We therefore also specify the basis of asymptotic expansions collecting the relevant mathematical definitions. In the following, R is a domain in the complex plane and x 0 a point inR, the closure of R. The functions f (x), g(x), and the sequence of functions ϕ 1 (x), ϕ 2 (x), ϕ 3 (x), . . ., are analytic in R. We furthermore abbreviate the latter sequence of functions as {ϕ n (x)}. For the purpose of our considerations we always assume that {ϕ n (x)} is an infinite set so that n runs over all natural numbers.
We start by quoting the definition of the well known O-relation and the lesser known "little" o-relation of two functions. While the O-relation is a statement on two functions approaching a certain point at the same type of "speed" (e.g. linear or quadratic), the "little" o-relation states that one of the two function approaches that point much faster. It is important for the definition of asymptotic expansions.
Definition A.7 (O-relation). It is said that f = O(g) as x → x 0 if there exists a constant A > 0 and a neighborhood U of x 0 so that |f (x)| ≤ A|g(x)| for all x ∈ R ∩ U .
Definition A.8 (o-relation). It is said that f = o(g) as x → x 0 , if for any ϵ > 0 there exists a neighborhood U ϵ of x 0 such that |f (x)| < ϵ|g(x)| for all x ∈ R ∩ U ϵ . When g(x ̸ = x 0 ) ̸ = 0 in some neighborhood of x 0 , the condition is equivalent to lim x→x 0 f (x)/g(x) = 0.
Asymptotic expansions approximate a function as a sum of other functions, which in order to be well-defined must form asymptotic sequences [29]. The latter are specified by the following definitions, where the special property of uniformity plays an important role in this article.
Definition A.10 (Uniform Asymptotic Sequence). A sequence of functions {ϕ n (x)} is an uniform asymptotic sequence as x → x 0 in R if ϕ n+1 = o(ϕ n ) uniformly in n. This means that for any ϵ > 0 there exists a neighborhood U ϵ of x 0 such that |ϕ n+1 (x)| < ϵ|ϕ n (x)| for all x ∈ R ∩ U ϵ and all n.
The requirement of an asymptotic sequence to be uniform excludes e.g. that the functions ϕ n have singularities, but also zeros which approach x 0 for increasing n. The sequence of power terms {x n } is a uniform asymptotic sequence as x → 0 in R as one can easily check from the property that x n+1 /x n = x independently of n.
We are now ready to state the definition of an asymptotic expansion, where it is at this point not relevant whether the asymptotic sequence involved is uniform in n or not. Corollary A.11.1. The coefficients of the previously defined asymptotic expansion of f (x) as x → x 0 with respect to the asymptotic sequence {ϕ n (x)} can be determined through the following recurrence formula: c n = lim x→x 0 [f (x) − n−1 i=1 c i ϕ i (x)]/ϕ n (x). Furthermore, this implies that the coefficients are unique.
where the term cos(nπ/2) accounts for B n being zero for large odd n and the correct overall sign. We also quote Stirling's formula Γ(n + 1) = n! n→∞ ≍ √ 2πn(n/e) n , which is used several times.