The strong coupling from an improved τ vector isovector spectral function

We combine ALEPH and OPAL results for the spectral distributions measured in τ → π − π 0 ν τ , τ → 2 π − π + π 0 ν τ and τ → π − 3 π 0 ν τ decays with (i) recent BaBar results for the analogous τ → K − K 0 ν τ distribution and (ii) estimates of the contributions from other hadronic τ -decay modes obtained using CVC and electroproduction data, to obtain a new and more precise non-strange, inclusive vector, isovector spectral function. The BaBar K − K 0 and CVC/electroproduction results provide us with alternate, entirely data-based input for the contributions of all exclusive modes for which ALEPH and OPAL employed Monte-Carlo-based estimates. We use the resulting spectral function to determine α s ( m τ ), the strong coupling at the τ mass scale, employing ﬁnite energy sum rules. Using the ﬁxed-order perturbation theory (FOPT) prescription, we ﬁnd α s ( m τ ) = 0 . 3077 ± 0 . 0075, which corresponds to the ﬁve-ﬂavor result α s ( M Z ) = 0 . 1171 ± 0 . 0010 at the Z mass. While we also provide an estimate using contour-improved perturbation theory (CIPT), we point out that the FOPT prescription is to be preferred for

I.

INTRODUCTION
Since the calculation of the term of order α 4 s in the Adler function [1], there has been a revived interest in the determination of the strong coupling, α s (m τ ), at the τ mass scale m τ , from non-strange hadronic τ decays. Two LEP experiments, ALEPH [2][3][4] and OPAL [5] conducted measurements of hadronic τ decays from which the inclusive non-strange vector (V ) and axial (A) isovector spectral functions were extracted with high accuracy as a function of the squared invariant mass s. A third experiment, CLEO, also used V and A inclusive spectral data in an early determination of α s (m τ ) [6], but these data have not been made publicly available. 1 Most determinations of α s since Ref. [1] have been based on the ALEPH data, the most recent of these using the 2013 version of this data [4,8,9], in which an earlier problem [10] with the data covariance matrix was corrected [4]. An exception is the determination of α s in Ref. [11], which was based on the OPAL data. This was later updated in Ref. [12] for changes in the exclusive-mode τ branching fractions (BFs) since the original 1998 OPAL publication. The determinations based on the ALEPH [8] and OPAL [12] data are consistent with each other, leading Ref. [8] to quote a weighted average as the best result taking the determination from both data sets into account. While this weighted average should be reliable, it was not based on a fit of the combined ALEPH and OPAL data, and compatibility of the two α s values does not test the compatibility of the two data sets directly. 2 Clearly, what one would like to do instead is to combine the two data sets to produce a single data set whose spectral functions and corresponding covariance matrices reflect locally, i.e., in an s-dependent manner, the combined constraints of the ALEPH and OPAL data. The process of combining the data sets tests for their compatibility, and the result is a data set with smaller errors than either of the two data sets alone. In this sense, the combined spectral functions will be the "best available" extracted from hadronic τ decays. A number of hadronic quantities useful for hadron phenomenology, such as α s , certain low-energies constants of chiral perturbation theory and operator product expansion condensates can then be determined from the V and A spectral functions.
In this paper, we begin this program by constructing the combined inclusive non-strange V spectral function, and using this to obtain the most precise value of the strong coupling that can be obtained from the combined ALEPH and OPAL hadronic τ -decay data in the V channel, supplemented with e + e − → hadrons cross-section data for some small but nonnegligible residual exclusive modes. The process of combining the two spectral functions involves several steps, and includes updating the normalizations of exclusive channels using updated BFs, before the two data sets are combined.
There are two reasons for limiting ourselves to the V spectral function in this paper. The first is that the V channel is dominated by the 2π and 4π decay modes, while the remaining channels, including those for which ALEPH and OPAL used Monte-Carlo (MC) input, play a much smaller relative role in the V than in the A channel. For OPAL, all V -channel modes apart from the dominant 2π and 4π channels, π − π 0 , π − 3π 0 and 2π − π + π 0 , are listed as having a MC source [5]. For ALEPH, MC simulations are used for the K − K 0 , KK2π and 6π V contributions with the simulation also used for the part of the ωπ − distribution not reconstructed in the 2π − π + π 0 mode [2][3][4].
A second important reason for focusing on the V channel is that the CVC (conserved vector current) relation 3 between the τ -based V spectral function and I = 1 hadronic e + e − cross section contributions allows almost all of the smaller, but still numerically relevant, contributions from exclusive modes other than the dominant 2π and 4π ones to be significantly improved using recent high-precision exclusive-mode e + e − → hadrons cross-section data. The use of CVC and electroproduction data paves the way for an almost fully experimental, and in this sense improved, determination of the V spectral function. Such CVC improvements are, of course, not possible in the A channel, which, in addition, receives larger relative contributions from higher-multiplicity exclusive modes for which exclusive-mode spectral function contributions and covariances were not provided by OPAL and ALEPH.
For these reasons, we postpone a discussion of the A case to a future work. As far as the determination of α s is concerned, we found, in previous work, that, while the addition of the A channel to the analysis provided a nice consistency check [8,11,12], it did not help reduce the error on α s . We thus consider the determination of α s using an updated version of the V spectral function alone to be of interest.
There are thus two parts to the work reported in this paper. In the first part, we update the determination of the inclusive non-strange V spectral function. This itself involves two steps. First, we combine, and hence update, the results for the contributions from the exclusive modes for which ALEPH and OPAL data is publicly available (the dominant 2π and 4π modes), using the method employed to combine exclusive-mode e + e − → hadrons R(s) data from different experiments and described in detail in Ref. [13]. Second, we use recent results for τ → K − K 0 ν τ , together with CVC and recent exclusive-mode e + e − → hadrons cross section results, to improve the determination of the contributions from the remaining modes, which in this paper we will refer to as the "residual" modes. In the second part of the paper, we apply the strategy developed in Refs. [8,11,12] to extract α s (m τ ) from the improved inclusive V non-strange spectral function obtained in the first part.
The determination of α s from V and/or A spectral functions makes use of finite energy sum rules (FESRs), which allow us to relate α s at the τ scale through the Adler function, calculated in QCD perturbation theory, to integrals of the spectral function from threshold to the τ mass [14][15][16][17][18][19][20][21][22]. Even so, with the strong coupling at scales around the τ mass being rather large, these sum rules are "contaminated" by non-perturbative effects. These non-perturbative effects are clearly visible in the experimental data, since the shape of the experimental inclusive spectral function does not agree with perturbation theory. As a result of the presence of resonances, the experimental spectral functions oscillate around the perturbative predictions, with these oscillations remaining visible for s close to m 2 τ . Several methods have been designed for dealing with these non-perturbative effects. The oldest method (the "truncated OPE" strategy) employs weight functions assumed to suppress the effects of the observed "duality violating" resonant oscillations and assumes the reliability of a truncation in dimension of operator product expansion (OPE) contributions required to make the analysis practical [22,23]. A newer method, the "DV-model strategy," instead takes quark-hadron duality violations (DVs), i.e., collective resonance effects, into account, and requires only very mild assumptions about the behavior of the OPE [11,24]. While it is not straightforward to ascertain the reliability of estimates of non-perturbative effects, self-consistency tests show that the truncated OPE strategy leads to unreliable results, with a theoretical uncertainty arising from the neglect of higher-order OPE terms and of II.
THEORY OVERVIEW In Sec. II A, we briefly recapitulate the use of FESRs to extract α s from spectral function input and define our theoretical framework. The choice of sum rules employed in our fits is discussed in Sec. II B. In Sec. II C we argue that fixed-order perturbation theory (FOPT) should be favored over contour-improved perturbation theory (CIPT) [28,29], if the goal is to convert our result for α s (m τ ) to a value at the Z mass to be compared to MS values of α s (M Z ) obtained from other sources. We also explain how we estimate the systematic error associated with the necessary truncation of perturbation theory beyond order α 4 s .

