Muon g-2 estimates : Can one trust Effective Lagrangians and global fits ?

Previous studies have shown that the Hidden Local Symmetry (HLS) Model, supplied with appropriate symmetry breaking mechanisms, provides an Effective Lagrangian (BHLS) which encompasses a large number of processes within a unified framework; a global fit procedure allows for a simultaneous description of the e+ e- annihilation into the 6 final states $\pi^+ \pi^-$, $\pi^0 \gamma$, $\eta \gamma$, $\pi^+ \pi^- \pi^0$, $K^+ K^-$, $K_l K_s$ and includes the dipion spectrum in the {\tau} decay and some more light meson decay partial widths. The contribution to the muon anomalous magnetic moment $a_{th}$ of these annihilation channels over the range of validity of the HLS model (up to 1.05 GeV) is found much improved compared to its partner derived from integrating the measured spectra directly. However, most spectra for the process $e^+ e^- \to \pi^+ \pi^-$ undergo overall scale uncertainties which dominate the other sources, and one may suspect some bias in the dipion contribution to $a_{th}$. However, an iterated fit algorithm, shown to lead to unbiased results by a Monte Carlo study, is defined and applied succesfully to the $e^+ e^- \to \pi^+ \pi^-$ data samples from CMD2, SND, KLOE, BaBar and BESSIII. The iterated fit solution is shown to be further improved and leads to a value for $a_{\mu}$ different from $a_{exp}$ above the 4 $\sigma$ level. The contribution of the $\pi^+ \pi^-$ intermediate state up to 1.05 GeV to $a_{th}$ derived from the iterated fit benefits from an uncertainty about 3 times smaller than the corresponding usual estimate. Therefore, global fit techniques are shown to work and lead to improved unbiased results.


Introduction
As is well known, the Standard Model is the gauge theory which covers the realm of weak, electromagnetic and strong interactions among quarks, leptons and the various gauge bosons (gluons, photons, W ± , Z 0 ). In energy regions where perturbative methods apply, the Standard Model (SM) allows to yield precise estimates for several physical effects, sometimes with accuracies of the order of a few 10 −12 . In contrast, in energy regions where the non-perturbative regime of QCD is involved, getting similar precision may become challenging. This is the case for the low energy part of the photon hadronic vacuum polarization (HVP); this HVP plays a crucial role in determining the theoretical value for the muon anomalous moment a µ , one of the best measured particle properties.
Fortunately, getting precise estimates in the low energy hadron SM sector is not completely out of reach as exemplified by the Chiral Perturbation Theory (ChPT) [1,2] which is rigorously the low energy limit of QCD, valid up to 400 ÷ 500 MeV but lets the resonance region outside its scope. Lattice QCD (LQCD) is also a promising method under rapid development which already allows to perform precise computations at low (and very low) energies [3]. Interesting LQCD estimates for the HVP's of the three leptons have already been produced [4,5] which clearly show that LQCD reaches results in accord with expectations; this is especially striking for a µ with, however, still unsatisfactory uncertainties [4].
So, much progress remains to be done before LQCD evaluations can compete with the accuracy of the experimental measurements already available [6,7] or, a fortiori, with those expected in a near future at Fermilab [8,9] or, slightly later, at J-PARC [10]. Since lattice QCD is intrinsically an Euclidean approach, it is intrinsically unable to account for the existing rich amount of low energy hadronic data in the non-perturbative time-like region, i.e. from thresholds to 2 ÷ 3 GeV. Therefore, other methods, able to encompass large fractions of the physics from this important energy region, are valuable.
A natural approach to this issue is provided by Effective Lagrangians which cover the resonance region. Such Effective Lagrangians should be constructed so as to preserve the symmetry properties of QCD as already done by standard ChPT, however only valid up to the η mass region. As it includes meson resonances, the Resonance Chiral Perturbation Theory (RχPT) [11] is an appropriate framework to study e + e − annihilations from their respective thresholds up to the intermediate energy region.
It has been proven [12] that the coupling constants occuring at order p 4 in ChPT are saturated by low lying meson resonances of various kinds as soon as they can contribute. This emphasizes the role of the fundamental vector meson nonet (V) and confirms the relevance of the Vector Meson Dominance (VMD) concept in low energy physics.
On the other hand, it has been proven [11] that the Hidden Local Symmetry (HLS) Model [13] and RχPT are equivalent provided consistency with the QCD asymptotic behavior is incorporated. It thus follows that the HLS model is also a motivated and constraining QCD rooted framework. As the original HLS Model only deals with the lowest mass resonances, it provides a framework for the e + e − annihilations naturally bounded by the φ mass region -i.e. up to ≃ 1.05 GeV.
The non-anomalous [14] and anomalous [15] sectors of the HLS Model open a wide scope and can deal with a large corpus of physics processes in a unified way. However, as such, HLS cannot precisely reach the numerical precision requested by the wide ensemble of high statistics data samples collected by several sophisticated experiments on several annihilation channels. In order to achieve such a program, the HLS Lagrangian must be supplied with appropriate symmetry breaking mechanisms not present in its original formulation [13].
This was soon recognized by the HLS Model authors who first proposed the mechanism to break SU(3) symmetry [16] named BKY according to its author names. Its success was illustrated by several phenomenological studies based on the BKY breaking scheme [17,18,19]. It was also soon extended to SU(2)/Isospin symmetry breaking [20]. However, in order to account simultaneously for all the radiative decays of the light flavor mesons, the additional step of breaking the nonet symmetry for light pseudoscalar mesons was required; based on the heuristic formulation of the V P γ couplings by O'Donnell [21] which includes nonet symmetry breaking in the pseudoscalar (P ) sector in a specific way, a global and successfull account of all V P γ and P γγ couplings has been reached [22]. The BKY SU (3) breaking and this nonet symmetry breaking included within the HLS Model was shown [23] to meet the requirements of Extended Chiral Perturbation Theory [24,25]. Finally, introducing the physical vector meson fields as the eigenstates of the loop modified vector meson mass matrix provided a mixing scheme of the ρ 0 − ω − φ system which together with the V − γ loop transitions implied by the HLS model at one loop 1 leads to a satisfactory solution [27] of the long-standing τ − e + e − puzzle [28,29,30,31].
Therefore, the approach just sketched is a global framework aiming at accounting for the largest possible ensemble of data spectra collected in the largest possible number of low energy physics channels. As this global model is an Effective Lagrangian constructed from the (P and V ) fields relevant in the low energy regime of QCD and because it is consistent with the symmetries of QCD, one naturally expects their low energy results to be consistent with the Standard Model.
It was then shown, that the Effective Lagrangian constructed from the original HLS Model supplemented with the breaking schemes listed above was able to provide a satisfactory simultaneous description of the e + e − annihilations into the π + π − , π 0 γ, ηγ, π + π − π 0 final states and of the dipion spectrum in the decay of the τ lepton [32,33]. This tended to indicate that the τ − e + e − puzzle just referred was related to an incomplete incorporation of isospin symmetry breaking effects within models.
Slightly extending these breaking schemes, one is led to the Broken HLS (BHLS) Model [34] which provides a fully consistent picture of all examined e + e − annihilation cross sections 2 , the τ dipion spectrum and, additionally, some light meson decay information with a limited number of free parameters to be extracted from data. An interesting outcome of the BHLS based fit framework was a novel evaluation of the dominating low energy piece of the HVP, leading to an improved estimate of the muon anomalous magnetic moment at more than 4σ from its measured value 3 [6,7].
Introducing the dipion spectra collected in the ISR mode confirmed that the muon g − 2 departs from expectation by more than 4σ [35] . One should note that the high statistics ISR 1 See also [26] where the role of the ρ 0 − γ mixing is especially emphasized. 2 Specifically the 6 e + e − annihilation channels to π + π − , π 0 γ, ηγ, π + π − π 0 , K + K − , K 0 K 0 , each from its threshold up to 1.05 GeV, i.e. including the φ signal region. 3 One should note that the BHLS evaluation for the muon HVP is the closest to the central value preferred by the Lattice QCD study [4]. dipion spectra recently published by the KLOE [36,37,38], BaBar [39,40] and BESSIII [41] Collaborations are strongly dominated by overall scale (i.e. normalization) uncertainties; additionally the KLOE and BaBar normalization uncertainties are energy dependent. However, sizeable overall scale uncertainties raise an important issue related with their possibly biasing the physics quantity values extracted from their spectra. This issue has been identified in the reference work of G. D'Agostini [42] where a very simple case is proposed which illustrates that biasing effects can be dramatic 4 . Of course, for a key quantity like the muon g − 2, the problem should be explored and possible biases identified and fixed. The way out is already mentioned in [42] and further emphasized in other studies [45,46,44]; the exact solution exhibits a delicate issue as the removal of the bias on some quantity supposes to know its exact value. Nevertheless, as already suggested in [42] and emphasized in [44], iterative methods can be defined and are expected to be bias free; this has been applied successfully to the derivation of parton density functions in [47].
The present work mostly aims at reexamining the results provided in [34,35] concerning the muon HVP using an appropriately defined iterative fit method adapted to the dealing with form factors or cross sections in such a way that fit results and derived quantities -like the HVP, but not only -could be ascertained to be bias free. In this way, one can positively answer the question raised in the title of this study at the methodological level.
The real issue of the physics model dependence can only be answered by having at disposal results derived from several independent model frameworks, all successfully (undoubtfully) accounting for the largest possible corpus of data. Indeed, the physics correlations relating the different physics processes encompassed within a given framework cannot easily accomodate a model-independent approach. Moreover, several issues within the global fit approach are related with the formulation of isospin symmetry breaking (IB) which can hardly be made model independent, especially in a global framework.
The paper is organized as follows. In Section 2, one briefly reminds the concern of using Effective Lagrangian global frameworks in order to strengthen the constraints on the parameters to be derived from global fits. As our HLS Lagrangian framework has a range limited upward to 1.05 GeV, the brief Section 3 reminds how the full HVP is derived from fit results and from additional information.
Section 4 is, actually, the center piece of the present paper as its purpose is to define the fit method when one should deal with samples affected by strong overall scale uncertainties. This firstly turns out to precisely define the χ 2 functions to be minimized, depending on the specific properties of the spectra considered and, secondly, to set up and justify the iterative procedure we propose 5 . Subsection 4.2 puts special emphasis on the specific χ 2 function associated with samples affected by overall scale uncertainties besides a more usual experimental error matrix. The iterative fit procedure to deal with biases is formulated therein.
Most of the ISR data samples exhibit s-dependent overall scale uncertainties, which are certainly a novel feature in our field; Subsection 4.3 defines an appropriate χ 2 function suitable for such a case. Finally, Subsection 4.4 reports on the main features of the iterative global fit method when fitting sets of data samples containing samples with overall scale uncertainties of various magnitudes compared to statistical errors. The conclusions reported here rely on a Monte Carlo study outlined and illustrated in Appendix A.
Section 5 reminds the data samples used within the BHLS procedure and reports for a (minor) correction affecting the amplitudes for the annihilation channels π 0 γ and ηγ. Section 6 reports on the updated results of the fits performed using only the scan data and discarding all ISR data samples; the effects of the iterative method is illustrated here and it is shown that the needed number of iterations in the global fit procedure does not exceed 1. The more general running is the subject of Section 7 where updated results are given to correct for coding bugs affecting some of the numbers given in our [34,35]. The properties of the recently published KLOE12 [38] and BESSIII [41] data samples are examined. The evaluation of the muon g − 2 based on the iterated fits of various combinations of data samples is the subject of Section 8, where the HVP slope at s = 0 is also computed within BHLS and compared to its value directly derived from experimental data. Finally, Section 9 is devoted to conclusions and remarks.