A. Finite energy sum rules
The sum-rule analysis starts from the current-current correlation functions where J µ stands for the non-strange vector (V ) current uγ µ d or axial (A) current uγ µ γ 5 d, and the superscripts (0) and (1) label spin. The decomposition in the third line is useful because Π (1+0) (q 2 ) and q 2 Π (0) (q 2 ) are free of kinematic singularities. With s = q 2 , the spectral function and the known analytical properties of Π (1+0) (z), application of Cauchy's theorem to the contour in Fig. 1 implies the FESR The sum rule is valid for any s 0 > 0 and any weight w(s) analytic inside and on the contour [14][15][16][17][18][19][20]. In this paper, we will always choose w(z) to be polynomial in z.
There is a cut on the positive real axis starting at s = q 2 = 4m 2 π (a pole at s = q 2 = m 2 π and a cut starting at s = 9m 2 π ) for the V (A) case. The solid curve shows the contour used in Eq. (2.3).
The flavor-ud V and A spectral functions can be experimentally determined from the differential versions of the ratios of the hadronic decay width induced by the relevant current to that for the electron mode. Explicitly [30], where S EW is a short-distance electroweak correction and With the exception of the known pion-pole part of ρ A;ud , the spectral functions ρ V /A;ex (s 0 ) and the theoretical representation of the contour integral on the right-hand side by I where the logarithmic z dependence of the OPE coefficients C 2k can be calculated in perturbation theory. For the k = 0 term, it is convenient to consider, instead of Π(z), the Adler function D(z) ≡ −z dΠ(z)/dz, which is finite and independent of the renormalization scale µ. Accordingly, the k = 0 contribution to the right-hand side of Eq. (2.3) can be expressed in terms of the Adler function via partial integration. The k = 0 contribution D 0 (z) to D(z) takes the form where the coefficients c nm are known to order α 4 s [1]. The independence of D 0 (z) on µ implies that only the coefficients c n1 are independent; the c nm with m > 1 can be expressed in terms of c n1 through use of the renormalization group. 5 In the MS scheme, c 01 = c 11 = 1, c 21 = 1.63982, c 31 = 6.37101 and c 41 = 49.07570, for three flavors [1]. While c 51 is not known at present, we will use the estimate c 51 = 283 of Ref. [32], with a 50% uncertainty. For the running of α s we use the four-loop MS β-function, but we have checked that using five-loop running instead [33,34] leads to differences of order 10 −4 or less in our results for α s (m τ ).
The C 2k with k ≥ 1 are different for the V and A channels, and, for k > 1, contain nonperturbative D = 2k condensate contributions. As in Refs. [11,12], we will neglect purely perturbative quark-mass contributions to C 2k , k ≥ 1, as they are numerically very small for the non-strange FERSs we consider in this paper. We will also neglect the z-dependence of the coefficients C 2k for k > 1. For a more detailed discussion of our treatment of the D > 0 OPE contributions, we refer to Ref. [11].
Perturbation theory, and in general the OPE, breaks down near the positive real s = q 2 axis [35]. If this were not the case, Eq. (2.3) would establish a direct correspondence between the OPE and the resonant behavior of the experimental spectral function, generally referred to as quark-hadron duality. We account for the breakdown of this duality by replacing the right-hand side of Eq. (2.3) by The imaginary parts 1 π Im ∆ V /A (s) can be interpreted as the duality-violating parts ρ DV V /A (s) of the V /A spectral functions, and represent the resonance-induced, oscillatory parts of the spectral functions not captured by the OPE.
In Ref. [24], we developed a theoretical framework for quark-hadron duality violations in terms of a generalized Borel-Laplace transform of Π(q 2 ) and hyperasymptotics, building on earlier work [37][38][39][40]. In the chiral limit, and assuming that for high energies the spectrum 5 See for instance Ref. [31].
6 becomes Regge-like in the N c → ∞ limit, we showed that the asymptotic form of ρ DV V /A (s) for large s can be parametrized as (2.12) up to slowly varying logarithmic corrections in the argument of the sine factor, and with γ ∼ 1/N c small but non-zero. 6 The parameters β V /A are directly related to the Regge slopes in the V and A channels, and the parameters γ V /A to the (asymptotic) ratio of the width and the mass of the resonances in those channels. While the framework of Ref. [24] is rather general, and the derivation of Eq. (2.12) is based on generally accepted conjectures about QCD (primarily Regge behavior), it does not provide a derivation from first principles. This introduces a certain model dependence in our analysis which, however, can be tested by fits to the data. Such tests, in particular, will provide information about the values of s above which this asymptotic form is likely to be sufficiently accurate. We emphasize, however, that modifications to the parametrization of Eq. (2.12) are constrained by the general framework of Ref. [24]. Equation (2.12) introduces, in addition to α s and the D ≥ 4 OPE condensates, four new parameters in each channel: δ V /A , γ V /A , α V /A and β V /A . This can be compared to the truncated-OPE approach in which DVs are neglected. Since resonance-induced oscillations are clearly visible in the experimental spectral data, and their dynamical effect is comparable in size to the α s -dependent dynamical effect of perturbative corrections to the (α sindependent) parton model contribution [25], this approach also assumes a model: one in which the parameters δ V /A are effectively set to infinity by hand. This is a stronger assumption, and one that has been shown to fail a number of subsequent, more stringent data-based tests [25,26].
In summary, as in Refs. [11,12], we will assume that Eq. (2.12) holds for s ≥ s min , with s min to be determined from fits to the data. This assumes of course that the s min for which this assumption is valid includes a region below m 2 τ , i.e., that both the OPE (2.7) and the DV parametrization (2.12) can be used in some interval below m 2 τ .

B. Choice of weight functions and the OPE
The logarithmic s dependence of the OPE coefficients C D (s) is calculable in perturbation theory. Because the running of α s becomes visible only at non-leading order in α s , this s dependence is an O(α 2 s ) effect in the chiral limit. Such effects were found to be safely negligible for D > 0 in the sum-rule analysis of the OPAL data reported in Ref. [11], and we will thus ignore them for D > 0 in the present analysis as well. With this simplification, a term in a weight w(z) proportional to the monomial z n picks out the D > 0 OPE contribution with D = 2k = 2(n + 1) in the sum rule (2.3). 7 The choice of a polynomial weight w(z) thus projects the sum rule on a finite number of D > 0 terms in the OPE.
In this paper, we will consider the weights w(z) = w n (z/s 0 ) with w 0 (y) = 1 , (2.13) w 2 (y) = 1 − y 2 , w 3 (y) = (1 − y) 2 (1 + 2y) , w 4 (y) = (1 − y 2 ) 2 , 7 where the subscript indicates the degree of the polynomial. These weights explore OPE terms with D ≤ 10, and form a linearly independent basis for polynomials up to degree four without a linear term. The weight w 0 (y) projects only the D = 0 term of the OPE (i.e., pure perturbation theory), the weight w 2 (y) projects, in addition, D = 6, w 3 (y) projects D = 0, D = 6 and D = 8, while w 4 (y) projects D = 0, D = 6 and D = 10. As the OPE itself diverges as an expansion in 1/z, it is safer to include sum rules with low-degree weights such as w 0 (y) and w 2 (y) in the analysis, and check for consistency among sum rules with different weights. We note that w 3 (s/s 0 ) = w T (s; s 0 ), cf. Eq. (2.6).
None of these weights contain a term linear in z, and thus the D = 4 OPE term does not contribute to the sum rules with these weights. This choice is motivated by the results of Ref. [43], in which a renormalon-model-based study suggested that perturbation theory is particularly unstable for sum rules with weights containing such a linear term. 8 The weights w 2,3,4 (y) are "pinched," i.e., they have zeroes at z = s 0 , and thus suppress contributions from the region near the timelike point z = s 0 on the contour, and hence also the relative importance of integrated DV contributions [45,46]. The weight w 2 (y) has a single zero at z = s 0 (a single pinch), while the weights w 3 (y) and w 4 (y) are doubly pinched, i.e., have a double zero at z = s 0 .

C. Perturbative uncertainties and FOPT vs. CIPT
It has become common practice to consider different resummations of perturbation theory in order to obtain insight into the effect of neglecting terms beyond those explicitly included in evaluating the D = 0 (i.e., perturbative) contribution to the right-hand side of Eq. (2.3). The two most commonly used resummation prescriptions are fixed-order perturbation theory (FOPT), in which the scale µ in Eq. (2.8) is chosen to be fixed at µ 2 = s 0 , and contourimproved perturbation theory (CIPT) [28,29], a partial resummation obtained by choosing µ 2 equal to −z point by point along the contour. In the CIPT prescription, the coupling is run along the contour |z| = s 0 , using the four-(or five-)loop beta function; as a result, only terms with m = 1 survive in Eq. (2.8). The two prescriptions lead to significantly different values of α s , with the difference being comparable to the combination of all other errors [8,12], but more significant, since the results obtained with different prescriptions using the same data are highly correlated.
The choice of such a prescription is entangled with the choice of renormalization scheme, since the choice of scheme also affects higher orders in perturbation theory. Both FOPT and CIPT prescriptions are usually considered to be MS schemes, but clearly, if the choice between FOPT and CIPT is rephrased as a choice of scheme, these two schemes are different. Since α s is not a physically measurable quantity, ideally, one would like to choose a scheme, and a prescription, which corresponds to the scheme chosen to quote other determinations of α s (such as that from Z decay itself), so that a direct comparison is possible at M Z .
In our case, the experimental quantities from which we determine α s are the spectral integrals I (w) ex (s 0 ), in which s 0 is varied between s 0 = s min and s 0 = m 2 τ . Apart from DVs, which appear in a transseries beyond the OPE, these quantities are then expressed in terms of the OPE, which at D = 0 is parametrized by α s (µ) and the ratio of scales s 0 /µ 2 , while at D ≥ 4 also the condensates enter. This leads to the theoretical representations I (w) th (s 0 , s 0 /µ 2 , α s (µ)), which are then to be equated with I (w) ex (s 0 ). Since I (w) ex (s 0 ) is an observable with a single physical scale, s 0 , the natural choice is to choose the scale µ equal to the scale of the observable, i.e., to choose µ 2 = s 0 . 9 This corresponds directly to the way the scale is chosen for other determinations of α s . For instance, the hadronic Z-decay rate can be expressed perturbatively in terms of α s (µ) in the MS scheme, and in that case, the natural choice of scale is µ = M Z . This leads us to conclude that, in the case of hadronic τ decays, the prescription most directly comparable with other determinations is FOPT. We emphasize that we do not claim to know which prescription, at a given order, gives the best approximation to the QCD answer for I (w) ex (s 0 ). This may depend on the weight w, and on the order in perturbation theory [43]. The point is to choose a scheme that corresponds most closely to the scheme employed in other determinations of α s at the Z-mass. Of course, it is necessary to estimate the systematic uncertainty inherent in the truncation of perturbation theory, and we will return to this point at the end of this subsection.
Before we do this, we present an additional reason for using FOPT as the prescription to be used in the FESR determination of α s from hadronic τ decays. It is well known that perturbation theory for the Adler function, and for the quantities I (w) th (s 0 ), is not convergent. Attempts to go beyond perturbation theory using Borel resummation techniques lead to ambiguities in the Borel sum that necessitate the introduction of D > 0 terms in the OPE [47]. The most well-known example is that of the leading infrared renormalon in massless QCD. This leads to the D = 2k = 4 term in Eq. (2.7), which removes the ambiguity associated with this renormalon. The appearance of D > 0 terms in Eq. (2.7) is thus intricately connected to resummations of perturbation theory.
In Ref. [48], it was shown that, in general, the FOPT and CIPT series lead to different Borel sums, with different analytical properties. Moreover, it was pointed out that this different analytical behavior of the Borel sums for I (w) th (s 0 ) appears to invalidate the correspondence between infrared renormalons and the D > 0 terms in the OPE in the case of CIPT. The Borel sum for the CIPT series does not allow for the usual renormalon ambiguities that in the case of FOPT are in one-to-one correspondence with terms in the OPE. Thus, there is a mismatch between the use of CIPT, and the representation of Π (1+0) (z) by the OPE, Eq. (2.7). This mismatch between the Borel sum and the OPE does not happen in FOPT.
This observation casts strong doubts on the consistency of using the OPE (2.7) in the case of CIPT. Again, this does not settle the issue of whether the (Borel sum of) FOPT or CIPT provides a better approximation to QCD. But it does imply that it is theoretically inconsistent to apply the OPE in the form (2.7) if one uses CIPT in the evaluation of the D = 0 perturbative contributions, and therefore casts doubts on all CIPT-based extractions of α s . We are thus led to the conclusion that FOPT should be taken as the preferred choice for analyzing hadronic τ decays using the OPE. While we will provide determinations of α s (m τ ) employing both FOPT and CIPT (ignoring the OPE subtleties in the case of CIPT) we will quote the FOPT value as our final result, to be compared with other MS determinations of α s . We will also give a CIPT value based on fits to I (w) ex (s 0 ), allowing the reader to compare to earlier CIPT results and assess the impact of changes in the input inclusive V spectral function.
We will determine the systematic error associated with the use and truncation of perturbation theory for our FOPT value of α s (m τ ) as follows. First, as already indicated in Sec. II A, we will vary the estimated value for c 51 = 283 by plus or minus 50%. This interval for c 51 includes all estimates available in the literature by a rather wide margin. The variation by ±50% was proposed in Ref. [32], and the estimate of Ref. [1] falls well inside this interval. Subsequent estimates of both the central value and the uncertainty of c 51 are also generously covered by the ±50% variation [49,50]. This choice of range, of course, provides an estimate only for the impact of the uncertainty in the value of the unknown coefficient c 51 , and might constitute an underestimate of the total perturbative error. At present, the perturbative expansion for the Adler function has been calculated to order α 4 s . An alternate, potentially more conservative, estimate of the perturbative error can thus also be obtained by omitting the O(α 4 s ) contribution altogether, i.e., by setting c 4m = 0 (as well as c nm = 0 for all n > 4) in Eq. (2.8). Finally, since our FOPT determination employs a range of s 0 varying between s min and m 2 τ , a third sensible estimate of the perturbative error can be obtained by considering, instead of just µ 2 = s 0 , also the alternate choices µ 2 = s min , µ 2 = m 2 τ and µ 2 = 2s 0 for the scale µ. We will take the largest of the variations in α s obtained by applying all three methods above as our best estimate for the systematic error associated with the necessary truncation of perturbation theory. We do not combine the errors obtained by using these three methods, as this would correspond to a double-counting of the estimated perturbative uncertainties. We will, however, add an independent measure of the uncertainty associated with the use of perturbation theory based on the comparison of the central values obtained from fits with different weights, as explained in more detail in Sec. IV C.

III. DATA
In this section, we construct an updated version of the inclusive, non-strange V spectral function, ρ ud;V (s), using publicly available ALEPH [2][3][4]51] and OPAL [5] τ data for the contributions of the dominant 2π and 4π exclusive modes, recent BaBar τ -decay results [52] for the contribution of the K − K 0 mode, and e + e − → hadrons cross section data as input to CVC evaluations of the contributions of the remaining exclusive modes.
In Refs. [2][3][4][5], the inclusive V and A invariant-mass-squared distributions were constructed as sums of (i) the measured (and publicly available) distributions for the main exclusive modes in the channel and (ii) the sum of small contributions from the remaining "residual" exclusive modes. The publicly available exclusive-mode distributions are normalized to then-current exclusive-mode BFs. While the accompanying inclusive-sum correlation matrices include the contributions from then-current exclusive-mode BF uncertainties and correlations for the main exclusive modes, the exclusive-mode correlation matrices are provided with the BF-uncertainty-induced contributions omitted, allowing subsequently improved BF information for these modes to be incorporated at a later time.
The V channel modes for which ALEPH and OPAL provide exclusive-mode distribution and correlation information are π − π 0 , π − 3π 0 and π − π + π − π 0 . For the A channel, ALEPH provides distribution and correlation information for only the two 3π modes, 2π − π + and π − 2π 0 , while OPAL provides this information, in addition, for one of the three 5π modes, 2π − π + 2π 0 . While all exclusive-mode distributions are normalized to then-current values of the corresponding exclusive-mode BFs, the s-dependences of some of the residualmode distributions are taken from Monte Carlo. For OPAL, this is true for all but the π − π 0 ω(→ non-3π) residual A mode. For ALEPH, the use of the Monte Carlo simulations is explicitly identified as entering the 6π, K − K 0 , KKππ and π − ω(→ non-3π) contributions to the V channel and the KKπ and KKππ contributions to the A channel. Since the information made publicly available by ALEPH and OPAL does not include the individual BF-normalized residual exclusive-mode distributions used by the collaborations in determining their final inclusive-sum results, it is not possible to update those residual-mode contributions to reflect subsequent improvements in our knowledge of the exclusive-mode BFs and/or new information on the s-dependence of exclusive-mode distributions for which Monte Carlo was previously employed. Improvements to the residual-mode contributions must, therefore, come from other sources. The dominant V and A channel modes (those for which both ALEPH and OPAL exclusive-mode distributions are available) represent 98.0% of the inclusive V channel BF and 94.2% of the continuum inclusive A channel BF.
Improving the treatment of residual-mode contributions is much easier for the V channel than for the A channel. The reason is that, strongly motivated by the drive to improve the determination of the Standard-Model hadronic vacuum polarization contribution to the anomalous magnetic moment of the muon, there has been an intensive program of collider and B-factory experiments aimed at determining the e + e − → hadrons cross-sections for all exclusive modes contributing to the inclusive R(s) ratio in the region below s 4 GeV 2 . A sizeable fraction of these exclusive modes can be uniquely classified as either I = 0 or I = 1 using G-parity. The CVC relation between the bare cross section for the electroproduction of the neutral member, , and the contribution of the charged isospin partner, X − , to ρ ud;V (s), then allows the cross section results for those exclusive residual modes with I = 1 to be used to determine the corresponding exclusive residual-mode contributions to ρ ud;V (s). Eq. (3.1) is valid up to isospin-breaking (IB) corrections which, in the absence of narrow interfering resonances, should be of the order a percent or so, and hence numerically negligible on the scale of the experimental errors on the already small residual-mode contributions. This strategy, of using CVC to improve the determination of an otherwise poorly determined V exclusive-mode τ spectral function contribution, was pioneered by ALEPH [51], which used the BaBar Dalitz-plot-analysis separation of I = 0 and I = 1 contributions to the e + e − → KKπ cross sections [53] to determine the V part of the τ → KKπν τ distribution, and hence the separation of that distribution into its V and A components. The CVC relation allows us to dramatically improve the vast majority of the residual-mode contributions to the V spectral function. This is especially helpful in the case of contributions from highermultiplicity modes, whose τ -decay distributions lie at higher s, increasingly close to the τ kinematic endpoint, and with, as a result, increasingly reduced statistical precision.
The main V residual mode for which such a CVC improvement is not possible is K − K 0 , where the e + e − → KK cross sections contain both I = 0 and I = 1 contributions, and it is not possible to identify only the I = 1 component. Fortunately, for this channel, BaBar [52] has recently published a rather precise determination of the unit-normalized τ → K − K S ν τ number distribution, allowing the residual-mode K − K 0 contribution to ρ ud;V (s), to be determined directly, without the use of CVC.
Using CVC and the recent BaBar τ K − K 0 results, 99.95% by BF of the inclusive V spectral function can be determined directly from experiment. The remaining 0.05% represents only 2.4% by BF of the already small sum of residual-mode contributions. CVC improvements are, of course, impossible for A channel residual-mode distributions. This, and the larger relative role played by residual-mode contributions in the A channel, are the primary reasons for our focus on the V channel in this paper.
The rest of this section is organized as follows. First, in Sec. III A, we specify the sources of external input employed in our update of ρ ud;V (s). Next, in Sec. III B, we outline the procedure used for combining the publicly available data from ALEPH and OPAL for the dominant 2π and 4π exclusive-modes, following closely that described in Ref. [13] for combining exclusive-mode e + e − cross sections from different experiments. Details of our updates of the individual residual exclusive-mode contributions are provided in Sec. III C. Finally, the resulting updated version of ρ ud;V (s), is presented in Sec. III D.
ALEPH quotes results in the form of exclusive-mode, BF-normalized contributions, dB X (s)/ds, to the differential BF distribution dB(s)/ds. The corresponding contributions to where B X is the BF for exclusive mode X, 1 is the corresponding experimental unit-normalized number distribution, and B e is the τ − → e − ν τνe BF. ALEPH results for dB X (s)/ds are updated by rescaling to the current value of B X . Current values of the external parameters, B e , S EW , V ud and m τ are then used to obtain the updated results for ρ X ud;V (s). OPAL quotes results in the form of the exclusive-mode contributions ρ X ud;V (s). These were obtained from the experimentally measured unit-normalized number distributions via Eq. (3.2), using then-current values of the exclusive-mode BFs, B X , and the external inputs B e , S EW , V ud and m τ . The underlying unit-normalized distributions are reconstituted using the values for the exclusive-mode BFs and external inputs quoted by OPAL, and converted to equivalent updated versions of the ρ X ud;V (s) using current values for these inputs. We employ the following values for the external parameters appearing in Eq. (3.2): for B e , the lepton-universality-improved HFLAV 2019 [54] result B e = 0.17814 (22); for m τ and |V ud |, the PDG 2020 [55] results m τ = 1.77686 (12) GeV and |V ud | = 0.97370 (14); and, for S EW , the result S EW = 1.0201(3) [56].
Specifics of the experimental inputs used in the determination of the residual-mode contributions to ρ ud;V (s) are detailed in Sec. III C below.