Effective Lagrangian Frameworks And Global Fits
As reminded in the Introduction, it is a common approach to rely on the Effective Lagrangian (EL) method to cover the low energy region where QCD exhibits its non-perturbative regime and where the quark and gluon degrees of freedom are replaced by hadron fields. Each EL of practical use generally depends on parameters originating from the starting Lagrangians (like the pion decay constant f π or the universal vector coupling g) and on parameters generated by the unavoidable symmetry breaking effects (like quark mass differences); all such parameters are determined from data with various precisions.
Needless to say that any (broken) Effective Lagrangian provides amplitudes expected to account simultaneously for several different processes. This has a trivial consequence which, nevertheless, deserves to be stressed : All the Effective Lagrangians predict physics correlations among the different physical processes they can encompass : H ≡ {H i , i = 1, · · · p}. Therefore, having plugged from start the physics correlations inside the (broken) Lagrangian, the amplitudes derived herefrom should allow for a global, simultaneous and constrained fit of all available data samples covering all the channels in H. Provided the global fit is clearly successful, the parameter central values and uncertainties returned can be considered as the optimal values accounting for all the processes in H simultaneously. Therefore, one can consider that the fit information -parameter central values and error covariance matrix -exhausts the experimental information contained in the data samples covering all the processes in H.
From now on, one specializes to the Broken HLS (BHLS) model as defined and used in [34]. All data samples used in the global fit procedure defined in this paper have already been listed and analyzed in this Reference 6 ; this will not be repeated here. As for the π + π − annihilation final state, which is a central piece of HVP studies, this Reference dealt with only the available scan data which are dominated by the samples from CMD2 [52,53] and SND [54]. The samples collected in the ISR mode by Babar [55] as well as the former KLOE data samples (KLOE08 [36] and KLOE10 [37]) have been considered in [35]. Preliminary results including also the most recent KLOE sample (KLOE12) [38] have been given in [56,57]. The BESSIII spectrum [41], published by mid of 2015, is also included within our analysis.

Estimating the Muon Non-Perturbative HVP
The issue raised in this paper is whether Effective Lagrangian methods really improve the evaluation of the dominating non-perturbative part of the HVP [34,35] compared to a direct integration of experimental data (see [26,58,59] for instance). As we are working within the original HLS framework [13], what is discussed is the HVP fraction associated with the π + π − , π 0 γ, ηγ, π + π − π 0 , K + K − , K 0 K 0 intermediate states -covered by BHLS -up to ≃ 1.05 GeV; this represents more than 80% of the total LO-HVP .
Basically, the leading order (LO) non-perturbative QCD contribution to the muon HVP is estimated separately for each intermediate hadronic state H i via : and the total non-pertubative HVP component is the sum of all the possible a µ (H i ). The function K(s) in Eq. (1) is a known kernel [31] enhancing the threshold regions (s H i ) for any channel H i and σ H i (s) is the undressed cross section 7 for the e + e − → H i annhilation; s cut is an energy limit above which perturbative expansions are supposed to become valid. BHLS permits to evaluate the 6 integrals {a µ (H i ), i = 1, · · · 6} up to s φ ≃ 1.05 GeV. As the energy interval [s φ , s cut ] contribution to a µ (H i ) is beyond the BHLS energy range of validity, it is estimated using customary methods (like those defined in [58,59,60], for instance), as also the full contributions of the channels outside the present BHLS scope, like the 4 pion final states. As already stated, these pieces represent altogether about 20% of the muon LO-HVP contribution to a µ .
As can be checked by looking at the cross section formulae given in [34], most parameters to be fitted appear simultaneously in the 6 different cross sections {σ H i (s), i = 1, · · · 6} and each annihilation channel H i comes in with several experimental data samples 8 . Therefore, for instance, the data samples covering any of the π 0 γ, ηγ, π + π − π 0 , K + K − , K 0 K 0 annihilation channels play as additional constraints on the π + π − cross section and are treated on the same footing than the π + π − annihilation data themselves. On the other hand, the constraints carried by the dipion τ decay spectrum data [49,50,51] influence the fit and allow to reduce the BHLS parameter uncertainties in a consistent way 9 . This explains why the global fit method is expected to improve each a µ (H i ) contribution compared to more traditional methods -those from [26,58,59] for instance -as these ignore the inter-channel correlations revealed by the BHLS Effective Lagrangian and validated by satisfactory global fits. Of course, inter-channel 7 Final state radiation (FSR) effects also contribute and are estimated as in [31]. 8 An experimental data sample is defined as the measured spectrum m and all the uncertainties which affect it. 9 So also do the decay partial widths of the form P → γγ and V → P γ (or η ′ → ωγ) extracted from the Review of Particle Properties (RPP) [61] and implemented within BHLS.
correlations are a general feature of Effective Lagrangians, and not particular for the BHLS implementation.
As any method, the BHLS based global fit method carries specific systematics which have been examined in great details in [35]. It is worth remarking, to avoid ambiguities, that the isospin breaking effects specific of the τ dipion spectra are introduced in the dipion spectrum [35] as commonly done in the literature [62,63,64,65,66,67,68,69] (see also [26]); they are totally independent of the isospin breaking schemes involved in the BHLS Lagrangian and, actually, come supplementing these [35].

Can One Trust Global Fit Results ?
The global fit method previously used in [34,35] defines a so-called VMD strategy which can be phrased in the following way : • 1/ If the physics correlations predicted by a given Effective Lagrangian Model are supported by the experimental data they encompass, they can be considered as exact at the accuracy level reported for the data.
• 2/ Whenever the description -global fit -provided by a given Effective Lagrangian is satisfactory, the model cross sections, the fit parameter values and the parameter error covariance matrix exhaust reliably the physics information contained in the fitted data samples.
In the present case where the BHLS model is concerned, and focusing on the muon LO-HVP, Statement # 2 means that the improvements for the 6 accessibles a µ (H i ) derived from Eq. (1) by integrating from s H i to 1.05 GeV/c are legitimately valid and conceptually supported.
On the other hand, Statement # 1 does not mean that the importance of the word "Effective" is forgotten, as clear from the italic sentence it carries : Its validity might have to be revised if the experimental context evolves towards a degraded account of the data 10 .
Obviously, a VMD strategy heavily relies on the statistical methods used to analyze and fit the data; thus, one should ascertain that all aspects of the data handling are taken into account as they should. In particular, all features of the experimental uncertainties should be implemented canonically within the minimized global χ 2 and in the fitting procedure. Indeed, as remarked in [45,70], incorrect fit results are more frequently due to an incorrect dealing with the experimental errors (and correlations) rather than to the minimization procedure itself. Therefore, special care is requested in dealing with experimental uncertainties and in choosing the appropriate χ 2 expression adapted to each data sample.
It is the purpose of this Section to address this issue and check whether the procedure defined in [34,35] fulfills this statement; this will lead us to complement the fitting procedure by an iterative method. 10 However, if an ensemble of data is internally conflicting within a given Effective Lagrangian framework, as the fit results can be affected in an unpredictable way, some action has to be taken. The simplest solution is certainly to discard the faulty data samples; however, as suggested by [46], a down-weighting of the outlier contributions to the minimized χ 2 might also be considered. This could be a way to reconcile the preservation of the fit information quality with the use of all available samples.

The Basic χ 2 /Least Square Method
Usually, performing a fit -global or not -requires to minimize a χ 2 function 11 relating the differences between the measurements (m = {m i , i = 1, · · · n}) and the corresponding model (theoretical) expections (M( a) = {M i ( a), i = 1, · · · n}) weighted by the error covariance matrix V provided together with the data spectrum. Leaving aside for now possible global (additive or mutiplicative) systematic uncertainties, the error matrix V provided by experimental groups gathers the statistical and systematic errors and, thus, is not necessarily diagonal. The vector a denoting the unknown internal model parameter list, minimizing : with respect to a allows to derive its optimum value a 0 . When several independent data samples are to be treated simultaneously, the minimized χ 2 is a sum of terms like Eq. (2), one for each data sample. As reminded in [45], if the model M( a) is linear in the parameters 12 and if the error covariance matrix is correct, the estimated parameter vector a 0 has unbiased components and this estimator a 0 has the smallest variance. As illustration, in the case of a straight line fit (M = q + px), V. Blobel [45] produced the residual plots for the model parameters using several kinds of error distributions for the generated data points (each with the same standard deviation) and showed that these plots are always gaussian distributions, as expected from the Central Limit Theorem. Of course, the probability distribution is flat only if the error distributions are gaussian, i.e. if the effective χ 2 function is actually a real χ 2 .
When analyzing (a collection of) actual spectra obtained by various groups, nothing better can be done and the derived fit solution faithfully reflects the whole data information on which it relies : It corresponds, at worst, to the least square solution and, at best, to the minimum χ 2 solution, depending on the functional nature of the true experimental error distributions.