B. Combining the ALEPH and OPAL 2π and 4π data
We begin with the updated versions of the ALEPH and OPAL exclusive-mode distributions, ρ X ud;V (s), with X = π − π 0 , π − 3π 0 and π − π + π − π 0 , obtained as outlined in the previous subsection. Ideally, we would like to combine the results for each of these exclusive modes separately, first combining the ALEPH and OPAL unit-normalized number distributions (which are independent of the BFs) and then multiplying the resulting combined exclusive distributions by the corresponding BFs. It turns out, however, that this is not possible. The reason is that the correlation matrices for the π − 3π 0 and π − π + π − π 0 distributions for both experiments have zero eigenvalues, i.e., 100% correlations between different bins, and hence are not invertible. This prevents us from combining the ALEPH and OPAL data for the individual 4π modes in the manner described below. If, however, we first sum the contributions from all three modes, we find that the correlation matrices for the resulting three-mode-sums are well behaved for both ALEPH and OPAL. We thus combine the ALEPH and OPAL exclusive-mode results by first summing, for each experiment separately, the contribution to ρ ud;V (s) and the corresponding covariance matrices from π − π 0 , π − 3π 0 and π − π + π − π 0 , using updated versions of the exclusive-mode BFs, and then combining those results using the method outlined below.
Following Ref. [13], we choose a number of clusters, distributed over the interval 0 < s ≤ m 2 τ . We assign a number of consecutive ALEPH and OPAL data points to each cluster m, m = 1, . . . , N cl , with N cl the total number of clusters; N m will be the total number of data points in cluster m. If the collective ALEPH and OPAL data points are parametrized by pairs (s i , d i ), where d i is the ALEPH or OPAL data point for the spectral function assigned to the s-value s i , we define weighted cluster averages where the sum is over all data points in cluster m and σ 2 i is the variance of d i , i.e., the σ 2 i are the diagonal elements of the covariance matrix C ij for the spectral-function data points d i . The set of s (m) then constitutes the values of s at which the combined spectral function ρ (m) will be defined. The values of ρ (m) will be determined by linear interpolation, minimizing where N = N cl m=1 N m is the total number of (ALEPH and OPAL) data points, and the piece-wise linear function R(s; ρ) is defined by where ρ is the vector of fit parameters ρ (m) , m = 1, . . . , N cl . At the boundaries, we extrapolate: Minimizing χ 2 (ρ) yields the linear equations which can be solved for the ρ (m) , with the cluster covariance matrix C mn given by The procedure for combining exclusive spectral functions followed in Ref. [13] is more complicated than the one outlined above. The inclusion of uncertainties in the BFs in the covariance matrices can lead to a bias in the fit [57], and the method employed in Ref. [13] adjusts for this bias [58]. However, for this to work, we would need to combine each channel separately, because multiplication by the exclusive-mode BF is needed in each channel to turn the normalized distribution into the corresponding contribution to ρ ud;V (s). This path is not available to us, because, as explained above, only the sum of the π − π 0 , π − 3π 0 and π − π + π − π 0 spectral distributions can be combined. This sum incorporates three different branching fractions, one for each exclusive channel. The effect of the BF uncertainties is, however, very small numerically. We have checked that their inclusion changes central values of the combined three-mode contribution to ρ ud;V (s) by less than 0.5%, while the errors on this contribution are about one percent larger with than without inclusion of BF uncertainties, and the cluster values s (m) are essentially unaffected. It is thus safe to ignore potential bias issues associated with the incorporation of BF uncertainties. The same is true for the impact of the uncertainties in the external normalizing factors V ud , S EW , and B e . The errors on V ud and S EW are completely negligible, while the error on B e , at about 0.1% is small enough that it does not lead to a discernible bias.
In order to obtain an optimal combined data set, one needs to choose the clusters judiciously. Clearly, the maximum number of clusters one can choose is equal to the sum of the total number of ALEPH and OPAL data points. However, such a choice is not useful. One would find a χ 2 per degree of freedom (dof) much smaller than one, but this is not what one expects if one averages two independent experiments. Since both ALEPH and OPAL measured the spectral function over the same range in s, the number of clusters should be chosen not larger than the number of data points in each of the experiments, and, given that these two experiments are independent, 10 one expects a χ 2 /dof of order one. The goal is thus to choose a set of clusters not larger than the number of data points in either experiment that leads to a value of χ 2 (ρ)/dof ≈ 1. Of course, narrower clusters should be used in regions where the spectral function changes rapidly, such as around the ρ-meson peak.
In addition to the "global" χ 2 (ρ) function defined in Eq. (3.4), we have also looked at χ 2 (m) , the "local" χ 2 function for each cluster, since the local χ 2 (m) values may reveal discrepancies in the data sets that are hidden in the global χ 2 . The local χ 2 (m) is defined as in Eq. (3.4), but with both the data points (s i , d i ) and the data covariance matrix restricted to those data points contained in cluster m. We then evaluate all χ 2 (m) on the solution ρ (m) , m = 1, . . . , N cl , obtained by minimizing the global χ 2 (ρ). Clearly, the global χ 2 (ρ) is not equal to the sum over all clusters of the local χ 2 (m) , because the full data covariance matrix C contains entries correlating data points in different clusters. If, for a cluster k, χ 2 (k) /dof> 1, this indicates a fluctuation or a local discrepancy between the ALEPH and OPAL data.
At this stage, we will have obtained a partially-inclusive combined spectral function and associated covariance matrix for the sum of the π − π 0 , π − 3π 0 and π − π + π − π 0 modes. We still need to add the residual-mode contributions to obtain our final, updated version of ρ ud;V (s). The determination of the residual-mode contributions is detailed in the next subsection.
Note that, where the e + e − → hadrons cross sections used to infer, via CVC, the corresponding contributions to ρ ud;V (s), are given in dressed form in the original publications, these have been corrected for vacuum polarization effects to obtain the corresponding bare cross sections required as input to the CVC relation. Statistical and systematic errors on the cross sections are those reported in the relevant references. Additional information, if any, provided by the collaborations is specified below.
All results for τ exclusive-mode BFs quoted below are obtained using basis-mode BF and correlation information from the 2019 HFLAV compilation [54].
We now turn to a more detailed discussion of the determination of the residual exclusivemode contributions to ρ ud;V (s). The discussion is organized mode by mode, in the order of decreasing residual-mode BF. Readers interested only in the final result for the inclusive spectral function may skip these details and jump directly to Sec. III D below.
The π − ω(→ non-3π) contribution to ρ ud;V (s) is obtained using CVC, BaBar results [59] for the e + e − → π 0 ω cross sections, and the 2020 PDG value, 0.107(6), for the ω → non-3π BF. The BaBar cross sections are in good agreement with, and have significantly smaller errors than, those reported by SND [60]. The BaBar results produce a CVC prediction of 0.0188 (19) for the τ → π − ων τ BF, in excellent agreement with the HFLAV 2019 result, 0.01955 (65). The contribution to ρ ud;V (s) implied by the BaBar cross sections has been rescaled by the ratio of the HFLAV to the CVC BF to normalize it to the HFLAV 2019 τ BF. The BF corresponding to the resulting π − ω(→ non-3π) contribution to ρ ud;V (s) is 0.00209(14).

The KKπ contribution
Determining the KKπ contribution to ρ ud;V (s) is less straightforward since the measured distribution in τ → KKπν τ is a sum of V and A contributions, while the e + e − → KKπ cross sections are sums of I = 0 and I = 1 contributions. The V and A contributions to τ → KKπν τ cannot be separated without an angular analysis, which has not been carried out to date. BaBar [53], however, has succeeded in using a Dalitz-plot analysis to separate the I = 0 and I = 1 parts of the e + e − → KK * → KKπ cross sections, which, with the smaller e + e − → π 0 φ → π 0 KK contributions, dominate the e + e − → KKπ cross section at CM energies below m τ . ALEPH [51] has previously used the I = 1 V cross sections extracted in this analysis, together with CVC, to determine the V component of the τ → KKπν τ BF. Following the ALEPH strategy, we obtain the sum of the contributions from the three KKπ states to ρ ud;V (s) using CVC, the I = 1 e + e − → KK * → KKπ and e + e − → π 0 φ cross sections measured by BaBar [53], standard vacuum-polarization corrections to convert these to the corresponding bare cross sections, and the 2020 PDG value for the φ → KK BF. 11

The KKππ contributions
No direct experimental determination of the KKππ contributions to ρ ud;V (s) is currently available. Even were τ → KKππν τ distribution results publicly available, no obvious strategy exists for splitting this distribution into its separate V and A parts. The BF situation is, moreover, incomplete for τ → KKππν τ decays, with HFLAV listing BFs for only two of the five possible KKππ τ modes [54].
The experimental situation is more complete for e + e − → KKππ, with cross sections available for all six KKππ final states. These cross sections are, however, at-present-unknown admixtures of I = 0 and I = 1 contributions, with no known method for separating the I = 0 and I = 1 components. This precludes a CVC determination of the KKππ contribution to ρ ud;V (s). The unseparated I = 0 + 1 cross sections, and the resulting full I = 0 + 1 six-mode sum, are, however, rather accurately known. In what follows we rely on the results for this sum obtained in Ref. [13] as part of the recent dispersive determination of the hadronic vacuum polarization contribution to the anomalous magnetic moment of the muon, and provided to us by the authors [73].
In Ref. [2][3][4]51], ALEPH employed a maximally conservative approach to the KKππ contribution to ρ ud;V (s), assigning 50 ± 50% of the V + A τ → KKππν τ distribution to ρ ud;V (s). With current HFLAV 2019 values, the sum of the two currently known τ → KKππν τ BFs (those for τ → π − π 0 K 0K 0 ν τ and τ → π − π 0 K − K + ν τ ) is 0.0004154(1207). The ALEPH choice would thus correspond to a V contribution to the all-modes τ → KKππν τ BF sum, from these two modes only, of 0.000208(60)(208), where the first error is 50% of the error on the HFLAV V + A sum and the second represents the assigned 100% uncertainty on the separation of the V + A sum into V and A components.
In contrast, if one makes the analogous maximally conservative assessment and assigns 50 ± 50% of the six-mode sum of I = 0 + 1 e + e − → KKππ cross sections to I = 1, the CVC relation yields a sum of the contributions from all KKππ modes to ρ ud;V (s) which corresponds to an all-modes V τ → KKππν τ BF sum of 0.000154(5)(154), where the second error reflects the assigned, maximally conservative 100% I = 0/1 separation uncertainty. Since this constraint on the V KKππ contribution is stronger than that resulting from the alternate maximally conservative assessment based on the two measured τ → KKππν τ 12 This cannot be compared to the corresponding HFLAV 2019 version for this three-mode sum since the BF for τ → π − 5π 0 ν τ has not yet been measured. The HFLAV 2019 results for the remaining two τ → 6πν τ BFs, in addition, have non-negligible A contributions which must be subtracted in order to identify the purely V contributions. The "wrong current" A contributions to the G-parity-positive 6π states are the result of IB η → 3π decays, which cause A 2π − π + η and π − 2π 0 η states to populate the experimental 3π − 2π + π 0 and 2π − π + 3π 0 distributions. Using HFLAV 2019 results for the unsubtracted 6π BFs (excluding K 0 contributions) and the two 3πη BFs, together with 2020 PDG values for the η → π + π − π 0 and η → 3π 0 BFs, one finds for the V contributions to the τ → 3π − 2π + π 0 ν τ and τ → 2π − π + 3π 0 ν τ BFs the results 0.000113(10) and 0.000077(28), respectively. The CVC prediction for the 3-mode τ → 6πν τ V BF sum thus corresponds to a V contribution of 0.000090 (46) to the BF of τ → BFs, we employ the CVC assessment and assign as the KKππ contribution to ρ ud;V (s), 50 ± 50% of the result obtained by applying the CVC relation to the full six-mode I = 0 + 1 e + e − → KKππ cross section sum. The CVC determination has the additional advantage that it includes contributions from all KKππ modes, even those for which the corresponding τ BFs are currently unknown. 7. The (3π) − ω(→ non-3π) contributions Contributions to ρ ud;V (s) from the last of the V residual modes considered by ALEPH, (3π) − ω(→ non-3π), can be obtained from the corresponding (3π) − ω(→ 3π) contributions using the known values of the ω → 3π and ω → non-3π BFs. While the relevant τ -decay distributions have not yet been measured, BaBar [68] has determined both the e + e − → π − π + π 0 ω(→ π − π + π 0 ) contribution to the e + e − → 2π − 2π + 2π 0 cross sections and the "wrong-current" I = 0 ω(→ π − π + π 0 )η(→ π − π + π 0 ) component of that contribution [74]. We subtract this wrong-current contribution to obtain the I = 1 contributions to the π − π + π 0 ω(→ π − π + π 0 ) cross sections, use 2020 PDG versions of the ω-decay BFs to obtain the corresponding I = 1 contributions to the π − π + π 0 ω(→ non-3π) cross sections, and the CVC relation to determine the corresponding contributions to ρ ud;V (s).
3π 0 ω(→ 3π) contributions are also, in principle, present in the e + e − → π − π + 4π 0 cross sections. The preliminary SND results for the latter [70] do not include an assessment of the e + e − → 3π 0 ω(→ 3π) substate contribution. The following argument, however, shows these contributions, though not measured, must be small enough to be ignored in our CVC determination of the residual-mode (3π) − ω(→ non-3π) contribution to ρ ud;V (s).
The results of Ref. [75] show that the contribution from e + e − → π 0 ηω(→ 3π) saturates the e + e − → ηπ − π + 2π 0 cross section below s = m 2 τ . We thus take the sum of e + e − →  η2π − 2π + and e + e − → π 0 ηω cross sections as input to the CVC relation, obtaining, as a result, the sum of η(4π) − and π − ηω(→ non-3π) contributions to ρ ud;V (s), which we identify by the short-hand label ηωπη4π in what follows. The resulting contribution to ρ ud;V (s) corresponds to a very small, 0.0000017(2), result for the associated τ BF sum. This provides further support for the expectation that contributions from additional higher-multiplicity V residual modes not included in the present analysis will be entirely numerically negligible in the region below s = m 2 τ .
D. The inclusive V non-strange spectral function.
The main decision to be made when combining the data into clusters is the choice of the clusters themselves. One possibility is the basis of the strategy used in Ref. [13]. In this strategy, small groups of data points consecutive in s are assigned to clusters, after which each cluster is assigned an s value according to Eq. (3.3). The algorithm described in Sec. III B is then applied. The choice of clusters can then be varied to find the combination which has both a χ 2 /dof close to unity and small errors on the sum-rule integrals I (w) ex (s 0 ). A choice of too few data points per cluster could lead to an erratic point-to-point behavior that would not reflect any gain of information, while a choice of too many data points per cluster could lead to a loss of information. However, we should keep in mind that we are concerned with only two datasets for the 2π + 4π contribution to the spectral function that will be combined, with one data set (ALEPH) being more precise than the other (OPAL). It is then reasonable to consider a combination largely based on the ALEPH energy bins, such that the majority of clusters contains at least one ALEPH data point. The cluster sizes near the ρ peak will be narrower and widen with increasing s, since the ALEPH bin widths increase with s in the region above the ρ peak. In order to construct the full covariance matrix C to be used in the fit of Eq. (3.4) we first combine the covariances for the 2π and 4π dB/ds distributions from ALEPH and OPAL, assuming at this stage no correlations between the two data sets. Then BF errors, as well as their correlations, are added into the full covariance matrix. This introduces correlations between ALEPH and OPAL. With the full covariance matrix and a choice of clusters in hand, the fit to Eq. (3.4) can be carried out.
Our final combination contains 68 clusters (a number not too far below the 79 bins of the 2013 ALEPH data set) and yields a χ 2 per degree of freedom close to the unity: with a good p value of 15%. The result of this fit, together with the ALEPH and OPAL data sets is shown in Fig. 2, where the error bars are the non-inflated errors obtained from Eq. (3.8), and the blue band represents the inflated errors obtained through local-χ 2 inflation. 13 The local p value for each cluster is shown in Fig. 3. In the rest of this paper, we will choose to work with non-inflated errors. One expects fluctuations in p value when combining these two data sets, and the local p value is never unacceptably small. 14 There 13 By "local-χ 2 inflation" we mean rescaling the errors on those clustered data points with local χ 2 /dof greater than 1 so the modified local χ 2 /dof values become equal to 1. 14 Only one p value, equal to 0.0077, is smaller than 1%. is thus no reason to assume that values of the local χ 2 /dof signal discrepancies between the ALEPH and OPAL data. As we will see in Sec. IV, inflating the errors has no effect on α s . We explored different choices of the clusters, with the number of clusters N cl ranging from 52 to 79 (the latter equaling the number of bins in the ALEPH data set). We found no significantly better fits with N cl > 68. This is not surprising, because 68 is not much less than the total number of ALEPH data points, which form the more precise data set. In fact, we found equally good fits with N cl < 68. However, we wish to have sufficiently many values of I (w) exp (s 0 ) available to probe stability of the sum-rule fits to I (w) exp (s 0 ) (cf. Sec. IV), and thus do not want to choose N cl too small.
With the combined 2π + 4π spectral function in hand, the residual modes we need to add to obtain the inclusive spectral function, in order of decreasing BF size, are π − ω(→ non-3π), K − K 0 , ηπ − π 0 , KKπ − , 6π, KK2π, (3π) − ω(→ non-3π) and ηωπη4π. In order to add these modes, a linear interpolation to the cluster s (m) values is performed individually for each mode. In the left panel of Fig. 4 the individual contribution to the spectral function for each residual mode is shown, while the right panel shows the cumulative effect, beginning with the π − ω(→ non-3π) contribution, then adding K − K 0 , then ηπ − π 0 , and so on. From these figures we also see that the residual modes (3π) − ω(→ non-3π) and ηωπη4π with the smallest BFs already give negligible contributions to the inclusive spectral function total. This observation supports the conclusion already noted above that omitted contributions from yet-higher-multiplicity modes can be safely neglected in the region up to s = m 2 τ relevant to the current analysis.
Finally, the inclusive spectral function ρ ud;V (s) is given by the sum of the contributions from the combined 2π + 4π and interpolated residual modes. Figure 5 shows the individual contributions as well as the sum. Notice that the residual modes give only a small contribution, which, moreover, is located toward the end of the τ -decay spectrum. The final spectral function is displayed in Tab. 1 with inflated and non-inflated errors for the 2π + 4π contribution. 15 Having obtained the combined inclusive spectral function, we can compare the moments I (w) ex (s 0 ) which result to those obtained using the ALEPH and OPAL versions of ρ ud;V (s). We expect these to be consistent with one another, of course, but also expect the errors for the combined case to be the smallest. In Table 2 we compare values of the I (w) ex (s 0 ) for w = w 0 , w 2 and w 3 at s 0 = s * 0 ≈ 1.5 GeV 2 and s 0 = s * * 0 ≈ 2.9 GeV 2 . With the binning of the three data sets being somewhat different, the central values cannot be directly compared, because the spectral moments in the table have been computed at slightly different s 0 values on each line. However, the errors can be compared, as they vary only very slowly with s 0 . The table thus indicates the gain in precision for spectral moments obtained from using the combined spectral function, instead of the ALEPH or OPAL spectral functions. It should also be borne in mind that the errors quoted for the ALEPH and OPAL entries in this table do not include additional, difficult-to-quantify systematic uncertainties associated with the use by ALEPH and OPAL of MC for the s-dependences of some of the residual exclusive-mode ρ ud;V (s) contributions. This additional systematic is absent from the evaluations of the spectral moments using our updated ρ ud;V (s) since the s-dependences of the numerically relevant contributions from all but the very small residual KKππ mode (where, instead, maximally conservative experimental constraints are used) are now based on direct experimental input.