Iterative Treatment of Global Scale Uncertainties
In the Subsection just above we have briefly summarized the traditional method which applies when the handled spectra are not significantly affected by (correlated) global uncertainties. These can be of either kinds : additive (offset error) or multiplicative (scale/normalization error). As no offset error issue is reported for the spectra we analyze within BHLS [34,35], we skip this case and let the interested readers refer to suitable references [42,45,46]. In contrast, multiplicative (global scale) uncertainties are reported for most experimental spectra; when they are non-negligible compared with the other (more standard) kinds of errors, they should be specifically accounted for within the global fit procedure. This is of special concern for the important e + e − → π + π − data samples collected in scan mode [52,53,54], and even more for those collected using the Initial State Radiation (ISR) mode by KLOE [36,37,38], BaBar [39,40] or BESSIII [41]; furthermore, the normalization uncertainties reported for each of the ISR data samples have all a peculiar structure which deserves each a specific treatment -this is the subject of the next Subsection.
A constant global scale uncertainty, as those affecting the data samples from CMD2, SND or BESSIII, can be written β = 1 + λ, where λ is a random variable with range on ] − 1, +∞[. As E(λ) = 0 and E(λ 2 ) = σ 2 with σ << 1, the gaussian approximation for λ is safe [45,46]. A data sample subject to such a global scale uncertainty provides an individual contribution to an effective global χ 2 glob. which should a priori be written : where m, M, V and a carry the same definitions as in Subsection 4.1 while λ and σ have just been defined. As for A, even if intuitively one may prefer A = m, the choice A = M( a) has been shown to drop out any biasing issue 13 [42,45,70].
Assuming that the unknown scale factor λ is solely of experimental origin -and, then, independent of the model parameters a -the solution to ∂χ 2 /∂λ = 0 provides its most probable value λ 0 [34]. After substitution, Eq. (3) becomes : which exhibits a modified error covariance matrix W and only depends on the (physics) model parameters. More precisely, the single recollection of the scale uncertainty λ is the occurence of its variance σ 2 in the modified covariance matrix W . However, Eq. (4) clearly points toward a difficulty if the model is not numerically known beforehand as the modified covariance matrix becomes a-dependent when setting the unbiasing choice A = M. In this case, the parameter error covariance matrix provided by the χ 2 minimization might be uneasy to interpret.
The way out is to define iterative procedures; this is allusively stated in [42], but explicitly considered in [44] as solution to the so-called "Peelle's Pertinent Puzzle" 14 [43], provided a good starting approximate solution is known beforehand; however, defining such a tool might be a delicate task if the underlying model is non-linear, as quite usual in particle physics. Such a procedure has already been followed and successfully worked out in [47] in order to derive through a minimization procedure the parton density functions from several measured spectra. When dealing with samples of form factor and/or cross section data, other appropriate iterative methods should be defined.
The starting step of the iteration implies choosing some initial value for A, say A = A 0 . Without further information, the best approximation one can choose is obviously A 0 ≡ m, the experimental spectrum itself. Quite interestingly, this turns out to start iterating with λ = 0 (σ = 0 in Eq. (4)), i.e. β = 1, a unit scale factor; this makes the connexion with the iterative method followed in [47].
Then, the minimization of the χ 2 in Eq. (4) with A = A 0 ≡ m is performed using the MINUIT procedure [71] which yields the (step # 0) solution 15 M 0 via the fitted parameter vector value a 0 . The next step (# 1) consists in minimizing Eq. (4) using A = M 0 ≡ M( a 0 ) which 13 This does not mean that the choice A = m necessarily leads to a significantly biased solution as shown below. 14 Peelle's reference is no longer of common access, but its main content -which closely resembles the D'Agostini issue raised in [42] -is reproduced in [44]. 15 The analysis method in [34,35] actually stops there; the present analysis aims at going beyond.
is easily implemented in the procedure and, at convergence, MINUIT provides the step # 1 solution M( a 1 ). This stepwise procedure 16 . is followed until some convergence criterium is met. As in each minimization procedure the covariance matrix is constant, the interpretation of the parameter error covariance matrix is canonical. The convergence speed of the iterative procedure cannot be guessed ab initio but may be expected fast, referring to the fit of the parton density functions where the convergence is essentially reached at the first iteration [47]. This is confirmed by the Monte Carlo studies reported in Appendix A.
Nevertheless, one may infer that the number of iteration steps is smaller for a starting guess for A close to the actual model than for an arbitrary choice; clearly, as the choice A = m (the experimental spectrum) should be the closest to the actual model, one may think that it should minimize the number of interations needed to reach convergence. Additionally, this choice does not imply any a priori assumption on the parameter vector to be fitted.
Among the data samples one deals within the BHLS based global fit method, most have been collected in scan mode, essentially at Novosibirsk, and carry a constant scale uncertainty merging several effects. This is especially the case for the e + e − → π + π − data samples collected by the CMD2 [52,53] and SND [54] detectors; this also covers the case of the BESSIII data sample [41].
In order to simplify and unify the notations in the following discussion, it is suitable to perform the change of random variable λ = σµ. Then, the statistical properties for λ propagate to E(µ) = 0 and E(µ 2 ) = 1 and, defining in addition B = σA, Eq. (3) above becomes : The condition ∂χ 2 /∂µ = 0 provides the most probable value for µ : and, substituting this into Eq. (5), one gets : Stated otherwise, from the point of view of the physics model, the minimization procedure keeps track of the scale dependence by a modified covariance matrix which, in turn, influences the fit. A faithful graphical comparison of data and model -like the usual fit residual plotsshould take into account the fitted scale, as illustrated in [35] for instance.

Global Scale Uncertainties Effects in ISR Experiments
With the advent of the Φ factory in Frascati, of the J/ψ factory in Beijing and of the B factories at SLAC and KEK, the possibility opened to get large data samples for the various e + e − annihilation channels in the region of interest of the BHLS model, namely, from the thresholds to the φ meson mass energy region ( √ s ≤ 1.05 GeV). The production mechanism involved is the emission of a hard photon in the initial state [72], the so-called Initial State Radiation (ISR) phenomenon. This ISR production mode has been used to collect high statistics data samples for the e + e − → π + π − channel covering the low energies by the KLOE [36,37,38], BaBar [39,40] and BESSIII [41] Collaborations. However, it is a common feature of the KLOE and BaBar (ISR) data samples to carry non-trivial error structures. Beside a non-diagonal statistical error covariance matrix (V ), they exhibit a large number of (statistically independent) bin-to-bin correlated uncertainties, most of these being additionally s-dependent. As far as we know, this seems to be a première in particle physics and how this is dealt with inside minimization procedures deserves to be clarified and explicitely stated (see also [35]).
Let us consider a given experimental data sample E, a spectrum m function of s, for which the (given) statistical error covariance matrix is V ; the information provided for the bin-to-bin correlated uncertainties defines several independent scale uncertainties λ α (α = 1, · · · n scale ) and should be understood as follows : each of the scale uncertainty λ α is a random variable of zero mean and carrying a s-dependent standard deviation σ α (s) as tabulated by each experiment. It is clearer to make the change of (random) variables λ α = σ α (s)µ α (α = 1, · · · n scale ) and assume that all the random variables µ α fulfill E(µ α ) = 0 and E(µ α µ β ) = δ αβ .
Then, the other notations being identical to those previously defined, the χ 2 in Eq. (5) generalizes to : where implicit sum over repeated greek indices is understood. One has defined B α (s) = σ α (s)A(s), A being the s-dependent vector already defined. A is iteratively redefined as emphasized in the previous Subsection. Using the minimum χ 2 conditions ∂χ 2 /∂µ α = 0 and the independence conditions of the various sources of scale uncertainty ∂µ α /∂µ β = δ αβ , the most probable values for the µ α 's can be derived [35]. A recursion can be defined and allows to derive 17 from Eq. (8) : in close correspondence with Eq. (7). A specific feature of Eq. (9) deserves to be noted. As each experimental group reports separately on each identified independent source of (scale) uncertainty, these should indeed be fitted separately as stated just above to go from Eq. (8) to Eq. (9). More precisely, for the experiment E, we are not using the quadratic sum (σ E (s)) 2 = α [σ α (s)] 2 for its partial χ 2 , which would have given σ E (s i )σ E (s j )A i A j inside the full error covariance matrix instead of what is shown in Eq. (9). Stated otherwise, the various sources of normalization uncertainties are not summed in quadrature but really treated as statistically independent.

Numerical Tests of the Global Fit Iterative Method
As stated in the header of the present Section, if the physics correlations predicted by the Effective Lagrangian (here BHLS) are fulfilled by the data, the estimate of the model parame- 17 For clarity, defining Z = A or B, Z k denotes the quantity Z(s k ) for short. ters and the parameter error covariance matrix are legitimate tools to serve evaluating related physics quantities.
As in the previous studies relying on the HLS model, at the early stages [32,33] or more recently [34,35,56,57], the method is to minimize a global χ 2 expression taking into account the largest possible number of data samples and using appropriately all information provided by the experimentalists concerning all kinds of uncertainties which affect their data samples. The aim of Subsections 4.1 -4.3 was to detail how the χ 2 piece associated with each data sample should be constructed, depending on its reported error structure.
In contrast with previous references (including ours), the fit procedure will be adapted in the present study in order to examine and cure possible biases produced by having stopped the fit procedure at the A = m step instead of iterating further on as suggested in [42], explicitly proposed in [44] and performed in [47].
In order to check whether estimates based on global fit results can be trusted as, for instance, the muon HVP central value and its uncertainty derived from the fit information returned by MINUIT, some additional checks on the fitting method and its iterative aspect deserve to be performed, at least to control that, indeed : • The parameter pulls are centered gaussians of unit standard deviations.
One should also check that the fit probabilities distributions are uniformly distributed on [0, 1] when the measurements are indeed true unbiased gaussian distributions. This condition list can be supplemented with some examination of the effects due to nonlinear dependences upon the parameters to be fitted.
However, checking this list of properties obviously implies that the true parameter values are known, that the measurements are indeed sampled on truely centered gaussian distributions and that their errors are indeed the true standard deviations of the measured spectrum. Stated otherwise, this exercise goes beyond using actual measured experimental data samples as, then, truth is unknown : The global fit method -as any other method -should be evaluated using data samples generated by Monte Carlo techniques; in this case, the true parameter values and their uncertainties are known at the sample generation level and can reliably be compared to the fit results. The detailed study is transferred to Appendix A; the most involved results are summarized right now : • The effects of non-linear parameter dependence within models used to fit data spectra (see Subsection A.2.1) are likely to be marginal for the kind of experimental distributions we are dealing with. This should be related with the local minimum finding structure of the algorithms gathered within the MINUIT package.
• When scale uncertainties dominate the sets of spectra globally submitted to fit, using 18 A = m E gives a solution which can exhibit strong biases, but this solution is the start of an iterative procedure which leads rapidly to the unbiased solution to the minimization problem. The biases occuring at start of the procedure can be very large, but they are observed to practically vanish already at the first iteration step (the previously called M 1 solution).
• When performing a global fit of some data samples dominated by global scale uncertainties together with others where the statistical errors (e.g. affecting randomly each bin) dominate, the iterative method obviously works as well as just stated. In this case, however, the presence of some samples free from scale errors exhibits an unexpected pattern : Even if the data samples free from scale uncertainties are affected by enlarged statistical errors, they strongly reduce the biases generated by the A = m E choice. Stated otherwise, the effects of data samples where the normalization errors are dominated by the (random) statistical errors is to favor the smearing out of the biases in the parameter value estimations.
The properties just listed concerning the unbiasing of the fit parameters extend to the estimates of physics quantitites derived from using the fit result information (parameter values and error covariance matrix). Additionally, as the parameter pulls are observed as centered gaussians of unit standard deviation, the calculated uncertainties relying on Monte Carlo sampling of the fit parameter distributions should also be reliable. This is of special relevance for the evaluation of the various contributions to the muon LO-HVP discussed in Section 3.
The last item in the list just above has important consequences while working with real (and so, not really perfect) experimental data. However, even if the fraction of data samples free from -or marginally affected by -scale uncertainties may look large enough, it is nevertheless cautious to ascertain that the fit solution is indeed unbiased by performing one or two additional iterations. Indeed, the studies reported in Appendix A tell that, anyway, the iterated fit solutions are always unbiased.
Therefore, one may conclude from this Section and from the simulation studies reported in Appendix A that global fit methods can indeed be trusted. The single proviso is that iterating the fit procedure as explained above is mandatory or, at least, cautious.
The issue is now to examine how the results given in [34] and [35] are modified when iterating beyond the approximation A E = m E for all data samples significantly affected by scale uncertainties, constant (as, mostly, the spectra reported in [52,53,54] ) or s-dependent (as all the ISR spectra reported in [36,37,38,39]). Observing the stabilizing effect of the data samples dominated by statistical errors (like the γπ 0 and γη final states) is also methodologically relevant.

BHLS Global Fit Method : Present Status and Corrigendum
As stated several times above, the Effective Lagrangian Model we use is the broken HLS (BHLS) model developped in [34]. In this Reference, the BHLS model is also applied to all data samples collected in scan mode, by the various Collaborations which have run on the successive Novosibirsk e + e − colliders. These e + e − annihilation samples cover the π + π − , π 0 γ, ηγ, π + π − π 0 , K + K − , K 0 K 0 final states and have been discussed in detail in several previous studies [32,33,34]; for the sake of conciseness, we will not repeat this exercise here. As the BHLS model also covers the τ decays from the early stages of its formulation [27], the previous studies include the dipion spectra collected in the τ ± → π ± π 0 ν τ decay mode by ALEPH [49,73], Belle [51] and CLEO [50]. Also included within the BHLS fit procedure are some light meson decay partial widths not connected with the annihilation channels already listed, like A second step has been to extend the study in [34] to treat the high statistics ISR data samples for e + e − → π + π − ; this has been the purpose of the study in [35] where the KLOE08 [36] and KLOE10 [37] data samples collected by the KLOE Collaboration and the data sample produced by BaBar [39] have been examined. Since then, two new samples have been produced by the KLOE (KLOE12 [38]) and BESSIII [41] Collaborations 19 Except otherwise stated, all the fit results presented in this paper have been obtained using the Configuration B [34] (i.e. dropping out from the fit procedure the three pion data samples collected in the φ mass region).
The studies covered by [34,35,56,57] rely on minimizing a global χ 2 function summing up partial χ 2 's, each associated with a given data sample. For each of the ≃ 40 ÷ 50 data samples, the partial χ 2 was (canonically) constructed following the rules detailed in Section 4. However, as the fit was not iterated in the studies [34,35], it is worth checking to which extent the value of the muon HVP derived herefrom is changed by the iteration procedure.
For the present study, a few coding bug fixes have been performed and a piece missing in the expression for the e + e − → π 0 γ and e + e − → ηγ cross sections has been included. So, when different, the results in the present paper supersede those in [34,35].
As for the missing piece just mentioned : In the amplitudes γ * → γP 0 (Eq. (65) in [34]) and the cross section formulae e + e − → γP 0 (Eq. (68) in [34]), the non-resonant piece should be modified as follows : This implies that the single process which depends separately on the FKTUY [15] parameters c 3 and c 4 is the e + e − → π + π − π 0 annihilation. In this case both c 3 +c 4 and c 3 −c 4 combinations come in, while all others quantities only involve the c 3 + c 4 combination 20 . We apologize for the inconvenience.

BHLS Global Fit Method : Iterating with NSK Data Only
In this Section, we report on global fits using the data reminded in the preceding Section and discussed in [34]; as for the pion form factor data, we focus for the present exercise on using only the most recent scan data collected by CMD2 and SND [74,52,53,54], excluding the older data samples from OLYA and CMD [75].
partial χ 2 's are essentially given by expressions like Eq. (4). For the other data samples, we performed as in [34]. The first data column in Table 1 displays the results of the fit performed by setting A = m in the χ 2 associated with each experimental data spectrum generically named m. The form factor returned by this (A = m) global fit is named M 0 and is used to perform the first iterated (A = M 0 ) global fit; the results of this fit are shown in the data column #2; this iteration #1 global fit returns the solution named M 1 . The iterated #2 fit is then performed by setting A = M 1 in the χ 2 expressions of the pion form factor data samples, leading to another (M 2 ) solution; the fit results are displayed in the third data column in Table 1.
One clearly observes a quite tiny change in the first iteration : 0.2 unit in the χ 2 value of the π + π − data samples; also the global χ 2 changes by only 0.7 unit. When going from the first to the second iteration, the changes are almost invisible. This corresponds for experimental data to the effect reported in Subsection A.2.3 for our Monte Carlo data. As derived quantity, let us report on the leading order (LO) contribution a µ (ππ) derived by integrating Eq. (1) between 0.63 GeV/c and 0.958 GeV/c; using obvious notations, the previously reported fits yield : in units of 10 −10 . So, one observes a tiny effect while iterating once (0.3% for the central value) and no effect when iterating twice. In the present case, where the former data from [75] have been dropped out from the fit, the "experimental" estimate is a µ (ππ, [0.63, 0.958]) = 361.26 ± 2.66 (see Table 7 in [34]). Another way to account for the scale uncertainty is to set A = M( a) (which depends on the parameters under fit) and perform the fit. A starting value for A must be chosen (denoted A start ) but its value changes at each step of the minimization procedure. In this case, the fit convergence time is much larger than previously but the results are almost identical to those already obtained by iterating. The last 2 columns in Table 1 display the fit results starting with A start = M 1 and also those starting from the fit solution derived herefrom (denoted M x ). As for a µ (ππ, [0.63, 0.958]), the values derived in these last fits numerically coincide with the iterated cases displayed above.
Therefore, one may indeed conclude, as can be inferred from the Monte Carlo studies reported in Appendix A, that the HVP value reached without iterating is very close to the HVP derived from the once iterated solution. One also observes, as expected, that iterating only once already leads to the final result; indeed, from iteration #1 to iteration #2, the changes for a µ (ππ) are at the level of a few 10 −12 .
As for the fit quality reflected by the χ 2 values at minimum and the corresponding fit probabilities, the last line in Table 1 indicates that, whatever is the way one treats the vector A, they are all alike. This, once more, corresponds to expectations, as can be checked with the discussion in Subsection A.2.3 and expecially the properties of Figure 8. Nevertheless, it is useful to check that the twice iterated solution does not modify the result derived from the once iterated solution in a significant way.

BHLS Global Fit Method : Iterating Scan and ISR Data
It remains to introduce the other π + π − data samples collected at e + e − colliders using the ISR mechanism. Reference [35] has already done this work with the data samples then available using the method described in Subsection 4.3 without, however, iterating the procedure. The conclusion reached was that the KLOE08 [36] and BaBar [39] data samples have difficulties to accomodate -within the BHLS framework -the whole set of data samples covering the channels already reminded in Section 5. In contrast, the KLOE10 [37] data sample was found to fit well the BHLS expectations. Complementing preliminary works [56,57], we revisit here the issue with the two new data samples provided by KLOE (KLOE12) and BESSIII.

The τ +PDG Analysis
In Ref. [35], it has been shown that the BHLS fitter can be run without explicitly using definite e + e − → π + π − data samples besides the non π + π − channels. Indeed, on general grounds, one expects that some limited isospin breaking (IB) information specific of this annihilation channel can make the job together with the τ dipion spectra. It has been shown that the partial widths Γ(ω/φ → π + π − ) and Γ(ρ 0 → e + e − ), together with the products (V = ω, φ) Γ(V → π + π − ) × Γ(V → e + e − ) represent an amount of information sufficient to reconstruct -within BHLS -the pion form factor in the e + e − channel. Figure 1: The τ +PDG prediction (red curve) of the pion form factor in e + e − annihilations in the ρ − ω interference region. The various superimposed data samples are not fitted; also displayed are the average χ 2 distances of each of the e + e − → π + π − data samples to the common τ +PDG prediction.
Before going on, it deserves noting that the decay information used to run the τ +PDG method has been extracted from the Review of Particle Properties (RPP) [61] and that the above mentioned pieces of information are in no way influenced by the data collected by KLOE, BaBar or BESSIII; actually, they are almost 100 % determined by the data samples from the CMD-2 and SND experiments. On the other hand, the τ +PDG analysis is not influenced by the global scale issue which mostly motivates the present work.
We have performed the τ +PDG run using all annihilation data mentioned in the above Sections (configuration A [34]). The fit returns χ 2 τ /N τ = 82.1/85 = 0.97. The best fit solution allows to reconstruct the predicted invariant mass distribution of the pion form factor in the e + e − → π + π − annihilation; this prediction is expected valid over the whole BHLS range as shown by Figure 2 in [35]. It is worth showing here the mass range from 0.70 to 0.85 GeV; Figure 1 displays the τ +PDG prediction on this range together with the available π + π − data superimposed (and not fitted); we have calculated the χ 2 distance of each sample over its full range 21 . The average χ 2 per data point is indicated inside the corresponding pannel. Figure 1 indicates that the average χ 2 distances for the NSK (CMD-2 & SND), KLOE10, KLOE12 and BESSIII samples are small enough to claim a success of the τ +PDG method. One can conclude that they fulfill the consistency issue discussed in Section 2 with the full set of data and channels covered by BHLS. One should note that the description of the BESSIII sample (which is not a fit) is as good as the fit published by the BESSIII Collaboration [41]. For KLOE08 and BaBar, we reach the same conclusion as in [35]; nevertheless, one can now compare the behavior the twin 22 samples KLOE08 and KLOE12 : We have χ 2 KLOE08 = 4.8 while χ 2 KLOE12 = 1.2 clearly reflecting a better understanding of the error covariance matrix, while the central values are almost unchanged, as clear from Figure 1.
Stated otherwise, the issue met with KLOE08 and BaBar is confirmed but the two new data samples published since [35] are both found is good correspondence with expectations.

The Iterative Method : Global Fit Properties
The issue is now to report on the behavior of the global fits performed using the iterated method when the π + π − ISR and scan data are considered simultaneously; this complements the work already presented in Section 6 when using the scan data only. Except otherwise stated, the τ data samples are always included into the fit procedure. On the other hand, as the behavior of the global fit for data/channels other than π + π − does not differ sensitively from the information already displayed in Table 1, this will not be repeated. Table 2 displays our main results using the scan and ISR e + e − → π + π − annihilation data. They correspond to the iteration # 1 fit (denoted above A = M 0 ), however the previously called A = m or A = M 1 solutions gives almost identical fit quality results 23 .
The first data line displays the global fit properties with the indicated e + e − → π + π − data samples used each in isolation within the global BHLS context, together with all other data samples covering the rest of the encompassed physics (see Section 5).
One observes that the average (partial) χ 2 per data point χ 2 π + π − /N π + π − is of the order 1 or (much) better and the probability high when running with any of the KLOE10, KLOE12,

Fit Configuration
Iteration Method    Table 2: Global fit results as function of the e + e − → π + π − data sample content. Each entry displays the [χ 2 π + π − /N π + π − ] value returned by the global fit. The data samples involved can be tracked from the column titles, the following line giving the corresponding data point numbers [N π + π − ] in the range up to 1 GeV. The value flagged by * has been obtained using a BaBar sample truncated from the energy region [0.76, 0.80] GeV (250 data points).
NSK 24 and BESSIII data samples; as in [35] the picture is not as good for KLOE08 and BaBar.
Performing a global BHLS fit using the data samples from KLOE10, KLOE12, BESSIII, NSK and BaBar (amputated 25 from the energy region [0.76, 0.80] GeV) leads to results given at the entry lines flagged by "Fit Combination 1"; as the correlations between the KLOE08 and KLOE12 data samples are strong and their content not explicitely stated 26 , it is more cautious to avoid dealing with the KLOE08 and KLOE12 samples simultaneously. Despite the removal of the drop-off region in the BaBar π + π − spectrum, the global fit quality looks poorer.
The results obtained when using the KLOE10, KLOE12, NSK samples within the fit procedure are displayed at the Entry "Fit Combination 2" when BESSIII data are also included and "Fit Combination 3" when they are not; the data and fit corresponding to the "Fit Combination 2" are shown in Figure 2. Both Fit Combination 2 and Fit Combination 3 are clearly satisfactory. 24 NSK here denotes the collection of data samples from CMD2 [74,52,53], SND [54] (127 data points in total) as well as the former (82 data points) samples collected by OLYA and CMD [75]. The numbers in Table 2 given within square brackets include the contributions from these former samples. 25 We remind that this removal is motivated by a possible mismatch in the energy calibration in the ρ 0 − ω interference region between BaBar and the other π + π − data samples submitted to the same global framework. In contrast, when running with the π + π − BaBar sample in isolation, its full spectrum is considered. 26 Some work in this field seems ongoing [76]. The e + e − → π + π − data samples are those shown in the entry "Fit Combination 2" in Table 2.
The inset in the top panel magnifies the ρ 0 − ω peak region. The dowmost panels magnify the behavior in both distribution wings. See Section 7.2 for further comments.
Therefore, this proves that the scan data from CMD2 and SND are consistent with the KLOE10, KLOE12 and BESSIII data samples and that all these are fully consistent with the other data spectra introduced in the global fit procedure as indicated by the global fit probability. One should also remark that the systematic uncertainties provided for KLOE12 lead to a satisfactory global fit, in contrast with KLOE08, as already noted in the previous Subsection.
Except otherwise stated, the fit parameter values presented from now on are derived using the e + e − → π + π − data samples corresponding to the "Fit Combination 2" (see Table 2); the fit results are those derived after the first iteration and they do not differ significantly from the corresponding results at iteration # 2. The fit quality for the non-π + π − data samples are almost undistinguishable from the numbers already given in the second data column from Table 1; they are not repeated for the sake of brevity.

The Iterative Method : Updating The Model Parameter Values
Beside improving the fits by mean of the iterative method, the present work accounts for an error and a couple of bugs affecting our [34,35]. Moreover, the present work includes the new KLOE12 data sample within the fit procedure; this is not harmless as KLOE12 constrains the fit conditions more severely than the KLOE10 sample. Therefore, the present results update and supersede the corresponding ones previously given in [34,35].

The HLS-FKTUY Parameters
The non-anomalous HLS Lagrangian (broken or not) can be written : The unbroken expression for L HLS can be found in [13] and its broken expression (BHLS) is given in [34]. The covariant derivative which allows to construct both pieces of L HLS introduces the fundamental parameter g, known as universal vector coupling. The coefficient a HLS is a specific feature of the HLS model, expected close to 2 in standard VMD approaches; however, phenomenology rather favors a HLS ≃ 2.5 since the early applications of the HLS model to pion form factor studies [77,78,23].
On the other hand, the anomalous (FKTUY) sector [15] of the HLS model [13] consists of 5 pieces (see also Appendix D in [34]), each weighted by a specific numerical parameter not fixed by the theory. Using common notations [13,34] and factoring out, for convenience, the weighting factors, the FKTUY Lagrangian collecting all the anomalous couplings can be written 27 : 27 Actually, the erratum involved in Eq. (10) comes from having missed the contribution of the (c 4 − c 3 ) term displayed in Eq. (13) which actually turned out to impose c 4 = c 3 . As already stated, after correction, all the anomalous decay couplings and the amplitudes for e + e − → (π 0 /η)γ anihilations only depend on the combination (c 4 + c 3 )/2 and the single place where the difference (c 4 − c 3 ) occurs is the e + e − → π 0 π + π − annihilation amplitude. In [34,35] where (c 4 − c 3 ) was absent, its physical effect was absorbed by (c 1 − c 2 ) to recover good fit qualities; so (c 4 − c 3 ) and (c 1 − c 2 ) should carry an important correlation.
where P and V indicate the basic pseudoscalar and vector meson nonets and A the electromagnetic field. As L HLS , L F KT U Y depends on the universal vector coupling g.
At iteration # 1, the global BHLS fit returns : with correlation coefficients never larger than the percent level, except for < δg δa HLS >= −0. Our value for c + agrees with the estimates derived in [13] from the π 0 γγ * ) form factor (c + = 1.06 ± 0.13) and from the ω → π 0 γ partial width (c + = 0.99 ± 0.16) with a much smaller uncertainty due to the large amount of data influencing the (global) fit. After the bug fixing, c − is found small but non-zero with a large significance and (c 1 − c 2 ) becomes closer to 1. Using the full 25 × 25 parameter error covariance matrix returned by the global fit, we have computed separately c 4 and c 3 by a Monte-Carlo sampling. This gives c 3 = 1.124 ± 0.022 and c 4 = 0.789 ± 0.021.
Among the numbers displayed in Eq. (14), some are appealing : The nearness to 1 of the fitted c 1 − c 2 and c + parameters, their customary guessed value [13], should be noted and deserves confirmation with more precise data on the anomalous annihilations and light meson radiative decays than those presently available.

The Iterative Method : Pseudoscalar Meson Mixing and Decay Parameters
The BHLS symmetry breaking of the Lagrangian piece L A leads to pseudoscalar physical fields constructed as linear combinations of their bare partners. The mechanism involved is the BKY mechanism extended so as to account for both Isospin and SU(3) symmetry breakings [34]; it can be complemented by the pseudoscalar nonet symmetry breaking scheme generated by the t'Hooft determinant terms [79]. The main effect of these determinant terms is to provide the bare Lagrangian with a correction to the PS singlet kinetic energy term governed by a parameter λ expected small (see Eq. (7) in [34]).
The BHLS model connects to (Extended) ChPT [25,24], especially its two angle θ 0 and θ 8 mixing scheme; in particular, it relates these angles to the singlet-octet mixing angle traditionally denoted θ P , together with the BKY breaking parameters z A , ∆ A and to λ [34].
The upper part of Table 3 displays in its first data column our fit results in the general case. The fit value for θ 8 is in good agreement with other expectations [24] as well as that for θ 0 . The smallness of this has led us to impose θ 0 = 0 within fits which leads to the results shown in  Table 3: Some parameter values derived when leaving free θ P and λ (first data column) or when relating them by imposing θ 0 = 0 to the fit (second data column).
the second data column. The value for λ undergoes a severe correction compared with [34,35] and, presently, because of its large uncertainty, could be neglected without any real degradation in fit qualities. BHLS also allows for some additional contribution to the π 0 − η − η ′ mixing based on some possible aspects of Isospin breaking not already accounted for by the extended BKY scheme developped in [34]. This turns out to redefine the physical (observable) fields (right-hand side) in terms of the (BHLS) renormalized (left-hand side) fields by [80] : Inspired by [80], one can lessen the number of free parameters by stating : and fit ǫ 0 . Then, using the fit results (parameter central values and error covariance matrix), one can reconstruct the value for ε and ε ′ . The updated values are given in Table 3 still indicate a π 0 − η mixing much larger than the π 0 − η ′ mixing (a factor of 4).
Before closing this Subsection, we mention that the Monte Carlo sampling method allows to reconstruct the decay constant ratio f K /f π = 1.265 ± 0.009 which becomes f K /f π = 1.295 ± 0.002 when constraining the fit with θ 0 = 0.

The Muon LO-HVP : Evaluations From Iterated Fits
The main aim of the present study is to produce improved estimates of the muon LO-HVP [ 34,35] by means of the iterated global fit method expected to cancel out possible biasing effects which could affect the A = m (i.e. non-iterated) solution. The validity of the iterated method is supported by the Monte Carlo study outlined in Appendix A, which clearly indicates that the iterated method cancels out possible biases and returns, correctly estimated, the fit parameter uncertainties. Therefore, building on the conclusions collected in Subsection 4.4 one can produce bias free evaluations of the muon LO-HVP. The effects of iterating 28 from M 0 to M 1 -the solution derived using A = M 0 within the fit procedure -will be especially emphasized. To be complete, this update also takes into account the new KLOE12 [38] and BESSIII [41] π + π − data samples -which happen to be very constraining -and also corrects for some bugs. Therefore the present numerical results supersede the corresponding ones in [34,35]. a µ (ππ, [0.63, 0.

958] GeV)
The point at top of Figure 3 is the so-called τ +PDG [35] value for a µ (ππ, [0.63, 0.958] GeV) derived by switching off the contributions of the various e + e − → π + π − data samples from the minimized χ 2 , replacing them by decay information extracted from the Review of Particle Properties (RPP) [61] as emphasized in Subsection 7.1.
In order to get the other points displayed in Figure 3, one always uses all the channels covered by BHLS, including the τ spectra from ALEPH, CLEO and Belle. As for the e + e − → π + π − data samples, one uses each of the BaBar, KLOE08, KLOE10 and KLOE12 samples in isolation as indicated within the Figure (see also Table 2 and Subsection 7.2). The point flagged by CMD2+SND is obtained from a fit to the so-called [34] new timelike data from CMD2 and SND [74,52,53,54], leaving aside the older data from OLYA and CMD collected in [75] (see Table 2 and Section 6 above). As for the BaBar spectrum, for reasons already stated, the fit is performed on the spectrum amputated from the drop off region ( √ s ∈ [0.76, 0.80] GeV). Finally, as the published BESSIII spectrum ends up at 0.9 GeV, one cannot produce an experimental value on the interval [0.63, 0.958] GeV. As a general statement, Figure 3 clearly illustrates that the iterated (M 1 ) and the noniterated (M 0 ) solutions provide quite similar fit estimates for a µ (ππ, [0.63, 0.958] GeV). One should nevertheless remark that the agreement between both fit solutions and the numerical integral of the experimental data is less satisfactory for the data samples which exhibit poor fit qualities within the global framework (KLOE08 and BaBar) than for the others (KLOE10, KLOE12, CMD2+SND) as can be inferred from the "fit in isolation" properties displayed in Table 2. Finally, the weighted averages of the experimental results for KLOE10 and KLOE12 alone or together with all NSK data (the so-called new timelike data and the former samples [75]) are always well reproduced by the global fit and are supported by quite good probabilities (see Table 2).
Using the NSK+KLOE(10/12) sample configuration, the iterated BHLS global fit gives a slightly smaller central value (by ≃ 1.5 10 −10 ) while the uncertainty is improved by a factor ≃ 2. It is also worth pointing out the role of the τ spectra within the BHLS global fit framework. The following numbers illustrate how the constraints involved by the τ spectra allow BHLS to yield a more precise fit estimate for a µ (ππ, [0.63, 0.958] GeV). Comparing the direct integration result to the values derived from fits, one indeed gets at iteration # 1 : Finally, the downmost point in Figure 3 displays the result derived using all data samples (except for KLOE08 as there is not enough published information to account for its strong correlation with KLOE12); this estimate for a µ (ππ, [0.63, 0.958]) which benefits from a very small uncertainty has, however, a poor fit probability as clear from Table 2.  Table 4: The contributions to the muon LO-HVP from the various channels covered by BHLS from their respective thresholds to 1.05 GeV in units of 10 −10 at start and after iteration. The last column displays the direct numerical integration of the various spectra used within BHLS. The π + π − data samples considered are those flagged by "Combination 2" in Table 2.

Contributions To The Muon LO-HVP Up To 1.05 GeV
The LO-HVP's integrated from their respective thresholds up to 1.05 GeV are displayed in Table 4; the central value for a µ (ππ) includes final state radiation (FSR) effects. The first data column shows the results from the fit solution M 0 derived from fitting with A = m; the second data column displays the results corresponding to the solution M 1 derived by fitting with A = M 0 . These two data columns report on the fits performed using all annhilation channels encompassed by BHLS and the τ dipion spectra. Finally, the rightmost data column provides the direct numerical integration of the experimental spectra -actually only those feeding the BHLS fit procedure, including the KLOE10, KLOE12 and BESSIII data samples besides the scan data.
As for the π + π − channel, both fits -which include the τ spectra -provide central values in agreement with each other and with the direct estimate within the quoted error 29 . If the A = m solution were (inherently) exhibiting a bias, comparing the first two numbers in the first line of Table 4 indicates that this does not exceed ≃ 0.5×10 −10 -e.g. half a standard deviation. Therefore, real experimental data samples confirm the gain provided by a global fit procedure when samples with normalization errors small compared to their statistical accuracies are included; exploring this effect is the purpose of Subsection A.2.3 in Appendix A.
One should also remark that the unbiasing iterative procedure lessens significantly the uncertainty on a µ (π + π − ) compared with the A = m solution and, over the whole range of validity of BHLS (up to 1.05 GeV), one ends up with a factor of ≃ 3 reduction of the uncertainty compared to the direct numerical integration. The same kind of effect is reported in [47] concerning the spread of the parton density functions 30 .
Therefore, relying on the iterative procedure, one observes that the global fit does not produce significant shifts of the central values of the HVP contributions which could be attributed to the normalization (scale) uncertainties strongly affecting some data samples. Relying on the Monte Carlo studies outlined in Appendix A, this can be attributed to the large number of data samples where the statistical uncertainties dominate over the normalization uncertainty. Moreover, the uncertainty on the part of the LO-HVP derived from the BHLS fit (more than 80% of the total LO-HVP) is very small and even marginal.

The Muon g − 2 From BHLS Global Fit Procedure
In order to evaluate the muon LO-HVP from the fit results derived by means of the BHLS global fit procedure, the numbers given in Table 4 should be supplied with several additional contributions which cannot be derived from within the BHLS framework but should be estimated by other means. This covers the channels opened below 1.05 GeV but remaining outside the present BHLS scope 31 and, more importantly, all hadronic contributions covering the non-perturbative QCD region above 1.05 GeV should be estimated via the direct integration method. Table 5 summarizes these additional contributions to be combined with the BHLS results to derive the muon LO-HVP; in this Table, one reminds the information available by end of 2011 and used in our previous [34,35]. The data column flagged by "LO-HVP (2014)" is the update derived by taking into account the data samples more recently collected (and published up to the end of 2014); these are the e + e − → 3(π + π − ) data from CMD-3 [81], the e + e − → ωπ 0 → π 0 π 0 γ from SND [82] and several data samples collected by BaBar in the ISR mode 32 [83,84,85,86]. These data samples highly increase the available statistics for the annihilation channels opened above 1.05 GeV and lead to significant improvements. One thus should note the important improvement these provide for the LO-HVP contribution from the [1.05, 2.0] GeV region : its uncertainty is reduced by 25 %, while its central value is almost unchanged. Despite this improvement, the energy region [1.05, 2.0] GeV still remains Indeed, roughly speaking, the experimental cross section σ exp (s) is related with the underlying theoretical cross section σ th (s) by a relation of the form σ exp (s) = σ th (s) + δσ(s) and the δσ(s) correction depends on the normalization uncertainties which just motivate the iterative method! Actually, this δσ(s) is exactly the scale dependent term in Eqs. (5) and (8). Obviously it cannot be estimated without some fitting procedure. 30 In particular, Figure 5 in this Reference, is quite informative about the variety of correction kinds revealed by unbiasing procedures. 31 For instance the 4, 5 of 6 pion annihilation channels, or the ωπ 0 final state. 32 These cover the pp, K + K − , K L K S , K L K S π + π − , K S K S π + π − , K S K S K + K − annihilation final states.   Deriving the full HVP value also requires to account for the higher order effets. This includes the next-to-leading order contribution (NLO) taken from [26] ([−9.97±0.09]×10 −10 ) and the recently estimated next-to-next-to-leading order (NNLO) effects which happen to be non-negligible ([1.24 ± 0.01] × 10 −10 ) [87].
To compute the muon g − 2, one should also include the light-by-light (LBL) contribution (here taken from [88]), the QED contribution [89,90] and the electroweak contribution (EW) [31]. The next-to-leading order contribution to the LBL amplitude (NLO-LBL) has also been computed recently [91] but is clearly negligible ([0.3 ± 0.2] × 10 −10 ). Altogether, the numerical values we use (see Table 6) are rather consensual [92]. Significance (nσ) 4.65σ 5.00σ 4.38σ Table 6: The various contributions to 10 10 a µ . ∆a µ = a exp µ − a th µ is given in units of 10 −10 . For the measured value a exp µ , we have adopted the value reported in the RPP which uses the updated value for λ = µ µ /µ p recommended by the CODATA group [93]. By KLOE, one means that the KLOE10 and KLOE12 π + π − data samples are introduced in the BHLS fit procedure and in the directly integrated spectra.
The first data column in Table 6 reproduces (after our methodological update) the muon anomalous moment estimate coming from the corresponding BHLS global fit where only the scan data for the π + π − channel are considered while all ISR data are excluded. This supersedes the corresponding information in [34]. The sample combination preferred by the BHLS global fit gives the results displayed in the second data column; it exhibits a 4.9σ significance for a non-zero ∆a µ = a exp µ − a th µ . The evaluation derived by direct integration of the spectra used within the global fits are given in the third data column. The new data, as a whole, increase the discrepancy for ∆a µ which is always found above the 4σ level; effects of additional and not still accounted for systematics will be examined in the next Subsection. Figure 4 displays the results for ∆a µ derived using or not the τ data and various combinations of the available π + π − data samples introduced within the BHLS global fit procedure at first iteration. For comparison, one also displays in this Figure the evaluations produced by other authors and flagged by Dhea09 [29], DHMZ10 [58], JS11 [26] and HLMNT11 [60]corrected however for the recently calculated NNLO-HVP and NLO-LBL -contributions as included in Table 6. A priori, the Dhea09 estimate compares exactly to our evaluations using scan data only; the other results are derived using, beside the NSK samples, the BaBar, KLOE08 and KLOE10 samples. These may compare to the last couple of lines in Figure  4 where the scan data are supplemented with the BaBar (not truncated), KLOE (10/12) and BESSIII samples.
The following comments are in order here : • 1/ The difference between our estimates and those of other authors mainly concerns the estimated central value for ∆a µ . Also, our uncertainties are now reduced because of the global fit method, but also because of using much more data samples than other authors; this is clear by comparing the errors shown in Figure 4 with those given in [35].
When using only the scan data, the shift one observes should reflect the biasing effect certainly present in the experimental data (see footnote # 29) and corrected in our approach by the iterated fit method. When the ISR π + π − samples are also involved, the issue just reminded is amplified because the weight of samples with large overall scale uncertainties is much increased 33 . The effect of the BaBar data sample is no longer enough to balance the effect of the new data samples as clear by comparing the lines for "NSK+KLOE+BESSIII" with the lines for "Global (ISR+scan)" which also include the (full) BaBar sample. Nevertheless, one should note the large difference of the correponding probabilities.
• 2/ When a comparison between a ∆a µ estimate derived using the τ data and the corresponding one excluding these is possible, ours exhibits the smallest difference (1.12 × 10 −10 for NSK+KLOE+BESSIII, −0.7 × 10 −10 for the Global fit including all the π + π − data samples). This is certainly due to the vector meson mixing which defines the BHLS model. It is interesting to note that the JS11 [26] value, which is based on the γ − ρ 0 mixing by loop transitions 34 , is the closest to ours.
• 3/ Relying on the global fit properties, the BHLS model favors the "NSK + KLOE10 + KLOE12 +BESSIII + τ " as the largest consistent set of data samples. This leads to ∆a µ = (37.55 ± 4.12) × 10 −10 which exhibits a 5.σ significance 35 . Our estimate is 33 All ISR data samples are strongly dominated by overall scale uncertainties, additionally s-dependent. 34 Within the BHLS model too, the γ − ρ 0 mixing is mediated by loop effects. 35 If using the data from 2011 in Table 5, as in our previous studies, this significance is "only" 4.8σ. This compares more directly to the results from other authors displayed in Figure (4). The increased significance is a pure consequence of the recent improvements of the hadronic contribution from the [1.05, 2.0] GeV region. expected to be free of biases generated by the overall scale uncertainties which dominate the ISR π + π − data samples.  The various a th µ have been derived from the global fit using the indicated e + e − → π + π − data samples and including/excluding the τ dipion spectra as indicated. In red we display ∆a µ corresponding to the iterated solution and in green those corresponding to the A = m (non-iterated) solution. In blue results from other studies are given corrected by the recently evaluated next-to-next-to-leading order contribution [87]. See Section 8.3 for comments.

Additional Systematics On The BHLS Estimate For The Muon g − 2
A detailed study of additional systematics possibly affecting the BHLS evaluation of ∆a µ has been already performed in [35]. It concluded to an uncertainty of the LO-HVP central value for ∆a µ = a exp µ − a th µ in the range [−1.3 ÷ 0.60] × 10 −10 coming from π + π − contribution in the φ mass region, where BHLS is weakly constrained. An uncertainty coming from using the τ spectra has also been considered; it was argued that the best motivated evaluation of this is the difference between fitting with the τ spectra and without them in the most constrained configuration. Presently, this means that the BHLS preferred value (∆a µ = (38.58 ± 5.04) × 10 −10 ) could be underestimated by ≃ 0.9 × 10 −10 .
Another mean to detect systematics is to compare with the accurate ChPT predictions on the P -wave π + π − phase-shift [94] and also with the available experimental data from the Cern-Munich [95] and Fermilab [96] groups. These are shown in Figure 5. Included also are the predictions derived from the Roy Equations [97] and from the phase of the pion form factor fit performed in [26] (JS11).
As for the BHLS predictions corresponding to using NSK+KLOE(10/12), we display in this Figure the phase of the full amplitude and those corresponding to dropping out the isospin breaking (IB) effects due to the vector meson mixing 36 . The τ spectra are included within the fit procedure. Figure 5: P -wave π + π − phase-shift data and predictions from [94] (CGL) and [26] (JS11) together with the BHLS phase-shift. The insets magnify the various behaviors close to threshold. See Subsection 8.4 for further explanations.
The standard BHLS phase shift predictions are displayed in the left-hand side panel of Figure 5. One clearly observes a very good prediction of the phase-shift up to about 1.2 GeV, i.e. much beyond our fitting range (from threshold to 1.0 GeV for the ππ data). Indeed the Cern-Munich data are very well accounted for and the BHLS predictions are in accord with the other predictions. The inset, however, exhibits a (minor) issue for the full amplitude phase, a small bump of about 1 • close to threshold, absent from the IB amputated amplitude. This can be tracked back to a peculiarity of the broken HLS model which does not split up the HK (Lagrangian) masses for the ω and ρ 0 mesons and, consequently, the mixing angle α(s) does not exactly vanish at s = 0 (see Figure 6 in [32]); in contrast the other angles fulfill β(0) = γ(0) = 0. Indeed, one has : where [34] ǫ 1 (s) is the difference of the charged and neutral kaon loops and Π ππ (s) is the pion loop which both vanish at s = 0. This assumption has been checked with fits by imposing [m HK ω ] 2 = (1 + η)[m HK ρ ] 2 and choosing various fixed values for η; the right-hand side panel in Figure 5 displays the phase shift for η = 5% and, quite satisfactorily, its inset does not reveal a bump any longer. A non-zero (HK) mass difference η [m HK ρ ] 2 cannot be generated by the breaking mechanisms already implemented within BHLS. However, a breaking of the nonet symmetry in the vector meson sector (VNSB) enables such an effect; this turns out to modify the customary vector field matrix -actually U(3) symmetric -within the covariant derivatives of the HLS model [13] by a perturbation term proportional to the singlet vector field combination. The effect of VNSB has been derived from specific fit studies and indicates that ∆a µ might have to be lessened by about 1.4 × 10 −10 .
Therefore, in total, the BHLS favored result can be expressed, in units of 10 −10 as : where the three additional contributions play as shifts on the central value. Adding them up linearly, the maximum shift (−2.7 × 10 −10 ) may reduce the central value to 34.85 × 10 −10 which has still a 4.6σ significance. The effect of these additional systematics is to reduce potentially by ≃ 0.3σ all the significances displayed in Figure 4. These are not due to overall scale uncertainties already accounted for by the iterative method; they might be reduced by new annihilation data samples covering the region up to 1.05 GeV in all the physics channels in the realm of BHLS.

The HVP Slope At Origin In BHLS Fits
In the lattice QCD approach of calculating a had µ , extrapolation methods have been developed (see e.g. contributions to [98]) to overcome difficulties to reach the physical point in the space of extrapolations. The low Q 2 behavior of the euclidian electromagnetic current correlators on a lattice, which exhibits a discrete momentum spectrum, poses a particular challenge (see e.g. [99,100] and References just below). The analysis of moments of the subtracted photon vacuum polarization function Π(Q 2 ) was particularly advocated in variants in Refs. [101] and [102]. Recent lattice calculations [103,104,105,106] have been utilizing moment analysis techniques for a more precise evaluation of a had µ . The leading moment is given by the slope of the Adler function [107,108], the latter being given by : where R(s) is the hadronic spectral function 37 and s min the smallest threshold energy squared (s min = m 2 π 0 within BHLS). Then, defining : the HVP slope at the origin is given by : The constant P 1 can be directly estimated from data and partly from the BHLS fits. Therefore, one can proceed as done above with our evaluations of a had µ and derive the results gathered 38 in Table 7. Here, one observes that the difference between the experimental and the HLS values for the HVP slope are at the percent level (a 2 σ exp effect) and the uncertainty is scaled down by a factor of 10. However, to really feel the HLS improvement on the slope, one needs once more an improved hadronic spectral function at high energies.  A lattice estimate of the Adler function slope D ′ (0) has been presented in [109]. The result is P 1 = 5.8(5) GeV −2 , and has been compared with P 1 = 9.81(30) GeV −2 , a result estimated using a phenomenological toy-model representation [110] of the isovector spectral function. The lattice results too include the isovector part only and is missing higher energy contributions above 1 GeV.
In the study [102], the authors provide numerical values from fits to Lattice data based on Padé approximants (PA). For this purpose, they parametrize the HVP as : which thus leads to : The parameters corresponding to the results they consider as optimal are given in their Table 3.

Concluding Remarks
The present study was motivated by the question which gives its title to this paper. More precisely, the issue is whether the D'Agostini bias [42,46] prevents to derive unbiased physical results from global fits to experimental spectra affected by dominant overall scale uncertainties 40 .
Actually, several issues are merged together. First, the effective global χ 2 functions to be used in the minimization procedure should be appropriately defined. For the data samples where the statistical errors dominate the overall scale uncertainties, the construction of the associated partial χ 2 's is quite standard. The real issue starts when the data samples are dominated by overall scale uncertainties. For each of them, substantially, the canonical partial χ 2 has been reminded in Section 2 and writes [42,45,46] : leaving aside the so-called "penalty term" [46] proportional to λ 2 . The (partial) χ 2 being appropriately defined, another issue is the choice of the vector A.
In our former studies [34,35], beside the ≃ 40 data samples dominated by statistical errors which follow the traditional treatment, the data samples covering the e + e − → π + π − annihilation channel are all, sometime very strongly, dominated by overall scale uncertainties; this especially refers to the samples collected by the KLOE and BaBar Collaborations using the ISR production mode. Here, for each sample, we chose for A the experimental spectrum itself; this choice has been referred to as A = m all along the paper. The guess behind was that all scale uncertainties affecting the different experimental spectra independently of each other should smear out possible biases in the central values of the (common) theoretical form factor function parameters [35].
It happens that the results one can derive in this way from the BHLS global fit undergo very small biases (compared to the errors derived from the fit procedure); this is shown in the present study 41 . However, the guess just reminded was incorrect and the actual reason which explains the almost bias free results is following : As shown in the Monte Carlo study presented in the Appendix, there is no smearing out of biases if all the spectra submitted to fit undergo comparable strong scale uncertainties; however, this study also shows that, if some of the fitted spectra are dominated by (random) statistical errors rather than global scale uncertainties, the fit results can be strongly unbiased.
Nevertheless, a high level of unbiasing cannot be taken as granted as the real weight of the samples dominated by statistical errors within the full global fit procedure cannot be ascertained beforehand. Basically, the choice A = m potentially leads to biases of unknown magnitude; this has been shown by G. D'Agostini [42] with a simple example and more generally argued by V. Blobel [46]. These authors also showed that all biases vanish if, instead of A = m, one makes the choice A = M, the "true" spectrum. But this is just not possible within contexts like ours, where fits are performed just in order to derive the "true" spectrum from data. Fortunately, iterative methods allow to circumvent this difficulty by taking the path opened in [47] in order to derive the parton density function from data and correct for biases. The iterative method we propose has been tested with the Monte Carlo study reported in the Appendix and shown to produce unbiased results with a quite fast convergence speed; indeed, only one iteration is sufficient.
So, our main conclusion is indeed that global fit methods including a fast iterative procedure are expected to produce reliable pieces of information as, methodologically, the central values are unbiased and the estimate for the uncertainties reliable; this especially applies to the part of the muon leading order HVP derived from e + e − annihilation cross sections.
Having shown that appropriate global fit methods should lead to results which can be trusted, a related remark is worth being expressed. Iterative global fits allow to supply the BHLS Effective Lagrangian cross sections with reliable and unbiased numerical central values for the fit parameters and a good estimate of their error covariance matrix. Then, using these cross sections and the fit information, Eq. (1) is expected to provide an unbiased estimate for a µ (ππ) as the ingredients are unbiased.
On the other hand, when computing a µ (ππ) by directly integrating a dipion spectrum in order to derive its so-called experimental value, one has to plug into Eq. (1) the experimentally measured cross section σ exp. (s). However, as already noted in footnote # 29, or as can be inferred from the canonical χ 2 expression reminded just above, the experimental and model cross sections are related by : where the best estimate of the second term writes 42 δσ(s) = λσ theor. (s). As obvious from Eq. (6), the best estimate of the scale factor λ equally depends on the measured spectrum and on the "true" spectrum, which can be identified with its (iterated) fit solution. So, using again self-explanatory notations, Eq. (1) leads to : a µ (ππ, exp.) = a µ (ππ, theor.) + δa µ (ππ) and thus a µ (ππ, exp.) looks intrinsically biased for any sample subject to strong enough overall scale uncertainties. This issue is also reflected by the residual plots which are improved when plotting the corrected residuals [m − (1 + λ)M( a)] instead of the raw ones [m − M( a)], as can be seen in Figure 13 of [35]; this allows to infer that δa µ (ππ) is small but non-zero. It amounts to δa µ (ππ) ≃ 2 × 10 −10 in the case "NSK+KLOE10+KLOE12+BESSIII+τ " favored by the BHLS model.
As for the physics conclusions, the present paper updates and corrects the results derived by the global BHLS fit method which, following the considerings just summarized, has been completed with an iteration procedure in order to cancel out possible biases. One thus confirms that almost all of the existing data samples covering the annihilation channels with the π 0 γ, ηγ, π + π − π 0 , K + K − , K 0 K 0 final states and the dipion spectra in the τ ± → π ± π 0 ν decay accomodate perfectly the BHLS framework. In the line of our previous works, one also finds that among the data samples covering the e + e − → π + π − annihilation, the data samples provided by CMD2 and SND, the KLOE10 and now also the KLOE12 and BESSIII samples behave consistently with each other and with the other considered data covering the various channels entering the BHLS scope.
The present update, which also includes the recently published KLOE12 and BESSIII π + π − samples, supersedes our previous results; these are mostly given in Table 3 and in Eqs. (14). From a theoretical point of view, it is interesting to note the corrected values for the c i 's coefficients of the anomalous (FKTUY) terms of the HLS model [15,13] : The combinations c + = (c 4 + c 3 )/2 and c 1 − c 2 are found very close to the usually assumed value, i.e. 1; in contrast, c − = (c 4 − c 3 )/2 = −0.166 ± 0.021 is non-zero with a 8σ significance. Figure 3 displays the values for a µ (ππ, [0.63, 0.958]) GeV derived from iterating the fits with the various available data samples. One observes a strong reduction of the uncertainty compared to the corresponding experimental value (about a factor of 2.5) and there is a close agreement between central values for all samples (or combinations of samples) which yield a good fit probability. The difference between the central values for the starting fit and the iterated one tends to indicate that biases are limited; this should be a consequence of also dealing with a large number of samples where the overall scale uncertainties are dominated by random statistical errors, as argued in the Appendix. Figure 4 exhibits the values for the muon ∆a µ = a exp µ − a th µ when various combinations of e + e − → π + π − and τ ± → π ± π 0 ν samples are used in the iterated global fit procedure. The present study confirms that, within BHLS and because of its specific isospin breaking mechanisms, one does not observe any serious mismatch between fits with only e + e − annihilation data and fits where these are supplemented with the τ dipion spectra. The central values 43 for a µ (e + e − ) and a µ (e + e − + τ ) only differ by 2 units (NKS), 1 unit (NSK+KLOE+BESSIII+τ ) or 0.7 unit in the global fit of all data samples (including BaBar) as can be seen in Figure 4. Figure 4 displays the value for ∆a µ derived using all data samples except for KLOE08, which can be written : where an estimate of the magnitude of possible uncertainties coming from outside the BHLS framework is proposed. This exhibits a 5σ significance (which may reduce to 4.6σ -in the least favorable case -if the additional systematics are added linearly and assumed to play as a shift). One should note however that the fit probability is poor. The most probable value for the muon ∆a µ is obtained by using the CMD2, SND, KLOE10, KLOE12 and BESSIII samples -and the τ spectra; this leads to : This BHLS preferred estimate exhibits a 5.σ significance for a non-zero ∆a µ , which may reduce to 4.7σ if one takes into account, as just above, the possible additional systematics. This solution is associated with a 99% fit probability. As a summary, even complemented with an iterative procedure shown in the Appendix to remove biases, the BHLS approach favors a significance for ∆a µ above the ≃ 4.5σ level; this value is a lower bound obtained by including possible additional systematics added linearly. New data expected soon may further clarify the picture. The uncertainties now become sharply dominated by the region above 1.05 GeV, i.e. outside the BHLS scope.

A.1 The Test Method
In order to test the iterative method, one has developped a minimization code which deals with spectra generated from a given underlying function M true (s) where the parameters {a i } (which, of course, are known at the generation level) are fitted within the code. The "experimental" spectra feeding this code are generated using the true distribution smeared by introducing gaussian uncertainty distributions. Indeed, for the purpose of testing our analysis method, it is certainly the most appropriate to rely on "perfect" data samples, with perfectly known properties.
For sake of simplicity, at the generation level, any "experimental" spectrum E is chosen to carry 100 "measurements" m E i , performed at 100 equally spaced energy squared s i points (s i ∈ [0, 1] GeV 2 ), the same sequence for all spectra. The "measurements" are derived by smearing the theoretical values M true (s i ) in the following way : For each spectrum E, one assumes the "measurements" are sampled out from gaussian distributions in the following way : where ε i,E stat (0, 1) indicates the i th sampling on a gaussian distribution of 0 mean and unit standard deviation generating the statistical error; it varies independently from "measurement" to "measurement" and from spectrum to spectrum. ηM true (s i ) denotes the statistical error common to all m i , η being some fixed fraction of the order of a few percents, chosen the same for all the "measurements" in the spectrum E.
On the other hand, λ E = σε E scale (0, 1) is the scale uncertainty affecting specifically the spectrum E; as indicated by its definition, it is sampled out from a gaussian distribution of zero mean and σ standard deviation. The overall scale uncertainty affecting E is obtained via one sampling of ε E scale (0, 1) which, thus, carries the same value for all the "measurements" m E i in the spectrum E. Of course, when going from a spectrum E to another E ′ , another sampling of ε E scale (0, 1) should be performed. For specific tests, the overall scale uncertainty can be switched off (σ = 0).
One defines N rep replicas (generally 1000) of N exp (generally 5) experimental spectra constructed as shown in Eq. (25) and submitted to a global fit where the parameters entering M true (s) are just the parameters to be derived from the fit. The "true" statistical error covariance matrix V ij = [ηM true (s i )] 2 δ ij is practically approximated by V ij = [ηm E i ] 2 δ ij ; we have avoided the unessential complication of non-diagonal covariance matrix. The fit results derived for each replica are stored and then used to construct the statistical plots -true residuals and pulls -with the help of the known parameter "true" values.
Therefore, we are just in the conditions described in Subsection 4.2. One should note that the MINUIT code we have built performs the minimization of the N exp samples and runs sequentially to treat the N rep replicas within the same job.
So, for each replica, the global χ 2 minimized by our Monte Carlo MINUIT procedure is simply a sum of N exp terms like Eq. (4): When initializing the iteration procedure, one uses A E = m E , i.e. the spectrum E serves to construct its χ 2 E ; so A E differs from some other A E ′ by statistical fluctuations. When iterating, at first or higher order, they become identical as Obviously, each such run provides simultaneously all the information allowing to examine the statistical properties of the iterative method corresponding to a given theoretical choice M true (s). The computer code also allows an easy change of the functional form of M true (s) in order to examine the behavior of various kinds of non-linear parameter dependences.
The behavior of the fit parameters compared to truth is, of course, the subject of the analysis; however, those of "physics quantities" derived from them are as important. For this purpose, we chose to examine the ratios 44 : which has properties similar to those of the a µ (H i )'s, as the weighting factor K(s) in Eq. (1) is an unessential complication while looking for possible methodological biases of the iterative method.

A.2 The Test Results
The aim of the present Appendix is to report on numerical analyses performed in various configurations in order to examine how overall (global) normalization uncertainties and biases are related and whether non-linearities in the model parameters to be fitted lead to significant incorrect estimates of errors. As Reference [47] which is faced with the same kinds of issues as the present work, we do not plan to establish rigorously general theorems on these topics -assuming the scope of the issues would permit it. Nevertheless, one can think that studying methods by relying on Monte Carlo technics is an acceptable way to check its (practical) validity under common conditions. After all, the fact that Eq. (3) with A = M (the theoretical function) is considered free from biases is not weakened by the fact that the general (formal) proof of this property -if established -is not commonly referred to.

A.2.1 Analytical Shape of the True Distributions
In order to use confidently fit results derived using the iterative method, one should examine the effects of non-linear dependences upon the fit parameters within contexts similar to our physics distributions. The lineshape of the pion form factor as a function of s on a given interval can be qualitatively reproduced using polynomials, ratios of polynomials, exponential of polynomials, sums of a Breit-Wigner function with polynomials etc . . . with appropriate numerical parameter values.
We have applied the method outlined in Subsection A.1 to perform fits relying on an intensive use of the tools provided by MINUIT taking various kinds of functions M true (s), resembling -sometimes weakly -the pion form factor. Running in sequence MIGRAD/HESSE and MINOS, we did not observe significant departures (beyond statistical fluctuations) from equality between parabolic and MINOS errors; as the issue was to examine effects of non-linear parameter dependences this exercise was performed assuming statistical uncertainties only. Therefore, 44 Remind that 0 and 1 GeV 2 are the energy squared limits of the generated spectra. this led us to conclude that, for the kind of experimental distributions one deals with, non-linear effects are not generally significant. For instance, using : η = 3% and no scale uncertainty (to discard any need for iterating), the probability distribution was observed flat and the parameter pulls consistent with normal gaussians G(m = 0, σ = 1); the distribution of the ratio I for the 1000 replicas was also found well centered at 1 (actually its mean is 1.0001 and its standard deviation 1.62 × 10 −3 from a gaussian fit with χ 2 /N points = 8.9/11). So, except for pathological cases which may always occur, non-linear dependences do not look practically an issue. From now on, we limit ourselves to reporting on using M true (s) as given by Eq. (28). Moreover, for sake of succinctness, we may only mention the fit parameter residual and pull distribution properties qualitatively and concentrate on discussing the distribution of the ratios I which, in fine carries -summarized -the relevant information. Each value of I entering this distribution is computed from a MINUIT fit of N exp = 5 data samples and this is done for N rep = 1000 replicas to construct numerically its distribution.

A.2.2 Normalization Uncertainty and Iterative Method
We first examined the results derived by fit of spectra with data points generated as in Eq. (25) with a statistical uncertainty η = 3% and generating the scale uncertainty λ with σ = 5%; so η is smaller than σ. In this case, the interesting plots are gathered in Figure 6.
As one knows M true (s), one can construct the N exp partial χ 2 's with A = M true (s) (see Eq.(3)) and minimize their sum using MINUIT. In this case, no bias is expected [42,45,46] and this is indeed confirmed by the top left panel in Figure 6 where the distribution of the N rep values for I is displayed.
When, instead, one uses A = m (the data spectrum), the results are shown in the top right panel of Figure 6, where one observes a shift of the central value by as large as 20% ! Denoting the result of the corresponding fit by M 0 , one restarts fitting the same data by setting A = M 0 , this -first -iteration leads to the distribution shown in the bottom left panel of Figure 6 which looks identical to having used A = M true . Denoting the fit solution of this first iteration by M 1 , one restarts fitting the same data by setting A = M 1 , and get the step 2 solution M 2 which correponds to the bottom right panel of Figure 6, which clearly indicates no change for the I distribution.
So, one may conclude that the iterative procedure has already converged at the first iteration and so, we have M 1 = M true . This fortunate high convergence speed has also been observed by [47] and it is quite remarkable that this has allowed to recover from 45 a 20% bias! Fit residuals are observed unbiased and pulls consistent with normal centered gaussians for A = M true , A = M 0 and A = M 1 . As for the χ 2 probability distributions, for A = m, it exhibits a huge spike at 1, while it is consistent with flatness (mean ≃ 0.5 and r.m.s. ≃ 1/ √ 12) for all the other cases. This already indicates that starting with A = m (the measured data spectrum) and iterating only once allows to give up knowing the theoretical function M beforehand to drop out biases in physics quantity estimates. Moreover, as the parameter pulls are centered gaussians of unit standard deviations, the uncertainties derived from from the fit parameter error covariance matrix are reliable.

A.2.3 Effects of Subsamples Free from Normalization Uncertainties
In the specific problem of globally fitting a large number of experimental data samples, one is faced with as many as 40 to 50 spectra to be treated [34,35,57,56]. Within this ensemble of data samples, one observes several configurations concerning uncertainties : some samples have statistical errors dominated by scale uncertainties (the ISR collected data samples), while, in contrast, some others are reported with scale uncertainties marginal compared to statistical errors (the e + e − → γP data, for instance); sometimes, no specific information is reported concerning scale uncertainties, as for the the τ dipion spectra [49,50,51]. Figure 7: Effect of having 1 (top panels) or 2 (downmost panels) data sample(s) among the fitted N exp = 5 samples simultaneously fitted. Left plots report on fitting with A = M (the truth), right plots on fitting with A = m E (the measured spectra); in the former case no bias is observed, in the latter case, the bias happens to be much limited. See text for more details.
This makes interesting to examine configurations mixing samples of both kinds. In this paragraph, one summarizes the results obtained by running N rep replicas of ensembles of 4 data sets where, as before, the scale error is σ = 5% and the statistical error η = 3%, together with 1 data set with σ = 0% (no scale uncertainty) and η = 6% (twice worse statistical precision). This (4,1) combination will be supplemented with a (3,2) combination with the same characteristics. The main results are shown in Figure 7. Here we do not report on iterating the fit procedure, as obviously the results will follow the pattern shown in Figure 6. . The top panel plots correspond to the case when among the N exp = 5 fitted spectra, one is systematically free from normalization uncertainty; in the downmost panels 2 of the 5 fitted spectra are free from normalization uncertainty.
The top panels in Figure 7 display distributions of the ratios I in the (4,1) configuration. The left plot shows the case when the N rep replicas are fitted using A = M truth in the χ 2 expressions. In this case, the absence of any bias is confirmed by the gaussian fit result shown within this plot. While using A = m E , the top right panel exhibits a 1.3% bias. Therefore, the effect of a single spectrum free from scale uncertainty out of 5 is enough to lessen dramatically the observed bias : It reduces from 20% to 1.3%.
The downmost panels in Figure 7 display the corresponding results when fitting N rep replicas of (3,2) combinations. In this case, using A = m E in the minimized χ 2 expression, leads to an even smaller bias (0.5%).
So, even if they carry a poor statistical precision, having some spectra free from a (significant) scale uncertainty is quite helpfull to strongly limit the real magnitude of a possible bias for a derived quantity. It is a quite interesting property to observe that some spectra with degraded statistical quality supplementing other spectra dominated by scale uncertainties might be enough to avoid the need of an iteration procedure to unbias physics pieces of information.
As for the probability distributions, comparing of the corresponding left and right panels in Figure 8 clearly shows that the departures from uniformity (i.e. average=0.5 and r.m.s.=0.289) due to using A = m E are quite limited.
Nevertheless, when dealing with true experimental data (and thus unknown truth), one cannot take as granted that the number of samples with negligible scale uncertainties compared to statistical errors is sufficient to ascertain that biases are negligible. Therefore, in the practical case of the global fit of real experimental data performed within BHLS, secure results can only be ascertained by iterating until the change of a µ (H i ) is small enough.