IV. THE STRONG COUPLING
In this section, we turn to the determination of α s (m τ ) from the V non-strange spectral function (and associated covariances) obtained in the previous section. We first briefly outline our strategy in Sec. IV A, then present our results in Sec. IV B. Further analysis and discussion of these results is contained in Sec. IV C.

A. Strategy
The V spectral integrals I (w) exp (s 0 ) for successive values of s 0 are highly correlated. Very strong correlations also exist between I (w) exp (s 0 ) with different weights w. We find that it is possible to carry out standard χ 2 fits taking into account all correlations when we limit ourselves to the weight w 0 , while this is not the case when fits to combinations of two different weights over the same interval in s 0 are considered.
We will thus carry out two types of fits. First, since we need sensitivity to the DV parameters in Eq. (2.12) [11], we will always include the spectral integrals with the unpinched weight w 0 in our fits. Our most basic fit is a single-weight χ 2 fit to I (w 0 ) exp (s 0 ) over an interval s 0 ∈ [s min , s max ], where we will always choose s max equal to the largest of the cluster s (m) values. Since this is a non-linear fit, errors determined from the second-derivative matrix at χ 2 min are not necessarily very meaningful, and we instead determine errors by varying each parameter such that χ 2 = χ 2 min + 1. In all cases, we find that these errors are approximately symmetric, and so take the average of the negative and positive errors as our estimate for the error on each parameter. Because correlations are fully taken into account, we will also provide the p value of these fits, and use this to determine optimal values of s min .
In the second class of fit, we carry out combined two-weight fits to I (w 0 ) exp (s 0 ) and one of the I (wn) exp (s 0 ), with n = 2, 3 or 4. These fits serve as consistency checks on the fits to I (w 0 ) exp (s 0 ) alone. For these fits, the combined two-weight correlation matrix is too singular to allow a fully correlated fit to be carried out. Thus, as in Refs. [8,11,12], we carry out "block-diagonal" fits, using a quadratic form in which all correlations between different s 0  ex (s 0 ) for w = w 0 , w = w 2 and w = w 3 at two values of s 0 , for the combined, the ALEPH, and the OPAL versions of the non-strange V spectral function. We choose s * 0 to be the closest value larger than or equal to 1.5 GeV 2 for each case (combined, ALEPH, OPAL), and s * * 0 to be the closest value smaller than or equal to 2.9 GeV 2 for each case. Note that, because the values of s * 0 and s * * 0 are slightly different for the three cases, the central values cannot be directly compared. The errors can, however, be compared, as they vary slowly with s 0 .
values for each spectral integral are retained, but not those between spectral integrals with different weights. All correlations are fully included after the fit, when we obtain parameter error estimates for such block-diagonal fits by linear error propagation, as summarized in the appendix of Ref. [11]. To distinguish it from a true χ 2 function, we will refer to the quadratic form minimized in this type of fit as the fit quality Q 2 . The distinction is relevant since Q 2 is, in general, not a χ 2 -distributed quantity. While this makes it more difficult to characterize, in a quantitative manner, the quality of such a fit, the relative "goodness" of two such fits involving the same pair of weights, but different values of s min , can still be assessed by comparing the optimized results for their Q 2 values per degree of freedom. We emphasize again that the full data covariance matrix, including now also correlations between spectral moments with different weights, is taken into account in the error propagation.
A final complication is caused by the fact that the correlation matrices for I (w 3 ) exp (s 0 ) and I (w 4 ) exp (s 0 ) turn out to have very small eigenvalues, for the relevant values of s min , precluding even the block-diagonal fits with these moments if the full set of cluster values s (m) is used for the s 0 values in the fit. It turns out that this problem can be solved by "thinning" the data, as will be described in more detail in the following subsection. We emphasize again that these two-moment fits serve primarily as consistency checks on the fully correlated fits to I (w 0 ) exp (s 0 ).

B. Results
We begin with fits to I (w 0 ) exp (s 0 ), using Eq. (2.3) with the right-hand side of this equation replaced by Eq. (2.11) with w = w 0 = 1. For the reasons explained in Sec. II C we will focus on FOPT, although we will briefly quote values obtained using CIPT as well.
The results ofour FOPT fits to I (w 0 ) exp (s 0 ) are shown in Table 3 for a range of s min values lying between 1.3246 and 1.9779 GeV 2 . The dependence of the fit results on s min is also displayed in Fig. 6 for α s (m τ ), and Fig. 7 for the DV parameters. Our first task is to determine the range of values of s min from which to obtain our best estimate of α s (m τ ). This  Table 3 as a function of s min . The yellow area correspond to the average reported in Eq. (4.1); this average is computed from the data points indicated in red (see text). The thin vertical dashed line separates the regions in which the p values shown in Table 3 are smaller than 16% (to the left), from the region where they are larger than 45% (to the right).
selection is based on several observations. First, while all p values are acceptable, those for s min ∈ [1.4251, 1.9779] GeV 2 are very good. Second, all fit parameters, including, in particular α and β, become stable for s min ≥ 1.5490 GeV 2 , while the fits become distinctly less accurate for s min ≥ 1.8249 GeV 2 , as the number of data points available to the fit becomes smaller. 16 Finally, we observe that the p value abruptly drops for s min = 1.3886 GeV 2 and becomes systematically smaller if s min is lowered below this value, signalling the expected breakdown of the theory description for low s 0 . The results of our fits are based on the spectral function without error inflation in the combined 2π +4π channels. Inflating the errors only produces even higher p values while leaving the fit parameters essentially unchanged. For example, for the fit with s min = 1.5863 GeV 2 with error inflation we find a p value of 88% and α s (m τ ) = 0.3053 (66), to be compared with 76% and α s (m τ ) = 0.3056(64) given in Table 3.
Results for the parameters obtained from fits with nearby s min are, of course, highly correlated. 17 We will take the correlated average of the parameter values at s min = 1.5490 GeV 2 , 1.6136 GeV 2 , 1.6849 GeV 2 and 1.7752 GeV 2 as our best estimate for the value of each parameter, thinning out the seven points on the interval between 1.5490 GeV 2 and 1.7752 GeV 2 to lessen the impact of the very strong correlations between the parameter values at neighboring s min . With this strategy, we find the parameter values (statistical errors only) α s (m τ ) = 0.3077 (65) , (w 0 , FOPT) (4.1) δ = 3.51 (28) , γ = 0.57(17) GeV −2 , 16 For the fit with s min = 1.9779 GeV 2 , though the result for the DV parameter γ is compatible with being positive within errors, the central value is negative, which is not allowed theoretically, as it renders the sum rule (2.11) ill defined. We take this as a sign that, at such large s min , the limited range of s 0 available is not sufficient for a reliable fit. We omit these results from Figs. 6 and 7. 17 We have calculated the correlations between the parameter values obtained from different fits with the method discussed in App. A of Ref. [ Table 3 Table 3 are smaller than 16% (to the left), from the region where they are larger than 45% (to the right). The points used in obtaining these averages are those marked in red and purple in the shaded regions of Figs. 6 and 7, respectively. We note that taking a straight average of the seven values inside the yellow window in Fig. 6 between s min = 1.5490 and 1.7752 GeV 2 , and taking the smallest parameter error on this interval as the error on this average yields α s (m τ ) = 0.3080(64), a result almost identical to that in Eq. (4.1). 18 We will take the result shown in Eq. (4.1) as our central value, with the slightly larger error shown there.
The fit to I (w 0 ) exp (s 0 ) for s min = 1.5490 GeV 2 is displayed in the left panel of Fig. 8. The right panel of the same figure shows a comparison of the representation for ρ ud;V (s) obtained using the parameters of this fit with the combined experimental result obtained in Sec. III.
We have also carried out these fits using the CIPT prescription for the D = 0 perturbative 18 If instead we take a correlated average of five α s values taking every other point starting from s min = 1.5154 GeV 2 , we find a value α s (m τ ) = 0.3060 (62). Thus, if we enlarge the window in Fig. 6 by one point on each end of the yellow "plateau," we find that our fits are stable: the central value moves by about one-fourth of the error in Eq. (4.1), and the error is essentially unchanged. The p value for our correlated averages of values inside these windows (enlarged or not) are always good.  (25)  We do not show figures equivalent to Figs. 6, 7 and 8, as they would look very similar. As is the case in all τ -based α s determinations, the CIPT value is about 5 percent larger than the FOPT value. The statistical error on the FOPT and CIPT values are about 2.1 and 2.6 percent, respectively. We note that the DV-parameter values are not significantly different between the FOPT and CIPT fits.
In Table 4, we present the results for the block-diagonal simultaneous fits to I (w 0 ) exp (s 0 ) and I  Table 3 (left panel), and the resulting spectral function (right panel), multiplied by 2π 2 . The black symbols denote data points, the red solid curve the fit, and the green dashed curve the OPE part of the fit. These parameter values are in excellent agreement with those in Eq. (4.1); of course, c 6 is new.
In the case of simultaneous block-diagonal fits to I (w 0 ) exp (s 0 ) and I (wn) exp (s 0 ) with n = 3, 4, we find that the correlation matrices for the spectral moments with the doubly pinched weights w 3,4 have very small eigenvalues, leading to unstable fits with very large Q 2 values (equal to about 16 per degree of freedom for s min ∼ 1.6 GeV 2 ). The smallest eigenvalue in each such case is around 10 −10 , orders of magnitude smaller than the smallest eigenvalue for the set of I (w 2 ) exp (s 0 ) or I (w 0 ) exp (s 0 ) integrals, which are around 10 −6 and 10 −5 , respectively.  We find that if we "thin" the set of integrals used in the fit, starting at a given s min and including only every second, third, etc., of the available higher s 0 , the Q 2 /dof drops rapidly to a value below 1, and the fit stabilizes as we increase the degree of thinning. 20 Tables 5  and 6 show the results of these fits for the cases n = 3 and n = 4, always thinning by a factor 3. Comparing the results of tables 3, 4, 5 and 6, we find good consistency among all these fits. Using the simplified averaging procedure employed above for n = 2, we find the following parameter values (statistical errors only) (70) 0.3079(70) (FOPT) , (4.4) δ = 3.43 (34) 3.41(34) , γ = 0.63 (21) 0.64 (21) [GeV where the first, respectively, second, column corresponds to a simultaneous fit to I (w 0 ) exp (s 0 ) and I (w 3 ) exp (s 0 ), respectively, I (w 0 ) exp (s 0 ) and I (w 4 ) exp (s 0 ). The results for those parameters also determined in the earlier fits also show good consistency with the values obtained in those earlier fits, reported in Eqs. (4.1) and (4.3). In addition, the values for the OPE coefficient c 6 shown in Tables 5 and 6 are consistent with those shown in Table 4. This constitutes an additional non-trivial consistency check.
We end this subsection with a comment. For reasons already explained, we did not construct the axial equivalent of the new inclusive spectral function ρ ud;V obtained in Sec. III, and thus did not carry out simultaneous fits to the V and A spectral functions. This precludes us from testing consistency between vector and axial channels, and from carrying out tests based on the Weinberg sum rules, as we did in Refs. [8,11,12]. Here we point out that such tests were always successful in the separate analyses of the ALEPH and OPAL non-strange inclusive spectral functions. We also note that our most precise results for α s were always obtained from purely V channel fits.

C. Analysis
To finalize our result for α s (m τ ), an estimate is required for the error resulting from the use of the four-or five-loop-truncated perturbation theory. This is obtained following the approach outlined at the end of Sec. II C. We focus on the single-weight fit to I (w 0 ) exp (s 0 ) with s min = 1.5490 GeV 2 .
It turns out that among the various strategies for estimating this error discussed in Sec. II C, varying c 51 by ±50% around thecentral value c 51 = 283 yields the largest, and thus most conservative, estimate of the truncation error. Symmetrizing the slightly asymmetric result produces an uncertainty of ±0.0026 on α s (m τ ). Alternate error estimates based on removing order-α 5 s terms (i.e., setting c 5m = 0), or removing both order-α 4 s and order-α 5 s terms (i.e., setting both c 4m = 0 and c 5m = 0) lead to differences equal to or smaller than the differences obtained from the 50% variation in c 51 noted above. These observations apply to the perturbative representation for I (w 0 ) th (s 0 ), and do not necessarily apply to spectral moments with other weights. Since moments with different weights have different perturbative behaviors [43,44], we will take the difference between the values of α s in Eqs. (4.1) and (4.3) to reflect an independent source of perturbative error. We multiply this difference by a factor two to take into account the fact that one of the two weights entering the fit leading to Eq. (4.3), w 0 , was also used in obtaining the results quoted in Eqs. (4.1). This leads to an additional perturbative uncertainty of ±0.0028 on α s .
Combining the statistical error of Eq. (4.1) and the two perturbative uncertainties discussed above in quadrature, we obtain our final result for α s at the τ mass scale: α s (m τ ) = 0.3077 ± 0.0065 stat ± 0.0038 pert (4.5) = 0.3077 ± 0.0075 where the subscripts "stat" and "pert" refer to the statistical and the perturbative error, respectively. As we explained in Sec. II, the τ scale is sufficiently low that non-perturbative effects are expected to be potentially non-negligible. For I (w 0 ) exp (s 0 ) non-perturbative contributions are generated by DVs, corresponding to the second term on the right-hand side of Eq. (2.11). It is interesting to quantify these effects. Even though the moment is dominated by perturbation theory, we find that the non-perturbative part of I (w 0 ) th (s 0 ), which is the moment most sensitive to non-perturbative effects, oscillates with an amplitude typically of order 20% of the α s -dependent part of the perturbative contribution (obtained by subtracting the α s -independent parton-model piece) with varying s 0 .
The non-perturbative effect is thus small but significant, and this is not surprising. The non-perturbative part accounts for the oscillation seen in the spectral function in Fig. 8 (red curve), which cannot be accounted for by the OPE (green dashed curve). We believe that it is unlikely that any variation of the DV ansatz (2.12) that does an equally good job of fitting the data would lead to a variation in α s larger than the error we obtained in Eq. (4.5). 21 It is clear that the data show the existence of non-zero DVs and, while a first-principles derivation from QCD does not exist, the main features of a DV ansatz cannot be taken to be arbitrary. As already pointed out in Sec. II, a minimal set of assumptions, based on commonly accepted properties of QCD such as, e.g., Regge behavior, leads to the parametrization (2.12) [24].
In fact, we have quantitative information on this issue, from the fits involving I (wn) exp (s 0 ) with n = 2, 3, 4, because of the single pinch in w 2 , and the double pinch in w 3,4 , which suppress DVs at different rates. Comparing the values of α s (m τ ) in Eqs. (4.3) and (4.4) to the value in Eq. (4.1), we see that the central value of α s (m τ ) varies by no more than 0.0004, i.e., 0.13% of the central value, to be compared with the 2.3% relative error in Eq. (4.5). Such variations are much smaller than we would expect were the larger DV contributions to the w 0 sum rule to have been incorrectly represented by the DV ansatz Eq. (2.12).
Running the result of Eq. (4.5) to the Z-mass scale using the standard self-consistent combination of five-loop running [33,34] with four-loop matching [78,79] at the charm and bottom thresholds (2m c (m c ) and 2m b (m b ), respectively, with MS masses from the PDG [55]) we obtain the corresponding n f = 5 result α s (m Z ) = 0.1171 ± 0.0010 (n f = 5 , FOPT) . With five-loop running and four-loop matching the uncertainty due to the running is very small. If we perform the matching at m c (m c ) and m b (m b ) we find a shift of just 0.00009, which does not contribute to the final uncertainty.
Previously, we quoted a weighted average of the two τ -based values in Eq. (4.7), of the ALEPH-based and OPAL-based results, α s (m τ ) = 0.303 ± 0.009, as our best determination from τ decays. This value and the values shown in Eq. (4.7) are in good agreement with our new, more precise value in Eq. (4.5). A direct comparison with other recent determinations of α s from τ decays [4,9] is problematic because they are all based on the truncated OPE strategy, which was shown in Refs. [25,26] to be contaminated by uncontrolled systematic effects arising mainly from the neglect of unknown higher-order terms in the OPE in Refs. [4][5][6]9]. The values of Refs. [4,9] are also highly correlated, since they are based on the same general strategy and the same ALEPH data set. We note that the values of Refs. [4,9] are significantly larger than ours α s (m Z ) = 0.1199 ± 0.0015, from Ref. [4] and α s (m Z ) = 0.1197 ± 0.0015, from Ref. [9].

V. CONCLUSION
The determination of the strong coupling from hadronic τ decays has the potential to provide one of the most precise values among the many determinations from different methods that have appeared in the literature. It thus makes sense to aim for a determination from the combined experimental information available, and this is what we set out to do in this paper. This led us to construct a new non-strange vector, isovector spectral function, which is presented in Table 1 and Fig. 5.
In order to construct this spectral function, we combined the τ → π − π 0 ν τ , τ → 2π − π + π 0 ν τ and τ → π − 3π 0 ν τ experimental data available from the ALEPH and OPAL collaborations, using the method employed before in Ref. [13]. The sum of these contributions constitutes 98% of the spectral function as measured by branching fraction.
Details of the contributions from the remaining exclusive channels, a number of which were estimated using Monte-Carlo, were not provided by ALEPH or OPAL. We have replaced the estimates for these residual-mode contributions using recent τ results for the K − K 0 mode and the large amount of data now available, via CVC, from electroproduction experiments for the remaining residual modes, with conservative estimates of the systematic errors associated with this approach. As measured by the spectral moments shown in Table 2, this leads to a more accurate determination of the spectral function ρ ud;V (s), especially in the upper part of the τ kinematic range. This is a consequence of the fact that electroproduction data are not kinematically limited near the τ mass. We emphasize that the inclusive spectral function which results is a sum of s-dependent exclusive-mode contributions, all of which are now obtained from experiment and none of which require Monte-Carlo input any more.
One of the most important applications of this new combined data set is a determination of the strong coupling α s at the τ mass scale. We employed previously developed methods using finite-energy sum rules to extract a new estimate of the MS value of α s (m τ ) from these data, which, when evolved to the Z mass scale, produces a five-flavor result with an estimated precision of about 0.8%. Our final result is α s (m Z ) = 0.1171 ± 0.0010. We also revisited the question of how to best estimate the effect of truncating the perturbative expansion for the moments involved in these sum rules, arguing that, for a direct comparison with values obtained from other methods, the fixed-order resummation scheme is the appropriate one. Because the perturbative Adler function is known to a high order, α 4 s , and the QCD β function is known to an even higher order, α 5 s , we arrived at a systematic error reflecting the use of perturbation theory which is quite small, 1.2% at the τ mass. We did not carry out an analysis of the perturbative error for CIPT, because of the issues discussed in Sec. II C.
For the same reasons, we emphasize that averaging the FOPT and CIPT values, or taking their difference as a measure of the perturbative error, would be misleading.
In order to carry out the analysis, we had to rely on the DV ansatz (2.12), and this introduces a model element into our framework. We note in this regard that, while no derivation of Eq. (2.12) from QCD exists, there are strong theoretical arguments supporting this ansatz based on commonly accepted conjectures about the spectrum of QCD [24], cf. Sec. II A. The consistency of the results of the fits presented above, which employ spectral moments with varying degrees of DV suppression, moreover, supports the consistency of our approach, as discussed in more detail in Sec. IV C.
Our approach can be subjected to further tests when more precise τ -decay data become available. With the strategy employed above, more precise data for just the low-multiplicity decay channels τ → π − π 0 ν τ , τ → 2π − π + π 0 ν τ and τ → π − 3π 0 ν τ would produce further improvements to the inclusive non-strange vector isovector spectral function and allow for much higher precision tests of our framework. High-precision data for τ → π − π 0 ν τ are, in fact, already available [81], implying that similar data for the two four-pion channels would potentially have a significant impact on the determination of α s from non-strange hadronic τ decays.
Finally, the new spectral function can also be used in other applications. It can be used to put constraints on low-energy constants in chiral perturbation theory, with nextto-next-to-leading-order theoretical representations being available [82][83][84]. Additional lowenergy constants would become accessible were an updated non-strange axial vector spectral function to be obtained. Such an update is, however, more difficult, given the absence of any analogue of the electroproduction data used in improving many of the vector residual-mode contributions. We leave the consideration of the axial channel to future work.