Upgraded breaking of the HLS model: a full solution to the τ−e+e− and ϕ decay issues and its consequences on g−2 VMD estimates

The muon anomalous magnetic moment aμ and the hadronic vacuum polarization are examined using data analyzed within the framework of a suitably broken HLS model. The analysis relies on all available scan data samples and leaves provisionally aside the existing ISR data. Our HLS model based global fit approach allows for a better check of consistency between data sets and we investigate how results depend on different strategies which may be followed. Relying on global fit qualities, we find several acceptable solutions leading to ambiguities in the reconstructed value for (aμ)th. Among these, the most conservative solution is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$a_{\mu}^{\mathrm{had,LO}}[\mathrm{HLS\ improved}]=687.72(4.63) \times10^{-10}$\end{document} and (aμ)th=11 659 175.37(5.31)×10−10 corresponding to a 4.1σ significance for the difference Δaμ=(aμ)exp−(aμ)th. It is also shown that the various contributions accessible through the model yield uniformly a factor 2 improvement of their uncertainty. The breaking procedure implemented in the HLS model is an extension of the former procedure based on a mechanism defined by Bando, Kugo and Yamawaki (BKY). This yields a quite satisfactory simultaneous description of most e+e− annihilation channels up to and including the ϕ meson (π+π−, π0γ, ηγ, π+π−π0, K+K−, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$K^{0} \overline{K}^{0}$\end{document}) and of a set of 10 (mostly radiative) decay widths of light mesons. It also allows to achieve the proof of consistency between the e+e−→π+π− annihilation and the τ±→π±π0ν decay and gives a solution to the reported problem concerning the measured partial width ratio \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\varGamma(\phi\to K^{+}K^{-})/\varGamma(\phi\to K^{0} \overline{K}^{0})$\end{document}. Prospects for improving the VMD based estimates of aμ are emphasized.


Introduction
The muon anomalous magnetic moment a μ is a physics piece of information which has been measured with the remarkable accuracy of 6.3 × 10 −10 [1,2]. From a theoretical point of view, a μ is the sum of several contributions; the most prominent contributions can be predicted with a very high accuracy by the Standard Model. This covers the QED contribution which presently reaches an accuracy better than 1.6 × 10 −12 [3] or the electroweak contribution where the precision is now 1.8 × 10 −11 [4]. The light-by-light contribution to a μ is more complicated to estimate and is currently known with an accuracy of 2.6 × 10 −10 [5].
Another important contribution to a μ is the hadronic vacuum polarization (HVP). Perturbative QCD allows to compute a part of this with an accuracy of the order 10 −11 ; this covers the high energy tail and the perturbative window between the J /ψ and Υ resonance regions. For the region below this threshold, one is in the non-perturbative region of QCD where estimates of the hadronic VP cannot so far be directly derived from QCD, relying on first principles only. However, this may change in a future. Indeed, some recent progress in Lattice QCD [6][7][8] gives hope that reliable calculations of the HVP are now in reach in the next years. They would be an important complement to the standard approaches, as well as to the approach presented here.
One is, therefore, left with estimates numerically derived from experimental data. Indeed where K(s) is a known kernel [4] enhancing the weight of the low s region, close to the threshold s H i of the final state H i . Then, the total non-pertubative HVP can be estimated by a μ (H ) = a μ (H i ), where the sum extends over all final states H i which can be reached in e + e − annihilations.
The accuracy of a μ (H i ) is, of course, tightly related with the accuracy of the experimental data set used to perform numerically the integration shown above. When different data sets are available for a given annihilation channel H i , a combination of the corresponding a μ (H i )'s is performed by weighting each estimate with the reported uncertainties affecting each data set, using standard statistical methods (see [9], for instance). Possible mismatches between the various estimates are accounted for by methods like the S-factor technics of the Particle Data Group [10]. In this approach, of course, the accuracy of each a μ (H i ) is solely determined by all the measurements covering the channel H i only, without any regard to the other channels H j (j = i).
This method succeeds in providing very precise values for the relevant contributions. Summing up the nonperturbative HVP estimated this way with the rest, one obtains an estimate of a μ quite comparable to the BNL average measurement [1,2]. However, the prediction based on e + e − annihilation data or τ decay data [4,[11][12][13][14][15][16] exhibits a long-standing discrepancy; the exact value of this discrepancy has gone several times back and forth, depending on whether one trusts the τ data based analyses or the scan e + e − annihilation data, which are obviously more directly related with a μ (H ). With the advent of the high statistics data samples collected using the Initial State Radiation (ISR) method [17][18][19], a precise value for this-possiblediscrepancy has become harder to define unambiguously.
In order to get a firm conclusion concerning the numerical difference between the measured and calculated values of the muon anomalous magnetic moment a μ = (a μ ) exp − (a μ ) th , one should first understand why τ based and e + e − based analyses differ; one should also understand the differences between scan data and ISR data and possibly the differences between the various available ISR data samples, as the KLOE samples [17,19] and the BaBar sample [18] seem to lead to somewhat conflicting results.
Anyway, while all proposed values for (a μ ) th differ from the average for (a μ ) exp , the theoretical uncertainties start to be comparable to the experimental one. Therefore, it becomes interesting to look for a method able to reduce the uncertainty on (a μ ) th by simply using the existing data. It is also an important issue to have a framework where the properties of each data set can be examined.
In order to cover the low energy regime of strong interactions, the most common approach is to use effective Lagrangians which preserve the symmetry properties of QCD. At very low energies, Chiral Perturbation Theory (ChPT) represents such a framework. However, the realm covered by the usual ChPT is very limited (not much greater than the η mass); Resonance Chiral Perturbation Theory (Rχ PT) permits to go much deeper inside the resonance region; it thus defines a framework suited to study the non-perturbative hadronic VP (HVP).
It was soon recognized [20] that the coupling constants occurring at order p 4 in ChPT were saturated by low lying meson resonances of various kinds (vector, axial, scalar, pseudoscalar) as soon as they can contribute. This emphasized the role of the fundamental vector meson nonet and confirmed the relevance of the Vector Meson Dominance (VMD) concept in low energy physics. Soon after, [21] proved that the Hidden Local Symmetry (HLS) model [22,23] and the Resonance Chiral Perturbation Theory (Rχ PT) were equivalent. Therefore, one may think that the HLS model provides a convenient and constraining QCD inspired framework for an improved determination of the HVP. It is, therefore, quite legitimate to check wether the HLS model allows a better determination of the HVP than the usual method sketched above. The basic HLS model has an important limitation for HVP studies: The vector resonances entering the model are only those embodied in the lowest mass vector meson nonet. This certainly limits upwards the relevant energy range to 1.05 GeV, i.e. slightly above the φ(1020) meson mass; going beyond while staying within the standard HLS framework certainly entails uncontrolled uncertainties due to the contribution of higher mass vector meson nonets.
However, relying on the standard method, one can estimate the contribution of the region √ s ∈ [m π 0 , m φ ] to 83.3% of the total HVP and show that its uncertainty is also a large fraction of the total HVP uncertainty: 4 × 10 −10 when using only scan data or 2.7 × 10 −10 when using also the recent ISR data samples. For comparison, the uncertainty provided by the region above 1.05 GeV is 4 × 10 −10 . Therefore, any significant improvement on the knowledge of (a μ ) th in the region √ s ≤ 1.05 GeV is certainly valuable. The (basic) HLS model provides a framework where the interrelations between the various observed decay channels are made explicit. The point is that the use of an adequate model allows for a global fit strategy. All available cross-section data are used to constrain the model parameters, which in turn allows us to predict physical amplitudes. Therefore, if the model provides a statistically acceptable common solution to some set H ≡ {H i } of different processes, 1 each covered by one or several data sets, the fit results can serve to reconstruct reliably the a μ (H i Indeed, if a global fit of the set H of the various data samples is successfully performed, then the parameter values and their error covariance matrix summarize reliably all the knowledge of the set H, including the physics correlations among them. Then, all cross sections contained in H can be estimated with an information improved by having taken into account all the underlying physics correlations. With the present formulation, of the HLS model, the various a μ (H i ) can be reliably and accurately determined up to 1.05 GeV, just including the φ resonance region. All the rest should presently be estimated by the methods usual in this field.
One can substantiate the benefits drawn for using such a global model: • As the model is global, it implies algebraic relations between the various channels it encompasses. Therefore, the accuracy of the estimate for a μ (H i ) is determined by the statistics available for any channel H i and also by the statistics associated with all the other channels contained in H. For instance, the accuracy for a μ (π + π − ) is certainly determined by the available statistics for e + e − → π + π − but all other data, acting as constraints, also contribute to the accuracy for a μ (π + π − ). This is the role of the e + e − → ηγ or e + e − → π 0 γ annihilation data, but also those of the decay width for φ → η γ or of the dipion spectrum in the τ → ππν decay, etc.
Conversely, the accuracy for a μ (π 0 γ ), for instance, is not only governed by the statistics available for e + e − → π 0 γ , but also by those provided by the e + e − → π + π − or e + e − → ηγ data, etc.
Therefore, the improvement expected from a global model should affect simultaneously all the channels contained in H and contributing to a μ .
• As the breaking procedure is global, it affects simultaneously all physics channels related with each other by the Lagrangian model. A successful global fit thus implies that it is validated by the fit quality of the largest possible set of data samples. This high degree of consistency indicates that the breaking model 2 is not simply had hoc. • Any data set is certainly subject to specific systematics; however, taking into account that the study we plan relies on 45 different data sets covering 6 different annihilation channels, 10 partial width decays (taken from the Review of Particle Properties [10]) and some decay spectra, 3 one may consider the effects of correlated systematics reasonably well smeared out. Indeed, one may consider unlikely that the systematics affecting as many different objects can pile up.
Basically, what is proposed is to introduce the theoretical prejudice represented by one formulation of the VMD concept in order to constrain the data beyond genuine statistical consistency of the various data samples referring to the same physics channel. It has already been shown [24] that theoretical (VMD) relationships among various channels are highly constraining. The present work plans to better explore such a framework with a much improved modelling.
Conceptually, the idea to include some theoretical prejudice in order to reduce the uncertainties on a μ is not completely new. A method to complement the e + e − → π + π − data with the constraints of analyticity, unitarity and chiral symmetry has been initiated by [25][26][27] with the aim of improving the π + π − contribution to a μ , but this has not been finalized.
For the present study, we have found appropriate to discard the data collected using the Initial State Radiation (ISR) method [17][18][19]; indeed, because of the complicated structure of their systematics, they almost certainly call for a more complicated statistical treatment than the usual e + e − scan data. The use of ISR data will be addressed in a forthcoming publication.
The HLS model [22,23] complemented with its anomalous sector [28] provides a framework able to encompass a large realm of low energy physics. This anomalous sector will be referred to hereafter as FKTUY sector. The nonanomalous sector allows to cover most e + e − annihilation channels and some τ decays. Thanks to its anomalous sector, the same framework also includes the radiative decays of light flavor mesons with couplings of the form 4 V P γ and P γ γ and also several anomalous annihilation channels. Actually, up to the φ meson mass, the only identified channel which remains outside the scope of the HLS model is the e + e − → π 0 ω annihilation channel, due to the large effect of high mass vector resonances [29,30] presently not included in the HLS model.
However, in order to use the HLS model beyond rather qualitative studies and yield precise descriptions of experimental data, symmetry breaking procedures have to be implemented. A simple mechanism breaking the SU(3) flavor symmetry [31] has been introduced, followed by several useful variants [32][33][34]. Nonet symmetry breaking in the pseudoscalar sector has also been introduced by means of determinant terms [35]. This breaking procedure has been shown to describe precisely the radiative decays of light mesons [36,37] and to meet [38] all expectations of Chiral Perturbation Theory.
In order to account for the reported mismatch between the pion form factor in e + e − annihilation and in the τ decay, it has been proposed [39] to take into account loop effects. Indeed kaon loops produce a mixing of the neutral vector mesons which is a consequence for the K 0 − K ± mass splitting. These turn out to modify effectively the vector meson mass term by identified s-dependent terms.
Introducing the physical vector fields which correspond to the eigenstates of the loop modified vector meson mass matrix, provides a mixing mechanism of the triplet ρ 0 -ω-φ system. In this change of fields the charged vector mesons remain unchanged. With this s-dependent mixing of neutral vector mesons, the fit residuals to the pion form factor in e + e − annihilations and in τ decays did not exhibit any longer any mass dependence [39]; thus this mechanism provides an important part of the solution to the so-called e + e −τ puzzle. 5 However, this solution is only partial. Indeed, if the dipion spectrum lineshape in the decay of the τ lepton is clearly predicted [24,39] from e + e − data, there is still some problem with its absolute magnitude. This issue has been found to be cured by allowing (i) a mass (δm) and a coupling (δg) difference between the neutral and charged ρ mesons, (ii) a rescaling of the τ dipion spectra consistent with the reported uncertainties on the absolute scales of the various measured spectra [40][41][42]. The results returned by fits did not lead to a significant mass difference 6 but, instead, δg and the fitted scales of the experimental spectra were found highly significant [24].
However, the numerical values of these parameters (never more than a few percent) suggest that some unaccounted for isospin breaking effects have not yet been included.
On the other hand, the HLS model supplemented with the SU(3)/U(3) breaking reminded above accounts successfully -and simultaneously-for the measured cross sections in the e + e − → π + π − , e + e − → π 0 γ , e + e − → ηγ , e + e − → π + π − π 0 annihilation channels and for the additional set of 9 decay widths, especially the radiative decays of the form V P γ or P γ γ , needed in order to constrain more tightly the model. This has been proved in [46]. However, as it stands, the HLS model fails to account for the annihilation channels e + e − → K + K − and e + e − → K 0 K 0 simultaneously. This is obviously related to the puzzling issue thoroughly discussed in [47] concerning the branching fraction ratio φ → K + K − /φ → K 0 K 0 . The reported disagreement with theoretical expectations is found significant and amounts to a few percent. This also allows thinking that some isospin breaking effects are not yet fully accounted for. 5 A similar result has been obtained in [16] relying rather on a ρ 0 -γ mixing mechanism; it should be interesting to study a more general V -γ mechanism supplementing the ρ 0 -ω-φ mixing scheme. 6 The mass difference following from fit corresponds to δm = 0.68 ± 0.40 MeV is in accord with what is expected for the electromagnetic mass difference [43,44] of the ρ mesons [45] δm 0.81 MeV.
In the present paper, we define a symmetry breaking procedure which is nothing but an extension of the BKY mechanism referred to above, but including now breaking in the non-strange sectors. This mechanism is only an upgrade of the BKY mechanism and applies likewise to the two different sectors (the so-called L A and L V sectors) of the nonanomalous HLS Lagrangian. We show that the τ scale issue is solved by breaking the L V Lagrangian piece while the φ → K + K − /φ → K 0 K 0 puzzle yields its solution from applying the same mechanism to the L A Lagrangian piece. Stated otherwise, within the framework of the HLS model broken in this way, the e + e − -τ and the φ → KK puzzles appear as twin phenomena yielding parent explanations. Actually, equipped with this upgraded breaking mechanism, the HLS model provides a satisfactory description of all the physics information listed above, including now both e + e − → KK annihilations.
Having discarded the 3 existing ISR data samples, a priori 45 different data sets of scan data are relevant for our present analysis. At each step of our analysis, we have checked the consistency of the various data samples with each other by relying, as strictly as possible, on the information provided by the various groups without any further assumption. We have found that 2 among them have a behavior not in agreement with what can be expected from the rest (43 data sets). One could have attempted to use them by weighting their contribution to the global χ 2 (a sort of S-factor); however, for now, we have preferred discarding them. Therefore our analysis relies on 43 data sets-mostly produced by the CMD-2 and SND Collaborations-and 10 accepted partial width information, which represents already an unusually large set of data consistently examined and satisfactorily understood.
The paper is organized as follows. In Sect. 2, we briefly outline the basics of the HLS model and its various sectors. In Sect. 3, we define the upgraded breaking procedure which is a trivial extension to the u and d sectors of the BKY breaking scheme as redefined in [34]-the so-called "new scheme". Section 4 and Sect. 5 examine the consequence for the modified BKY breaking scheme on the two different parts (L A and L V ) of the non-anomalous HLS Lagrangian. In Sect. 6 we first remind the loop mixing scheme [39,46] of the vector mesons and, next, construct the pion form factor in e + e − annihilation and in the decay of the τ lepton. The condition F π (0) = 1 has some consequences for how parametrizing the Breit-Wigner amplitudes should be done for narrow objects like the ω and φ mesons. Other topics are also examined: the direct ωππ coupling and the φ → KK couplings. The anomalous sector is examined in Sect. 7 where we also provide the expressions for the e + e − → π 0 γ , e + e − → ηγ and e + e − → π + π − π 0 cross sections. The expression for the various couplings of the form P γ γ and V P γ are also derived; these are important ingredients for the set of radiative decays included into the HLS framework.
We have found it appropriate to summarize the main features of the HLS model under the upgraded breaking scheme which underlies the present study; this is the matter of Sect. 8. Section 9 is devoted to list the different data sets available for each physics channel; in this Section, our fitting method, previously defined and used in [24,39,46], is reminded.
At this point, we are in position to confront our model and the data. Section 10 examines the fit properties of the available e + e − → π + π − π 0 data and Sect. 11 reports on the simultaneous analysis of the e + e − → K + K − and e + e − → K 0 K 0 annihilation data. The analysis of the e + e − → KK channels allows us to show how the problem raised by both φ → KK decay widths is solved within the new release of the broken HLS model. Section 12 provides our analysis of the dipion spectrum in the τ decay in conjunction with all e + e − data. It is shown therein that e + e − data and τ data are fully reconciled; the precise mechanism solving this issue, somewhat unexpected, is exhibited.
The short Sect. 13 is devoted to examining the exact structure of the ωππ coupling and compare with similar results of other authors [48,49]. Similarly, another short Sect. 14 examines in some detail some properties of the π 0η-η mixing; it is shown here that the conclusions derived in [38] about the mixing angles θ 0 and θ 8 introduced by [50,51] are confirmed, together with their relationship with the traditional singlet-octet mixing angle θ P . In Sect. 15, one examines the fitted values of the parameters involved in the absolute scale of the FKTUY anomalous Lagrangian pieces and compare with existing estimates; this leads to the conclusion that the usual assumption c 3 = c 4 is consistent with data.
Section 16 is devoted to study in detail the consequences for our HLS model determination of the non-perturbative part of the photon hadronic vacuum polarization. This is found to yield much reduced uncertainties compared to estimates derived by the direct averaging of data.
The consequence for g − 2 are also examined with the conclusion that the theoretical prediction differs from the BNL measurement [2]. The significance of this difference is shown to stay in between 4.07σ and 4.33σ . This looks an important improvement, as we are still not using the ISR data.
Finally Sect. 17 provides a summary of our conclusions and the perspectives. A large part of the formulae have been pushed inside several Appendices in order to ease as much as possible the reading of the main text.

The HLS Lagrangian
The Hidden Local Symmetry Model (HLS) has been presented in full detail in [22] and, more recently, in [23]. One can also find brief accounts in [34,52].
Beside its non-anomalous sector, which allows to address most e + e − annihilation channels and some τ decays up to about the φ meson mass [39,46], the HLS Model also contains an anomalous (FKTUY) sector [28] which provides couplings of the form V V P , V P P P , γ P P P , V P γ or P γ γ among light flavor mesons. These are the key in order to incorporate within the HLS framework the radiative decays of the form V P γ or P → γ γ , or decays importantly influenced by the box anomaly like η/η → π + π − γ (see [37,53] for instance). It has been shown that, while implementing (U(1)) nonet symmetry and SU(3) symmetry breakings, one reaches a remarkable agreement with data [36,37].
The anomalous pieces of the HLS Model are also the key tool when dealing with annihilation processes like e + e − → π 0 γ , e + e − → ηγ or e + e − → π 0 π + π − as successfully shown in [46].
In order to be self-contained, and without going into unnecessary detail, let us briefly remind the salient features of the HLS Model relevant for the present purpose.
One defines the ξ fields by: where the scalar field σ is usually eliminated by means of a suitable gauge choice [22] (the so-called unitary gauge).
However the decay constant f σ still survives through the ratio a = f 2 σ /f 2 π which is a basic (free) ingredient of the HLS Model. The standard VMD Lagrangian corresponds to having a = 2. The pseudoscalar field matrix P : contains singlet (P 0 ) and octet (P 8 ) terms. By π 3 we denote the bare neutral pion field; the traditional naming π 0 will be devoted to the fully renormalized neutral pion field.
On the other hand, the usual η and η meson fields are (essentially) combinations of the η 8 and η 0 fields shown in (2).
The HLS Lagrangian is defined by: where the covariant derivatives are given by: with: exhibiting the Z, W ± boson fields together with the photon field A μ . The vector field matrix is given by: The quark charge matrix Q is standard and the matrix T + = [T − ] † is constructed out of matrix elements of the Cabibbo-Kobayashi-Maskawa matrix [23,39]. One should note that the neutral charge entries of the vector field matrix V are expressed in terms of the so-called ideal fields (ρ I , ω I and φ I ).
In the expressions above, one observes the electric charge e, the universal vector coupling g and the weak coupling g 2 (related with the Fermi constant by g 2 = 2m W G F √ 2). As the influence of the Z boson field is quite negligible in the physics we address, the Weinberg angle θ W plays no role.
We do not present here the anomalous sectors which can be found in the original HLS literature [22,23,28]. A summarized version, well suited to the present purpose, can be found in [46]) and will not be repeated.
The non-anomalous Lagrangian L HLS at lowest order in field derivatives can be found expanded in [34,52]. Its τ sector is explicitly given in [39,46].
The HLS Lagrangian fulfills a U(N f ) × U(N f ) symmetry rather than SU (N f ) × SU (N f ). The additional axial U(1) symmetry has several undesirable features [35,54], especially a ninth light pseudoscalar meson. This symmetry can easily be reduced by adding appropriate terms to the effective Lagrangian. Defining [22] the chiral field U ≡ ξ † L ξ R = exp 2iP /f π , this reduction is achieved by adding determinant terms [35] to the HLS Lagrangian. After this operation, one gets [38]: where μ has obviously a mass dimension and λ is dimensionless. In the following, the additional Lagrangian piece will not be modified while breaking symmetries. Actually, in the present work, one is only concerned by the perturbation of the pseudoscalar meson kinetic energy.

The BKY-BOC breaking of the HLS Lagrangian
The HLS Lagrangian above is certainly an interesting and attractive framework. However, without introducing suitable mechanisms for symmetry breaking effects, one cannot account for the experimental data at the level of precision required by their accuracy. There is no unique way to implement such a mechanism within the HLS model and, actually, several SU(3) breaking schemes exist. The basic SU(3) symmetry breaking scheme has been proposed by Bando, Kugo and Yamawaki (BKY) [31]. It has, however, some undesirable properties which have motivated its modification. A first acceptable modification has been proposed by Bramon, Grau and Pancheri [32,33] and another one in [34], where the various schemes have been critically examined. Following this study, we prefer using in the following the so-called "new scheme" variant defined in [34]; when referring to the BKY mechanism throughout this paper, we always mean the "new scheme" variant just mentioned. It will be referred to as either BKY or BKY-BOC. This breaking mechanism (BKY-BOC) has been examined in detail and its predictions-relying on fits to experimental data-have been found to meet the corresponding ChPT expectations [38] at first order in the breaking parameters. It has also been extensively used in several successful studies performed on radiative decays of light mesons [36] and on e + e − annihilation cross sections [46]. Up to now, the BKY mechanism was limited to SU(3) symmetry breaking effects; the issue now is to examine its extension to isospin symmetry breaking.
Briefly speaking, our variant of the BKY mechanism [34] turns out to define the broken non-anomalous HLS Lagrangian pieces by: where X A and X V are matrices carrying the SU(3) symmetry breaking associated with, respectively, L A and L V . These are written as: and departures of z A and z V from 1 account for SU (3) symmetry breaking effects in the L A and L V Lagrangian pieces. 7 A priori, these two parameters are unrelated and should be treated as independent of each other.
In order to extend to isospin symmetry breaking, we propose to generalize (9) above to: As isospin symmetry breaking is expected milder than SU(3) breaking, the additional breaking parameters are obviously expected to fulfill: In previous fits, performed with only SU(3) symmetry breaking, we got (see for instance [39,46]) |z A − 1|, |z V − 1| 0.5. Such ways of extending the BKY breaking mechanism have been already proposed within similar contexts [55,56]. We find appropriate to define: exhibiting the sum and difference of u A/V and d A/V . Indeed, the expressions for most physical couplings are simpler in terms of these rather than in terms of u A/V and d A/V . As clear from (8), the BKY breaking of the HLS Lagrangian exhibits a global character. It does not correspond to some systematic way of including specific breaking terms of given kind or order as done within ChPT. As the numerical values of the breaking parameters are phenomenologically derived from fits to a large set of experimental data, they account globally for several effects of different origin without any way to disentangle the various contributions. This remark is especially relevant for the breaking parameters corresponding to the u and d entries of X A and X V which are small enough that several competing effects can mix together 8 ; because of their relatively large magnitude, the SU(3) breaking effects can be more easily identified [36,38].

Breaking the L A Lagrangian piece
The pseudoscalar kinetic energy term of the full (broken) Lagrangian is carried by L A + L tHooft . In terms of bare fields, it is: 7 In the following X A and X V are named breaking matrices; this convenient naming should not hide that the true breaking matrices are rather X A − 1 and X V − 1. 8 This may include effects due to the quark mass breaking and to electromagnetic corrections. It may also absorb corrections to hadronic vertices which can hardly be derived from an effective Lagrangian.
which is clearly non-canonical. In order to restore the canonical structure, one should perform a change of fields. This is done in two steps, as in [38].

First step PS field renormalization
The first step renormalization turns out to define the (step one) renormalized pseudoscalar field matrix P R 1 in term of the bare one P by: The charged pion and both kaon terms in this expression actually undergo their final renormalization already at this (first) step. Indeed, the axial currents are defined as in [38] and are given by: in terms of the bare fields and of the Gell-Mann matrices T a , normalized such as Tr[T a T b ] = δ a b /2. The π ± , K ± and K 0 /K 0 decay constants are defined by the relevant axial current matrix elements closed on the renormalized PS meson fields: 0|J π ± /K μ |π ± R = if π/K q μ . As one chooses the renormalized (charged) pion decay constant to coincide with its experimental central value [57] (f π ± = 92.42 MeV), this turns out to impose q A y A = 1. At leading order in breaking parameters, this implies Σ A = 0. Then, the breaking matrix X A writes: depending on only two free parameters ( A and z A ). On the other hand, the kaon decay constant is: One thus yields a marginal change compared to the previous BKY breaking scheme [34,38] (dealing only with SU(3) symmetry) as one got [f K ± /f π ± ] 2 = z A . Anticipating somewhat on our numerical results, let us mention that the fits always return A (5. ÷ 6.)%, much larger than expected from solely an effect of the light quark mass difference [58]; this will be further discussed in Subsect. 11.2. As clear from (16), the X A matrix resembles the usual quark mass breaking matrix. However, the entry z A is essentially related with the ratio [f K /f π ] 2 1.5, while the corresponding entry in the quark mass breaking matrix is numerically 20. Therefore, the correspondence between these two matrices is only formal.
The following relationship defines some bare PS fields in terms of their (fully) renormalized partners: This produces changes going in opposite directions for the couplings involving the physical K ± and K 0 mesons compared to their bare partners. This has a clear influence on the cross section ratio σ (e + e − → K + K − )/σ (e + e − → K 0 K 0 ). On the other hand, one also gets at leading order in breaking parameters the following relationship between some bare PS fields and the (first step) renormalized PS fields: and is independent of A . As we get-at leading order-the same dependence as before [38], the diagonalization procedure for this term is known (see Sect. 3.1 in [38]). Let us only recall the results in the present notation set: where: It thus looks more appropriate to use v ( λ/2) rather than [38] λ as a breaking parameter, as it allows to work with simpler expressions. v is the first parameter in our model which exhibits the intricacy between U(1) and SU(3) breakings (λ and z A ). The canonical PS fields-denoted by the superscript R-are finally defined by (18), (19), (21).

4.3
The π 0 -η-η mixing The physically observed η and η are traditionally described as mixtures of the singlet and octet PS fields η 8 and η 0 involving one mixing angle named here θ P . Some authors, following [54,59] prefer now using mixtures of the uū + dd and ss wave functions. However, as these two approaches are equivalent, we prefer sticking to the traditional description.
Since [50,51], it is admitted that the most appropriate ChPT description of the η-η mixing involves two decay constants (F 0 and F 8 ) and two mixing angles (θ 0 and θ 8 ). However, [38] has shown how, within VMD, the usual octet-singlet mixing scheme connects with this new approach. In this reference, it was also shown that, relying on experimental data, the broken HLS model favors θ 0 = 0 with a very good accuracy; this led to a relation between θ 8 and θ P numerically close to θ 8 = 2θ P . Comparing accepted ChPT numerical values for θ 8 -like those in [51]with the one derived from θ P (determined in VMD fits) was found quite satisfactory. Moreover, it was shown that fits to experimental data lead to an algebraic relation of the form θ P = f (λ, z A ). We will check whether this relation for θ P still fits within the present form of our broken HLS model.
As in all previous studies in this series, one could have limited oneself to considering only the η-η mixing, decoupling this from the π 0 sector. However, it is a classical matter of Chiral Perturbation Theory to address the issue of π 0η mixing, as this is related with the (light) quark mass difference [60]. Therefore, it may look interesting to see if such a phenomenon could be exhibited from the experimental data we deal with. In this case, there is no reason not to address the issue of the relevance of a full π 0 -η-η mixing mechanism. We choose to parametrize this PS mixing using [61]: where π 3 R , η 8 R and η 0 R are the already redefined fields (see (21) above), the physically observable mesons being π 0 , η and η . In the smooth limit of vanishing and , one recovers the usual η-η mixing pattern with one (θ P ) mixing angle, while the pion field decouples. Even if one does not expect a large influence of and in the full data set collection we consider, it does not harm to examine their effects and, if relevant impose = = 0 to the model.
Finally, at leading order in breaking parameters, the pseudoscalar meson kinetic energy term is canonical when expressed in terms of the fully renormalized fields (those carrying a R subscript).
4.4 About the θ 8 , θ 0 and θ P mixing angles Let us define: These functions tend to unity when the U A (1) symmetry is restored (λ = 0). The property g 8 (v) 1 is the simplest way to justify the approximation done in our previous works to parametrize nonet symmetry breaking by the parameter x (see, for instance, the discussion in [38]). Using these functions, one can derive from (15) the following axial currents: The mixing angles θ 8 , θ 0 [50,51] yield the following expressions: One can easily check that g 8 (v)/g 0 (v) 1 − λ/2. A property to check from fits using the present form of the model is whether θ 0 is still consistent with zero [38]. In this case, θ P is still no longer a free parameter, but fully determined by λ and z A , i.e. by breaking parameters and θ P tends to zero when the symmetries are restored.
One should also note that the usual ChPT definition of the π 0 decay constant ( 0|J π 3 μ |π 0 = iq μ f π 0 ) provides f π 0 = f π ± , not influenced by our isospin breaking procedure. However, as will be seen shortly-and as already emphasized in [38] for the decays η/η → γ γ -this is not the quantity actually involved in the decay π 0 → γ γ .

Breaking the L V Lagrangian piece
The L V Lagrangian, is defined by (8)- (12). It yields the following vector meson mass term (m 2 ≡ ag 2 f 2 π ): while keeping only the leading terms in the breaking parameters Σ V and V (the K * mass term has been dropped out). As can be seen, the canonical structure of the mass term is broken by a ρÎ · ωÎ term. In order to restore the canonical form of the mass term, one should perform a field redefinition in only the (ρÎ − ωÎ ) sector. Interestingly, the requested transform is not a rotation but: makes the work when non-leading terms in Σ V and V are neglected. 9 In terms of the R 1 renormalized fields, one gets: having renamed for convenience φ I ≡ φ R 1 . A few remarks are worth being done: (1) Beside the two breaking parameters Σ V and V , one gets a third free parameter h V which governs the mixing of ρ R 1 and ω R 1 . Therefore, the exact content of isospin 1 (ρ I ) inside ω R 1 and of isospin 0 (ω I ) inside ρ R 1 should be extracted from data.
(2) The masses for ρ R 1 , ω R 1 and ρ ± remain degenerated at leading order in the breaking parameters as the needed R 1 change of fields results in a vector meson mass term independent of V . Therefore, even if one may legitimately think that breaking isospin symmetry inside L V could result into a non zero (Lagrangian) mass difference [24] δm 2 between the charged and the neutral ρ mesons, our breaking procedure rules out such a possibility at leading order in breaking parameters. Actually, electromagnetic corrections [43,44], presently neglected, generate such a mass difference. Such a term has been considered in [24] but found numerically insignificant; preliminary studies within the present framework leading to the same conclusion, we have given up considering explicitly a ρ 0 − ρ ± mass difference.
(3) The field transform (29) propagates to the vector meson kinetic energy by generating a term of the form V ∂ρ 0 R 1 ∂ω R 1 which breaks the canonical structure of the kinetic energy. This is a classical issue [62] known to imply the occurrence of wave-function renormalization factors [62] which are absorbed into the effective couplings defined by the Lagrangian vertices. In our case, they are certainly absorbed in our breaking parameters. This is exactly the same issue which arises in the electroweak Standard Model with the γ -Z 0 mixing. This has been investigated in detail within Z 0 lineshape studies (at the one plus two-loop level) and by the LEP experiments. The same issue also appears when treating the γ -ρ 0 mixing and has been discussed in [16].
The second step renormalization of the vector meson fields, which accounts for loop effects [39,46], is considered below.

The fully broken non-anomalous HLS Lagrangian
For definiteness, we name (abusively) from now on "HLS Lagrangian" the full expression given in (7), i.e. including the determinant terms. The HLS Lagrangian is explicitly provided in Appendix A, dropping out for conciseness all terms not relevant for the purpose of the present study.
At the present step-which does not still include the (second step) redefinition of the neutral vector fields [39,46]several remarks are worth being done: • The vector meson masses occurring in the Lagrangian fulfill m 2 ρ 0 = m 2 ρ ± = m 2 ω . Thus, no mass splitting is generated, except for the φ meson.
• The couplings ρππ undergo isospin breaking (Σ V ) but remain strictly identical for the charged and neutral ρ mesons. Instead, a direct ωππ coupling is generated; it is proportional to (1 − h V ) V . • The ρ 0 -γ and ρ ± -W ± transition amplitudes 10 may slightly differ, as h V V /3 should not exceed a few percent level. 10 Compare f ργ and f ρW as given by (93) and (94).
Therefore, non-vanishing δm 2 = m 2 ρ 0 −m 2 ρ ± and δg ρππ = g ρ 0 π + π − − g ρ ± π ± π 0 , as stated in [24], are not derived by extending the X A /X V breaking scheme to include isospin symmetry breaking. 11 Therefore, non-vanishing δm 2 and δg ρππ are not the way followed by the (broken) HLS model in order to account for the (slightly) different normalizations of the pion form factor in τ decays and in e + e − annihilations. The actual mechanism at work is emphasized below.

Loop mixing of vector mesons
As remarked in [39], pseudoscalar loops modify the vector mass matrix by s-dependent terms. In this way, the ρ, ω and φ squared masses become s-dependent through contributions at real s of analytic functions, 12 namely the KK loops and, for the ρ, also the charged pion loop. Conceptually, this turns out to remark that the inverse vector meson propaga- More important, however, is that this s-dependent mass matrix becomes non-diagonal, showing that, at one-loop order, the ρ, ω and φ (corresponding here to the R 1 renormalized vector fields) are no-longer mass eigenstates. Mass eigenstates can easily be constructed by standard perturbative methods [63] as shown in [39]; one observes that they become s-dependent.
This mass matrix can be written: where 13 : is treated as the unperturbed part of the squared mass matrix. The pion loop is weighted by the square of the ρ R 1 ππ coupling constant (see (92) in Appendix A) and has been included in the ρ R 1 entry as Π ππ (s) is not really small in the timelike region. Instead, as the ω R 1 ππ coupling is first order in V , the pion loop contribution to the ω R 1 entry should be neglected ( O( 2 V )). The values for these (Higgs-Kibble) masses can be found in (93); they fulfill m ρ = m ω . On the other hand, δM 2 (s) is given by: 11 As stated above, electromagnetic corrections contribute to generate a non-vanishing δm 2 without, however, a significant influence on the fit properties. 12 Actually, the anomalous FKTUY Lagrangian and the Yang-Mills terms contribute respectively with V P and V V loops; one can consider their influence absorbed in the subtraction polynomials of the P P loops [39]. 13 Entries are ordered respectively ρ R 1 , ω R 1 and φ R 1 . and contains only the perturbations generated by kaon loop effects. The kaon loop transition from a vector meson V to another one V has been denoted V V . One should note [39] that M 2 is an analytic function of s satisfying the (so-called) hermitian analyticity condition: The entries of these matrices are appropriately parametrized in terms of: where Π (s) denotes the amputated pion loop, while Π + (s) and Π 0 (s) denote, respectively, the amputated charged and neutral kaon loops; their analytic expressions can be found in the Appendices of [39]. 1 (s) 2 (s) do not contain symmetry breaking terms beyond the effects of the kaon mass splitting. The expressions for the entries in δM 2 (s) are given in Appendix B and show this dependence explicitly (see (97)). One can construct, as in [39], the eigensolutions of M 2 . These are the final (step two) renormalized vector fields denoted respectively by ρ R , ω R and φ R and are related with their R 1 partners by: where the s-dependent mixing angles are defined by: using the eigenvalues of M 2 (at first order): The ρ (s), ω (s) and φ (s) quantities, defined in (97), only depend on the kaon loop functions and on breaking parameters.
6.2 The pion form factor in e + e − annihilations and in τ decays The pion form factor in the τ ± decay to π ± π 0 ν τ can easily be derived from the Lagrangian piece L τ given in (94): where: (38) and the loop functions are: where π (s) and K (s) denote respectively the amputated charged pion and kaon loops, P W (s) and P ρ (s) being subtraction polynomials. In order to fulfill current conservation, these polynomials should vanish at s = 0. Here, as in former studies [24,39,46], identifying P ± P ∓ and P ± P 0 loops has been found numerically justified. If one compares with the corresponding formulae in [24] (Subsect. 2.1.1), one sees that δm 2 and δg-supposed to reflect different properties of the charged and neutral ρ mesons-have been deleted. As the loop functions vanish at s = 0, one clearly has F τ π (0) = 1. The pion form factor in e + e − annihilations is not as simply derived. One needs first to propagate the transformation in (34) into the Lagrangian equation (92) and collect all contributions to, respectively, ρ R , ω R and φ R . In this way, the V -γ couplings associated with the fully renormalized vector fields become: including the mixing angle contributions. Using the Lagrangian pieces given in (99), one can construct easily the pion form factor: where: The loop corrected V -γ transitions amplitudes F e V γ (s) are defined by: with the s-dependent loop terms Π V γ (s) being defined in Appendix C. All Π V γ (s) are requested to vanish at s = 0 because of current conservation. The inverse ρ propagator D ρ (s) is defined by (see (36)): As the ρ self-mass Π ρρ (s) vanishes at s = 0, one certainly has D ρ 0 (0) = −m 2 ρ . Concerning the ω and φ mesons, one can correspondingly write their inverse propagators as: and one can legitimately assume their self-energies to also vanish at s = 0. Then, should certainly be fulfilled. However (most of) the ω selfenergy cannot be computed in closed form and the 3-pion part of φ self-energy too. Therefore, convenient forms for their propagators should be considered. This issue is readdressed just below for both mesons. At this step, it is of concern to compare the properties of the isospin 1 part of F e π (s) with F τ π (s). The most important pieces of information are listed in Table 1. The difference displayed for the non-resonant term is tiny. One can see that there is no mass difference between the charged and neutral ρ mesons, nor different couplings to a pion pair. Instead, most of the difference is actually carried out by the transition amplitudes (see the fifth data line in Table 1) which are significantly s-dependent, as can be inferred from Figs. 6 and 7 in [46].
Finally, it is interesting to note that the renormalization factor introduced in couplings involving a kaon pair plays in opposite directions for charged and neutral kaon pairs. 6.3 The ωππ direct coupling and the condition F e π (0) = 1 As can be seen from (42), the fully broken HLS model reveals a total coupling of the ω to a pion pair given by: This expression illustrates that the ωππ coupling in our model is a priori a superposition of a direct isospin breaking term and of another piece generated by vector meson mixing through kaon loops. This kind of sharing has been emphasized several times [48,49]. The full data set we use should give the most precise and motivated estimate for these two pieces as this is still presently controversial [48,49,64].
On the other hand, the parametrization of the ω contribution to the pion form factor may pose a conceptual problem related with the condition F e π (0) = 1 which is worth addressing.
The pseudoscalar meson loops which enter the V V transition amplitudes (see (32), (35) and (97)) behave as O(s) near the origin. The running vector meson masses (see (36)) are such that λ ρ (s) − λ ω (s) vanishes at the origin, while the two other differences which come into (36) tend to a nonzero constant. Therefore, ab initio, the mixing angles are expected to fulfill: Even if clear in the previous publications using the loop mixing mechanism (Fig. 7 in [39] or Fig. 6 in [46] clearly show that α(0) −5%), this was not explicitly pointed out. Therefore, the s-dependent ωππ coupling generated by loop mixing 14 does not vanish at the origin. This has some consequences.
Indeed, using (40) and (42), and the vanishing properties of the functions Π V γ (s), Π V V (s) and β(s) at the origin, one gets: when keeping only the first-order terms in breaking parameters.
As already discussed at the end of the previous Subsection, it is motivated to assume the ω self-energy Π ωω (s) vanishing at the origin. Moreover, this allows to stay consistent with the so-called "Node theorem" [65,66]. Then, ρ . This provides the vanishing of the last bracket in the formula above and, thus, F e π (0) = 1, whatever the values for h V , V and α(0).
However, in most applications, for objects carrying such a narrow width as the ω and φ mesons, one generally uses approximate inverse propagators, e.g. either 15 : with values form V andΓ V either taken from the Review of Particle Properties or extracted from one's fits. Then, with either of these Breit-Wigner lineshapes, the condition F e π (0) = 1 is not necessarily fulfilled. From our model results, this condition is even violated at a few percent level. However, it is easy to check that either of: (remember that m 2 ρ = m 2 ω ) and: Table 1 Comparison of the pion form factor information in τ decay and in e + e − annihilation. Second column lists only isospin 1 related information. In the last entry of the rightmost data column, the upper sign refers to K + K − pairs, the lower to certainly cures this disease. This turns out to parametrize the ω self-energy Π ωω (s) with an ansatz which satisfies its vanishing at the origin. It is worth stressing that using standard Breit-Wigner lineshapes or their modified partners provides practically unchanged fit results. This is due to the fact that the ρ and ω masses (with tilde or not) are close to each other, and then, the factor s/m 2 ρ is very well approximated by 1 all along the sensitive region of the ω peak.
In order to substantiate the possible changes, we have run our code using BW a and BW a as inverse ω propagators. As a typical example of modification, one can comparem ω = 782.44 ± 0.06 MeV andΓ ω = 8.46 ± 0.09 MeV while fitting with BW a , andm ω = 782.49 ± 0.06 MeV andΓ ω = 8.36 ± 0.08 MeV when using instead BW a . For definiteness, in the fits presented in this paper, D ω (s) will be modified as just explained. As β(0) = 0, the pion form factor value at s = 0 is not sensitive to how the φ propagator is approximated.
Even if our choice is motivated, others are certainly possible as exemplified in [13,49]. Transposed to our model, the just mentioned choice would turn out to weight the full ω contribution to the pion form factor by a factor s/m 2 ω or s/m 2 ρ which restores F e π (0) = 1. The behavior of this choice is identical to ours, basically becausem ω and m ρ are very close to each other.

The charged and neutral kaon form factors
We give here the annihilation cross sections/form factors within the extended BKY-BOC breaking of the HLS Lagrangian. Cross sections and form factors are related through: for any meson pair P P . q P = s − 4m 2 P /2 is the P momentum in the center-of-mass system. The kaon form factors are given by: where the γ -V transition amplitudes F V γ have been already defined (see (43)). The V KK couplings can be read off from the corresponding Lagrangian pieces ( (100) and (101)).
The kaon form factors fulfill: However, it is easy to check that these conditions are both fulfilled, only if: Therefore a fixed width Breit-Wigner shape for the φ should be adapted as already discussed for the ω.

Parametrization of the φ propagator
As for the pion form factor, in order to fulfill F e K c (0) = 1 and F e K 0 (0) = 0, one should impose that the ω and φ inverse propagators at s = 0 are equal in magnitude and opposite in sign to their respective Lagrangian masses (m 2 Here again, this turns out to parametrize the full selfenergy Π φφ (s) by an ansatz vanishing at s = 0. For the twobody loops, this is well known [39]; however, the three-body loop is not known in closed form (as for Π ωω (s)).
However, in contrast with the case for ω, using: for the φ inverse propagator, instead of the usual (fixed width) BW a form, should be further commented, as √ z V m ρ significantly differs fromm φ normally fitted 16 (e.g. with BW a ).
Even if anticipating on our fit results, it is worth discussing this matter right now. As far as cross sections are concerned, the two kinds of fits provide almost identical results. In order to yield this result, almost all parameters vary within errors except for 17 z V , which could have been expected. However, it will be shown that this change has a marginal influence on all information of physics importance. Anyway, such kind of information is interesting as it provides a hint on the model dependence of numerical results. Therefore, it has been of concern to compare results obtained with either of BW a and BW a , when appropriate.
Before closing this Section, one may note that, at the φ peak location ( √ s 1020 Mev), the modified Breit-Wigner lineshape provides: which explains why the fit remains successful when using BW a and also why z V should change correspondingly, taking into account that m 2 ρ cannot much vary. The fit quality of the e + e − → KK cross sections will illustrate the validity of this parametrization of the φ propagator.

The Coulomb interaction factor
Beyond modelling, there is an important issue to discuss when dealing with the charged kaon form factor. In the decay φ → K + K − , and more generally as close to the KK threshold, one has to take into account the Coulomb interaction among the emerging charged kaons. This has been first addressed in [67] and recently readdressed (and corrected) in [47]. The net result of this effect is to multiply the charged kaon cross section by the Coulomb factor 18 : In [68], and later in [69], the cross section for charged kaons is multiplied by Z(s)/Z(m 2 φ ). This turns out to consider the Coulomb interaction as a breaking mechanism which affects the charged kaon sector and not the neutral 16 More substantially, with appropriate fits, one yields √ z V m ρ 925 MeV, while a direct fit yieldsm φ 1020 MeV! 17 In fits with BW a for the φ meson, one gets z V = 1.368±0.005, while with BW a the fit returns z V = 1.472 ± 0.001. 18 Actually, the full electromagnetic correction factor is more complicated, but the main effect comes from the Coulomb factor. One assumes that the kaon data which have been submitted to fit have been appropriately corrected for soft photon corrections, which allows to cancel out the term named C i in [47]. one; as the corresponding φ branching fractions are fit independently, this should not affect their results. One may just have to remark that this turns out to incorporate the Coulomb effects inside the corresponding estimates for the φ → K + K − branching fraction.
Up to well defined phase space factors generated by the kaon mass splitting, the partial width ratio φ → K + K − /φ → K 0 K 0 is the square of the corresponding s-dependent effective coupling ratio. Neglecting for each coupling corrections terms of order greater than 1, one can derive from (100) and (101): where the last equation follows from remarking (see Fig. 7 in [46]) that the mixing angle β(s)-defined by (35)-is negligibly small compared to √ 2z V in the φ mass region. Therefore, this mechanism proposes a way for this ratio to depart from unity.
In their throughout study of the φ → K + K − /φ → K 0 K 0 ratio, the authors of [47] examined this issue using several other mechanisms than this one and concluded that none of them was able to accomodate a coupling constant ratio smaller than one (in absolute magnitude). The global fit, based on the suitably broken HLS model, provides a new approach. In this framework, the determination of A is constrained by both e + e − → KK annihilation cross sections separately, and by some more light meson anomalous decays, which also depend on A .

The HLS anomalous sector
In order to treat radiative decays, i.e. the V P γ couplings, and some important annihilation channels (namely e + e − → π 0 γ , e + e − → ηγ and e + e − → π 0 π + π − ) within the HLS framework, one needs to incorporate the appropriate Lagrangian pieces. These are given by the Wess-Zumino-Witten (WZW) terms [70,71] which traditionally account for the triangle (AAP ) and box (AP P P ) anomalies, together with the FKTUY Lagrangian pieces [23,28]: where the four c i are constants left unconstrained by theory [28]. A closer examination of the FKTUY Lagrangian allows to identify five different pieces-listed in Appendix D-and one then remarks that the accessible physics is sensitive to the difference c 1 − c 2 and not to each of them separately. One is then left a priori with three unconstrained parameters [23]. When no breaking is at work, the amplitudes for the couplings 19 P 0 γ γ and P 0 π + π − γ at the chiral pointcomputed within the FKTUY-HLS framework 20 -coincide with those directly derived from the WZW piece in isolation [46]. Due to a sign error 21 in the FKTUY Lagrangian piece L AV P , it was asserted in [46] that the constraint c 3 = c 4 was mandatory in order to recover this property. Actually, this property is automatically satisfied [22,23]. In addition, we have verified that this property is maintained within our fully broken HLS model.
However, the condition c 3 = c 4 , which is fulfilled by VMD models [23] is successful and only turns out to reduce the freedom in fits. Nevertheless, one has examined relaxing this condition and found that our fit results are well compatible with the constraint c 3 = c 4 .

Breaking the anomalous HLS Lagrangian
At this step, the anomalous HLS Lagrangian can be written: with pieces listed in Appendix D. As for the non-anomalous HLS Lagrangian, each among these pieces may undergo specific symmetry breaking independently of each other. This may lead to plenty of free parameters as illustrated by M. Hashimoto [56] who implemented combined SU (3) and Isospin symmetry breakings in the anomalous sector.
A simpler mechanism has also been proposed for SU(3) breaking by Bramon, Grau and Pancheri [32,33]; however, this was insufficient to account for both K * [±,0] → K [±,0] γ decay widths. In [36,39] it was proposed to supplement it with a breaking of the vector field matrix resembling a vector field redefinition. Quite unexpectedly, this provides a (successful) parametrization for the K * radiative partial widths identical to those proposed by G. Morpurgo [72] within a completely different context. Interestingly, this combined mechanism leaves totally unaffected the other sectors of the L V V P piece we deal with; this is well accepted by all data considered [36,39]. This combined breaking mechanism has been studied in detail [46] for all pieces of L anom. with similar conclusions. 19 Here and in the following P 0 denotes either of the π 0 , η and η mesons. 20 E.g. using (55) and the V -γ transitions provided by the nonanomalous HLS Lagrangian. 21 We gratefully acknowledge B. Kubis (HISKP, Bonn University) for having kindly pointed out the issue.
The combined breaking mechanism, as presented in [46], has been examined by combining SU (3) and Isospin symmetry breakings using the complete data set discussed below within the minimization code underlying the present study. It was concluded that possible Isospin symmetry breaking effects-not propagated from the field redefinitions provided by non-anomalous HLS Lagrangian breakingprovide invisible effects. It was then decided to neglect this additional possible source of Isospin symmetry breaking, as the parameter freedom it gives is found useless.
Therefore, for sake of clarity, one only quotes the specific forms for the decay amplitudes K * [±,0] → K [±,0] γ , referring the interested reader to [46] for more information.
As a summary, our dealing with the anomalous sectorexcept for the limited K * sector-involves only 3 parameters: c 1 − c 2 and c 3 and c 4 ; former studies [24,37,46] remain valid, as the condition c 4 − c 3 = 0 is well accepted by the data, as will be shown shortly.

Radiative couplings
For what concerns the radiative decays of light mesons and the e + e − → P γ annihilation processes, one needs L AAP and an effective piece named L AV P defined below.
In terms of the final renormalized pseudoscalar fields and assuming the π 0 -η-η mixing defined in Sect. 4, one can write: At leading order in breaking parameters, the coefficients g P 0 γ γ are given by 22 : another way to define the neutral pion decay constant. The other equations also illustrate that the so-called octet and singlet decay constants as derivable from there have little to do with the standardly defined ones, i.e. from the currents in (25). This question has raised some confusion which motivated the study in [38].
In order to treat the e + e − → π 0 π + π − annihilation process the part of the L AP P P Lagrangian describing the socalled box anomalies is needed. This can be written: with: Equations (58) and (60) show how the triangle and box anomaly amplitudes behave under isospin, SU(3) and PS nonet symmetry breakings. One should especially note the intricacy of SU(3) and PS nonet symmetry breakings.
In order to derive the radiative decay couplings, an effective Lagrangian has been built up from L V V P and the non-anomalous Lagrangian in the same way as in [46]. This can be written in terms of the renormalized R 1 fields: The expression for the various coupling constants g V P γ can be found in Appendix E. In order to derive the physical couplings, one should first apply the transformation given in (34) and then collect the various contributions to each of the (neutral) ρ R , ω R and φ R . Concerning the AV P couplings, it is quite interesting to compare the expressions in (108)- (110) with the corresponding ones in [38,39], derived using an approximate expression for nonet symmetry breaking 23 (the x parameter in the quoted papers). Indeed, the three variants by which nonet breaking occurs (see (107)) are close together and can reasonably well approximated by x eff = 1 − v 1 − λ/2. 23 In order to restore the condition c 3 = c 4 , one should simply make in [39] the replacement c 3 → (c 3 + c 4 )/2.

Breaking the V V P and V P P P anomalous Lagrangians
The V P P P anomalous Lagrangian is given by: where one has limited oneself to display the V P 0 π + π − sector. The leading terms of the couplings occurring in this expression are given in Appendix F. The L V V P Lagrangian piece plays an important role in the annihilation process e + e − → π 0 π + π − . Its relevant part is: where: When going from R 1 -renormalized to the fully renormalized vector fields R, one has to take some care with attributing the s-dependence between the two neutral fields of each monomial in the last two lines of (63). This should be tracked for each R 1 field while applying (34).

7.4
The e + e − → P 0 γ annihilation cross sections Using the Lagrangian pieces given above, the transition amplitudes γ * → P γ can be written similarly to [46]: where Y = −α em N c /πf π has been factored out. q is the incoming photon momentum (q 2 = s), p the outgoing photon momentum (p 2 = 0) and N c = 3. The pieces provided by L AAP are 24 : using the g P 0 γ γ couplings defined in (58), where the (1 − c 4 ) has been factored out. The resonance contributions are gathered in K P 0 (s): where the H P 0 V i (s)-given in Appendix G-are the resonance couplings to P 0 γ and the F V R i γ (s) are the V -γ transition amplitudes defined in (43). The D V i (s) are the vector mesons inverse propagators already encountered. The cross sections can then be written: 7.5 The e + e − → π 0 π + π − annihilation cross section Following as closely as before the notations in [46], the amplitude for the γ * → π + π − π 0 is given by: where ε μ (q) (q 2 = s) is the (heavy) photon polarization vector. T sym is the symmetric part of the amplitude (in terms of the ρπ 'final' states), while T brk (denoted T ρ in [46]) breaks this symmetry. We have found appropriate to introduce separately the contribution T AV P (s) to the full amplitude gener- 24 The corresponding expressions given in [46] carry a missprint: Each of the right-hand sides of (41) is missing a factor of 2.
ated by the L AV P Lagrangian piece (see (105)); its first term is symmetric in terms of the ρπ 'final' states. One has 25 : where all parameters and functions have been already defined, except for the N i (s) functions which are given and commented in Appendix H. One has kept as much as possible the notations used in [46] in order to exhibit the effects of our additional isospin symmetry breaking effects by simple inspection. Finally, T AV P (s) identically vanishes when c 4 = c 3 .
The differential cross section writes: using the (x and y) parametrization proposed by E. Kuraev and Z. Siligadze [73] who provided the kernel function G(x, y) reminded in Appendix H. Note also that each of T sym (s), T brk (s) and T AV P (s) also depend on x and y.

Ugraded breaking of the HLS model: a summary
In the former studies performed along the present lines [24,39,46], roughly speaking, one incorporated nonet symmetry and SU(3) symmetry breaking in the pseudoscalar (PS) sector. In the vector meson sector, only SU(3) symmetry breaking was considered. However, some important effects can be already attributed to isospin breaking effects in the PS sector. Indeed, it is the non-vanishing character of the mixing "angles" α(s) and β(s) which induces s-dependent ρ-ω and ρ-φ mixings at the one loop level. This non-vanishing of the α(s) and β(s) functions proceeds from the kaon mass splitting which breaks the symmetry between the neutral and charged kaon loops and, then, allows to choose the analytic function 1 (s) as non-identically vanishing. Therefore, except for the ω-φ system which would mix anyway at one loop, the full loop mixing mechanism for vector mesons is the prominent consequence for this limited account of isospin breaking. 26 This quite limited breaking scheme, allows already for a good account [24,39,46] of the available data. However, within the realm accessible to the HLS model, two experimental issues remain unsolved: (i) The dipion spectrum lineshape in τ decays is consistent with expectations from e + e − annihilations [39,46], but not its absolute scale [24].
is found inconsistent with all reported expectations [47].
Obviously, this inconsistency propagates to the corresponding e + e − annihilation cross sections.
The first topic has been shown to get a satisfactorybut not perfect-solution by allowing some difference between ρ 0 and ρ ± meson properties to be fitted from data. If the effect of a non-vanishing δm 2 = m 2 ρ 0 − m 2 ρ ± was found small, those generated by a non-vanishing δg = g ρ 0 π + π − − g ρ ± π 0 π ± was found especially significant [24]. Moreover, some rescaling of the τ spectra, consistent with the reported experimental scale uncertainties remained unavoidable.
The second topic is experimentally addressed by considering [68,69] that the Coulomb interaction 27 plays as a symmetry breaking mechanism which modifies the SU(3) relationship g φK + K − = g φK 0 K 0 between coupling constants to . This approach, which turns out to consider the Coulomb interaction as some breaking effect, may look unsatisfactory; anyway, it does not fit with our breaking scheme.
These two issues motivated an upgrade of the breaking scheme of the HLS model in order to check whether an acceptable solution can be derived. The extension to isospin breaking of the BKY-BOC breaking mechanism is a priori an obvious candidate to examine. This has been done in the preceding Sections with several interesting conclusions, which can be summarized as follows: (j) One does not find any signal for a mass or a coupling difference between the ρ 0 and ρ ± mesons. 28 However the coupling difference between ρ-γ and ρ-W might be enforced with respect to [24,39,46] if the break-ing parameter product h V V is found significantly nonzero (see Table 1), (jj) Everything goes as if the universal coupling g remains unchanged in the anomalous sector, while one observes that g is effectively modified to g(1 + Σ V ) for the whole non-anomalous sector. Therefore, isospin breaking in the HLS model generates some mild disconnection between anomalous and non-anomalous processes which needs to be explored.
is found subject to isospin breaking in a novel way compared with the various possibilities examined in [47], Topics (j) and (jj) are both important for scale issues. Indeed, by disconnecting somewhat more than before the ratio of transition amplitudes ρ-γ and ρ-W , one allows the HLS model to get more freedom for the purpose to account for scale issues. More important, both τ and e + e − physics share the same universal coupling (g(1 + Σ V )), but it is no longer common with the scale of the anomalous processes which remains governed by g. Moreover, none among the anomalous couplings, all displayed in several of the Appendices, exhibits a dependence upon Σ V . Stated otherwise, the anomalous couplings-which fix the scales of the anomalous meson decay and annihilation processes-no longer constrain the non-anomalous process scales as sharply as formerly assumed [24,39,46].
Concerning the topics (ii) and (jjj), it should be stressed that the parameter A governing the change of this ratio is not involved only in the ratio. Indeed, each of the e + e − → K + K − and e + e − → K 0 K 0 cross sections should keep valid absolute scales separately. Moreover, as clear from Appendices E, F, G and H, and from (58) and (60) given above, this change of scale should also fit with all anomalous processes, including the π 0 → γ γ partial width, now within the partial width data sample submitted to the global fit. Before ending up this Section and this Part, let us remark that the upgraded breaking of the HLS model allows to address the question of the π 0 -η-η mixing in an unusually large context. Moreover, as seen in Subsect. 6.3, the exact structure of the ωππ coupling discussed several times in the literature [48,49,64] can also be examined within the largest possible data set.
A last remark is worth being emphasized. The scale treatment and the partial width ratio quoted in (i) and (ii), within the upgraded breaking of the HLS model show up as two different aspects of the same mechanism. Indeed, the former proceeds from applying the extended BKY-BOC breaking scheme to L V , while the latter follows from applying the same mechanism to L A .

The data sets and their handling
In this Section, we outline the data sets submitted to the global fit and the way correlated and uncorrelated uncertainties are dealt with. Nothing really new is involved here compared to what is already stated in [24,39,46], except for the data sets associated with the e + e − → K + K − and e + e − → K 0 K 0 cross sections. One may, thus, consider that this Section is, to a large extent, a simple reminder provided in order to ease the reading of the present paper.

9.1
The e + e − → π + π + data Four data sets have been collected recently in Novosibirsk at the VEPP2M ring. The first one [74,75], covering the region from about 600 to 960 MeV, is claimed to carry a remarquably small systematic error (0.6%). Later, CMD-2 has published two additional data sets, one [76]-covering the energy region from 600 to 970 MeV-is supposed to reach a systematic error of 0.8%, and a second set [77] closer to the threshold region (from 370 to 520 MeV) has an estimated systematic error of 0.7%. On the other hand, the SND collaboration has published [78] a data set covering the invariant mass region from 370 to 970 MeV. Except for the two data points closest to threshold which carry a sizable systematic error (3.2%), a reported systematic uncertainty of 1.3% affects this spectrum. These four data sets may be referred to in the following as "new timelike data" [39]. When dealing with these data sets, statistical and uncorrelated systematic uncertainties have been added in quadrature as usual. However, these four data sets also carry a common correlated systematic uncertainty estimated to 0.4% which affects all of them in the same way [79]. This is accounted for by modifying appropriately the covariance matrix as outlined in [39,46]-see also Subsect. 9.7 belowand by accounting for the data set to data set correlations. This is performed by treating these four data sets altogether, as if they were subsets of a single (merged) data set.
In order to be complete, we have also included in our fit all data on the pion form factor collected formerly by the OLYA and CMD Collaborations as tabulated in [80] and the DM1 data [81] collected at ACO (Orsay). These data will be referred to globally as "old timelike data". The systematic uncertainties carried by OLYA data (4%) and CMD (2%) contain an uncorrelated part which has been added in quadrature to the reported statistical errors. A common correlated part of the systematics, conservatively estimated [79] to 1%, has been dealt with appropriately. Instead, the accuracy of the DM1 data set being poor and its weight marginal, we did not find any need to go beyond the published uncorrelated errors. 9.2 The e + e − → (π 0 /η)γ data Since 1999, several data sets on the anomalous annihilation channels e + e − → π 0 γ and e + e − → ηγ have been made available by the CMD-2 and SND Collaborations. In our analysis, we only use the provided data points up to √ s = 1.05 GeV. The first one used is the data set from CMD-2 [82] on the ηγ final state (η → π + π − π 0 ) which carries a systematic error of 4.8%. CMD-2 has also provided [83] a second data set on the ηγ final state, tagged with the decay mode η → 3π 0 . The systematic uncertainty carried by this sample is estimated to 6.1% and 4.1% for, respectively, the energy regions below and above 950 MeV. More recently, CMD-2 has also published two more data sets [84] covering both the (π 0 /η)γ final states, tagged with the 2-photon decay modes, in the energy region from 600 to 1380 MeV. These are reported to carry a 6% systematic error.
The SND Collaboration has recently published [85] two different data sets for the ηγ final state with an estimated systematic uncertainty of 4.8%. The first one covers the energy region from 600 to 1360 MeV and the second from 755 to 1055 MeV. A sample covering the energy range from 600 to 970 MeV for the π 0 → γ γ decay mode was also published [86]. Finally, two data sets for both (π 0 /η)γ final states with 14 data points (from 985 MeV to 1039 MeV) from SND [87] are also available; these exhibit the much lower systematic error of 2.5%.
Altogether, these two Collaborations have provided 86 measurement points for the e + e − → π 0 γ cross section and 182 for e + e − → ηγ for √ s ≤ 1.05 GeV. Preliminary analyses [46] did not reveal any need to split up correlated and uncorrelated parts of the systematic errors for the (η/π 0 )γ data samples. Nevertheless, we have made a few checks by comparing fit results derived by adding in quadrature statistical and systematic uncertainties with fit results derived assuming the reported systematic error to be 100% bin-tobin correlated. We did not observe any significant difference. Therefore, when analyzing the e + e − → (π 0 /η)γ data, the reported statistical and systematic uncertainties have been simply added in quadrature as in [46].

9.3
The e + e − → π 0 π + π − data This channel is important as it provides a single place where the box anomaly sector [70,71] is present. Other physics channels involving the box anomaly in the η/η sectors exist (η/η → π + π − γ ) and may be relevant. However, the overall experimental situation is unclear [37,46], even if the Crystal Barrel data sample [53] may look secure. Therefore, we find preferable to wait for confirmation with new data samples which could come from BES and KLOE.
There are several published data sets for the e + e − → π 0 π + π − annihilation channel with various statistical and systematic uncertainties. We first included in our data sample the data sets collected by CMD-2 which consist of a measured sample covering the ω region [75] affected with a global scale uncertainty of 1.3% and two others which cover the φ region with a reported scale error of, respectively, 4.6% [88] and 1.9% [89]. The most recent CMD-2 data sample [90] also covers the φ region with a scale uncertainty of 2.5%.
SND has published two spectra covering altogether the region from 0.44 to 1.38 GeV, the former below 980 MeV [91], the latter above [92]. For both data samples, the correlated part of the systematic uncertainty has been extracted in order to be treated as a scale uncertainty (3.4% for [91] and 5% for [92], respectively); the uncorrelated parts have been added in quadrature with the reported statistical errors.
Former data sets are also considered which cover the region in between the ω and φ peaks where physics constraints are valuable. The most useful has been collected by the ND Collaboration with 10% systematics and can be found in [93], the latter is a small data sample from CMD [94] providing 5 measurement points with 15% systematics in the intermediate region. Concerning these two complementary data samples, we perform as in [46] and do not extract the correlated part of the systematics as the accuracy is poor enough that this could not lead to visible effects in global fits. Finally, there also exists a small data sample from DM1 [95] which has been used for illustrative purposes only [46].
The analysis of these data samples has been performed in [46]; however, as the N 5 term which contributes to the cross section (see (70)) was missing, the analysis is redone and the conclusions revisited. 9.4 The τ ± → π ± π 0 ν τ data In the collection of data samples submitted to global fitting, we also use the ALEPH [40], CLEO [42] and BELLE [41] data sets. When dealing with τ data, it is important to note that the relevant quantity, sensitive to the spectrum lineshape and to its absolute normalization is given by: where Γ τ is the full τ width, B ππ the branching ratio to ππν, and 1/N dN (s)/ds is the normalized spectrum of yields as measured by the various experiments. The data published by the ALEPH Collaboration correspond directly to the quantity shown in the left-hand side of (72). Instead, each of CLEO and BELLE has published separately the normalized spectrum of yields and the measured branching ratio B ππ . In the τ data handling, we have considered the reported uncertainties on these measured B ππ 's as bin-to-bin correlated scale uncertainties; these come into the various χ 2 associated with each data set in the way reminded in Subsect. 9.7. Stated otherwise, they are no longer fitted as previously done [24].
Following closely the experimental information provided by [40][41][42], the scale uncertainties have been estimated to 0.51% (ALEPH), 1.53% (Belle) and 1.74% (CLEO). On the other hand, a possible absolute energy scale uncertainty of 0.9% r.m.s. affecting the CLEO data sample [42] has not been found significant [24,39] and is not considered in the present study. All these experiments have provided their statistical and systematic error covariance matrices; these are the main ingredient of the χ 2 functions used in the fits.
As the HLS model relies on the lowest mass vector meson nonet only, it cannot access Γ τ which is therefore taken from the Review of Particle Properties [96]. Finally, our model provides [39]: with: and F τ π (s) is given in (37). Isospin symmetry breaking specific of the τ decay will be considered and taken into account as emphasized in Sect. 12.
Of course, the published τ spectra extend much beyond the validity range of the HLS model, as this presently stands. Therefore, when using it, we have to truncate at some s value. Consistency with the treatment of scan data would imply a truncation at 1.05 GeV. However, various studies [24,41] showing the behavior of fit residual clearly observe that ALEPH data on the one hand and Belle and CLEO data, on the other hand, exhibit inconsistent behavior starting in the 0.9÷1. GeV region. Therefore, we have preferred truncating the spectrum at 1. GeV, where the three spectra are in reasonable agreement with each other. 9.5 The e + e − → KK data Several data sets have been collected by the CMD-2 and SND Collaborations on both annihilation cross sections e + e − → K + K − and e + e − → K 0 K 0 . Here also, we have discarded the data points above 1.05 GeV. The oldest data sets, published by CMD-2 [88], provide the spectra for both the neutral and charged decay channels with a systematic uncertainty of 4%. Recently CMD-2 has reanalyzed four data sets for the neutral decay mode [97] getting small systematic errors (1.7%). More recently, CMD-2 has also published two scans of the charged mode spectrum [98] with a systematic uncertainty of 2.2%.
On the other hand, SND has published in 2001 several data sets [68]: 2 for the charged decay channel with a systematic error of 7.1%, 2 data sets in the neutral mode with K S → π 0 π 0 and 2 more with K S → π + π − , with respectively 4.2% and 4.0% systematics.
The quoted systematics are treated as correlated scale uncertainty as outlined in Subsect. 9.7 below.

The partial width data set
In order to work out the fit procedure and get enough constraints on the physics parameters of the model, an important input is the set of decay partial widths [39]. All decay modes of the form V P γ and P γ γ not related with the cross sections listed above should be considered. This covers the radiative partial widths ρ ± → π ± γ , η → ωγ and φ → η γ on the one hand and (η /η/π 0 ) → γ γ on the other hand. They have been extracted from the Review of Particle Properties [96]. The accepted values for radiative partial widths for K * ± → K ± γ and K * 0 → K 0 γ have also to be used [96].
As the currently available data on e + e − → π + π − stop slightly below 1 GeV, the phase of the φ → π + π − amplitude and its branching ratio as measured by SND [99] are relevant pieces of information, not included in the above listed annihilation data. 29 In contrast, the corresponding information for the ω meson is irrelevant as it is fully contained in the amplitude for e + e − → π + π − (see (41)) and is already part of the data sample.
With respect to former studies within the same framework, the only new piece of information included in the fit data set is the partial width π 0 → γ γ . Indeed, as can be seen from (58), the corresponding amplitude may constrain A as well as the e + e − → KK annihilation amplitudes. In fits involving all the above quoted annihilation channels, one has no longer to consider the leptonic widths (ρ 0 /ω/φ) → e + e − and the decay widths (ρ 0 /ω/φ) → (η/π 0 )γ as they are essentially extracted from some of the cross sections listed above which permanently enter our fit procedure.
Therefore, the additional decay information to be used as input to final fits represents in total 10 more pieces of information.
9.7 Outline of the fit procedure (the method) For all data sets listed above, one always has at one's disposal the statistical error covariance matrix. For scan data, this may include the uncorrelated part of the systematic errors; if not done at start, enough information is generally provided to allow one to perform this (quadratic) sum. In the case of τ data, the systematic error covariance matrix may 29 However, one might have to be cautious with these data. Indeed, as emphasized in [46]-see Sect. 13 therein-the single piece of information truely model independent is the product B ee B ππ . Therefore separate values for B ee and B ππ , given as "experimental" values in the various releases of the Review of Particle Properties, are actually model dependent to an unknown extent. be provided by the experimental groups (as ALEPH [40], for instance).
In this case, for each group of data sets (π + π − , π 0 γ , ηγ , π + π − π 0 , K + K − , K 0 K 0 , π ± π 0 ν) one computes the partial χ 2 : using matrix notations, and denoting by m and M the measurement vector and the corresponding model function vector. V is the error covariance matrix already referred to. The function to minimize is simply the sum of the χ 2 i . Actually, this is the procedure to estimate χ 2 i when the corresponding data sample is not subject to an overall scale uncertainty. If such a scale uncertainty takes place for some data set, one should perform a modification.
Let us assume that the data set i is subject to a scale uncertainty; this is supposed 30 to be a random variable ε(0, σ ) of zero mean (unbiased) and with r.m.s. σ , independent of s. Then any fit corresponds to getting one sampling of ε(0, σ ), named λ i . In this case, (75) should be modified to: where [100] A is traditionally the vector of the model values M and the other notations are obvious. One can solve for λ, which turns out to perform the change: in (75). The modified covariance matrix W depends on the vector A. As just stated, the best motivated choice for the vector A is the model function A = M. However, this implies a recursive determination of the modified covariance matrix, and, therefore, recalculating (or inverting) large matrices at each step of the minimization procedure (several hundreds of times for each fit attempt). It happens, however, probably because the experimental data we deal with are already accurate enough, that choosing A = m (i.e. the measurement vector of the corresponding experiment) does not sensitively affect the results and strongly improves the convergence speed of the minimization procedure [46]. Therefore, unless otherwise stated, we always perform this approximation. 30 In practical use, a data set # i, subject to a scale uncertainty λ i,0 is supposed to have been corrected in order to absorb a possible bias; this is the reason why the corresponding random variable is supposed unbiased, e.g. carrying zero mean. If not, (76) should be modified by performing λ i → λ i − λ i,0 .

The discarded data sets
There exists data sets which have been discarded for the present study. The most important are the three data sets collected using the Initial State Radiation (ISR) method by the KLOE [17,19] and BaBar [18] Collaborations. These suppose a specific statistical treatment as the structure of the reported systematic errors is much more complex than for any set of scan data. The method used in [46] for KLOE 2008 data [17] allows to deal with, but should be studied carefully with each ISR data set separately. In order to keep clear the message of the present study, we prefer avoiding using now data sets invoking delicate statistical methods. Therefore, the ISR data sets [17-19] will be treated in a forthcoming publication. Because of their high statistics, if well understood, these data samples may improve the physics results derived by using the model and the fit procedure presented in this study.
Other data sets could have been useful: • Those providing the pion form factor in the spacelike region close to s = 0 [101,102]. Indeed such data could severely constrain the pion form factor in the threshold region. This was illustrated in [39] where an archaic form of our model has been used. However, we gave up using them-especially [101]-because there is some suspicion concerning their estimated overall scale. Such a kind of data would nevertheless help in getting more precise information on g − 2.
• More data involving the box anomaly, especially in the η/η sectors may also help in constraining the model parameters. For instance, the dipion spectra in η/η → π + π − γ provide such information. Some available data collected in [37], especially those for η → π + π − γ provided by the Crystal Barrel Collaboration [53], might be considered sometime. However, new data sets on this subject, with larger statistics and better systematics should come from the KLOE and BES Collaborations, especially concerning the decay process η → π + π − γ . These are certainly more easy to handle than the e + e − → ηπ + π − annihilation data which in fine carry the same physics information.

The physics parameter set
It looks appropriate to give the list of the free model parameters to be fitted from data. The model parameters are of various kinds: • The basic HLS (4) parameters: the universal vector coupling g; the relative weight a of the Lagrangian pieces L A and L V , expected a 2 from most VMD models; finally the weights c 3 , c 4 and c 1 − c 2 of the anomalous FKTUY Lagrangian pieces to be added to the HLS Lagrangian in order to address the full set of data outlined in the above Subsections.
• SU(3) breaking parameters which modifies the physics content of the HLS Lagrangian (z A , z V and z T ), together with the parameter named λ which accounts for nonet symmetry breaking in the pseudoscalar sector. This amounts to a total of 4. • The isospin breaking parameters A , Σ V , V and h V which affect the non-anomalous HLS Lagrangian. These represent the Direct Isospin Breaking mechanism introduced in this paper through the BKY mechanism. • Some parameters [61] allowing the π 0 -η-η mixing. The η-η mixing angle θ P and the parameters named above and , which may account for, respectively, the π 0 -η and π 0 -η mixings. The last couple of parameters is not important for g − 2 estimates but may provide interesting physics information. One may anticipate on fit results by saying that the condition θ 0 = 0 is well accepted by the data as in previous analyses [38]; as a matter of consequence θ P can be (and will be) chosen as entirely fixed by the nonet symmetry breaking parameter λ (see (26)). One will also see that the pair and can be safely replaced by a single free parameter [61]. Therefore, the number of really free parameters accounting for the π 0 -η-η will be reduced to one. • Some subtraction parameters (8) involved in the mixing functions of vector mesons, in the ρ meson self-energy and in the γ -V transition amplitudes. • Some more parameters (4) describing the mass and width of the narrow ω and φ mesons. As a detailed description of the loop corrections to their inverse propagators is of little importance for the present purpose, there is no need to go beyond.
Stated otherwise, only the parameters A , Σ V , V and h V are new and all others have been already dealt with in previous releases of the present model [24,39,46].
One may be surprised to face a so large number ( 25) of parameters to be fitted from data. This only reflects that the number of physics pieces of information and of processes to account for is also exceptionnally large: more than 900 data points, six annihilation channels (π + π − , π 0 γ , ηγ , K + K − , K 0 K 0 , π + π − π 0 ), 10 radiative decay modes (V P γ , P γ γ including now the π 0 → γ γ partial width), the φ → π + π − decay and finally the dipion decay mode of the τ lepton.
All these pieces of information should get simultaneously a satisfactory description. It should be stressed that the parameter space is sharply constrained, as will be confirmed and illustrated by the present study. One should also stress that the π + π − , π 0 γ and ηγ cross sections, together with the decay modes referred to above, allow already a good determination of all fit parameters except for two: c 1 − c 2 and A . The former is derived from fitting the 3 pion cross section, the second from fitting both KK annihilation channels. Actually, in order to accurately determine Σ V , the dipion spectrum in the decay of the τ lepton also plays a crucial role.
This peculiarity leads us to a motivated critical analysis of the available π + π − π 0 , KK and τ data sets. As one plans to motivate a value for the hadronic contribution to g − 2, our dealing with the corresponding data should also be motivated.
As far as cross sections are concerned, it is already known from our previous studies that the π + π − , π 0 γ and ηγ annihilation cross sections are very well described within a simultaneous fit including also the decay data already listed. This can be seen in [39,46]; indeed Fig. 2 in [39] and Figs. 1 and 2 in [46] are indistinguishable from what is derived in the present study.
10 Reanalysis of the π + π − π 0 annihilation channel Taking into account the error described in Footnote 25, the analysis of the model description of the π + π − π 0 data is worth being redone. We take profit of this case in order to exemplify how the dealing with data sets is done.
The available 3-pion data sets can be gathered into 3 different groups: (i) The former data set collected by the Neutral Detector (ND) at Novosibirsk and published in [93]: we include in this group the few data points from [94]. These mostly cover the energy region in between the ω and φ peaks. (ii) A CMD-2 data set covering the ω region [75] together with a corresponding SND data sample [91] which actually extends up to 980 MeV. (iii) Several CMD-2 data sets covering the φ region and extracted from [88][89][90], accompanied by a data set from SND [92] starting at 970 MeV.
The small data sample from DM1 [95] is used for illustrative purposes and is not included in the fit procedures. It would not influence the fit results.
In fit procedures, it is very hard to run MINUIT normally because integrating the parameter dependent 3-pion cross section (see (70) and (71)) renders prohibitive the execution time. Therefore, we still use here the iterative method described and motivated in Sect. 10.3 of [46].
The choice of the 3-pion data sets considered in the global fit was performed in [46] relying on the data sets listed in (i). Indeed, the π + π − data used in the global fit serve to fix all parameters, except for the ω and φ mass and width parameters which are derived from having included the π 0 γ and ηγ cross sections; therefore, the ND data having a large lever arm (see downmost Fig. 5), they are alone able to determine accurately the value for c 1 − c 2 (see third line in Fig. 3).
Here one proceeds otherwise in order to learn more as each of the just above mentioned data set carries intrinsically a value for c 1 − c 2 . Nevertheless, the group of data sets needed in order to fix all parameters except for c 1 − c 2 has been enlarged: Beside the π + π − , π 0 γ and ηγ cross sections, we have included the τ decay information from ALEPH, Belle and CLEO. This will be justified later on. On the other hand, one assumes c 3 = c 4 which is justified in Sect. 15. Fits are performed by including either the CMD-2 data sets or SND data sets, each in isolation. On the other hand, separate (and independent) fits are performed in either of the ω and φ regions. Therefore, in these fits, the ω region fits are not influenced by the φ region information and conversely. Moreover, CMD-2 and SND data are not influencing each other. The data sets associated with the so-called ω and φ regions is not ours; it has been performed by the experimental groups who published the corresponding data sets separately.
It should be stressed, especially in the present case, that the notion of data set covers, as importantly, the data points, the full error covariance matrices (i.e. including the correlations reflected by the non-diagonal entries), and all the additional pieces of information provided by the experimental groups. Among this last kind of information, the global scale uncertainty included in the systematics should be suitably accounted for. As far as scan data are concerned, the statistical methods we use are the standard (text-book) methods briefly reminded in Subsect. 9.7.
The results of these fits are summarized in Fig. 1 and are commented on now. As a word of caution, it should be noted that the experimental errors shown in these plots are the quadratic sum of the reported statistical and systematic errors, neglecting all correlations. As the error bars do not (and cannot) take into account the correlations, they should only be considered as a visual indication of what is going on. The real distance of data points to its best fit curve is instead accurately reflected by the χ 2 values which, indeed, take appropriately into account all the reported pieces of information about the error covariance matrix.
Top left Fig. 1 shows the fit of only the CMD-2 data in the φ region; this provides a good fit 31 (χ 2 /npoints = 110/80) returning c 1 − c 2 = 1.21 ± 0.10. Top right Fig. 1 shows the case for the SND data in the φ region in isolation; the fit is much better (χ 2 /npoints = 26/33) but returns c 1 − c 2 = 2.18 ± 0.13. These two fit values for c 1 − c 2 differ by 10σ , clearly tagging an inconsistency between the CMD-2 and SND data sets in the φ region. 31 The numbers for χ 2 /npoints are the 3-pion sample contributions to the global χ 2 and its number of data points. One cannot provide the number of degrees of freedom as several hundreds of (other) data points are involved in each fit. Fig. 1 Best fits to e + e − → π + π − π 0 cross sections for data sets in isolation. Left column displays fits of the CMD-2 data, right column displays fits of the SND data. Top shows the φ region, bottom the ω region. The plotted data are extracted from [89,90] (CMD-2) and [92] (SND) for the φ region and from [75] (CMD-2) and [91] (SND) for the ω region. The empty circles (bottom right plot) are superimposed on the SND fit results and are not used in the fit displayed in this Figure On the other hand, one has performed likewise for the ω region in isolation. One then gets for CMD-2 data a large χ 2 /npoints = 26/13 with c 1 − c 2 = 1.29 ± 0.04 (bottom left Fig. 1). A closer examination of these data shows that an important part of this relatively large χ 2 is due to only the first data point which falls right on the vertical axis in this Figure. Instead, the SND ω region data yield χ 2 /npoints = 48/49 and c 1 −c 2 = 1.12±0.06 (bottom right Fig. 1). These two fit values for c 1 − c 2 differ by 3σ ; then, one may consider that the CMD-2 and SND data sets in the ω region are in reasonable agreement with each other.
One should note from fitting the SND ω data set, the important effect of correlations: In the bottom right Fig. 1, the large distance of the (SND) data points to their fitting curve is compensated by the correlations in such a way that χ 2 /npoints remains quite reasonable. The high level of compensation can be checked by computing the "diagonal" part 32 of the χ 2 which reflects the visual impression provided by the bottom right Fig. 1; one gets χ 2 diag = 554! 32 Denoting by V the full covariance matrix constructed as explained in Subsect. 9.7, the (full) χ 2 can be split up into its diagonal part χ 2 diag = i V −1 i,i ( i ) 2 and its non-diagonal part In addition, one has found instructive to plot the CMD-2 data together with the SND ones and the fit performed to the SND data in isolation. Thus, the bottom right Fig. 1 illustrates that the correlations reported by SND allow a reasonable reconstruction of the cross section valid for both the SND and CMD-2 data sets.
For information, the fit performed using only the ND data 33 yields χ 2 /npoints = 25/37 and c 1 − c 2 = 1.30 ± 0.06, in good accord with the previous fit result c 1 − c 2 = 1.17 ± 0.07, derived under comparable conditions (see second data column in Table 3 of [46]); the difference between these two estimates for c 1 − c 2 can be attributed to the influence of the τ data samples.
The various estimates for c 1 − c 2 derived from our fits are gathered in Fig. 3 using obvious notations. Using the fit values for c 1 − c 2 , as tag of consistency, this plot clearly shows that the φ region SND data set behaves differently from the other three-pion data sets.
From these considerations, one can conclude that: where i is the difference of the i th measurement and the corresponding value of the theoretical cross-section. 33 As reminded above, this data set covers the region in between the ω and φ peaks.

Fig. 2
Simultaneous fit of e + e − → π + π − π 0 cross section on the φ region data from [89,90] (CMD-2) and [92] (SND) • In the ω region, there is a good agreement between CMD-2 and SND data from within the filter of our model. • In the φ region, at minimum χ 2 , one can get a reasonable description of both CMD-2 and SND data, but with much different values for the fit parameters as reflected by their c 1 − c 2 values.
Therefore, one observes a qualitative difference between all CMD-2 data and the SND data in the ω region, on the one hand, and the SND data in the φ region, on the other hand. One has pushed a little further the analysis by two more series of fits: • One has simultaneously submitted to fit the CMD-2 and SND data but only in the φ region. One gets the result shown in Fig. 2. The fit might look reasonable (χ 2 /npoints = 176/113) and returns c 1 − c 2 = 1.94 ± 0.07, close to the SND value, as can be seen from Fig. 3. From this series of fit, one can conclude that it is possible to fit simultaneously the CMD-2 and SND data in the φ region and get a reasonable solution. However, mixing the ω and φ regions returns, in the case of SND, an unacceptable solution, even if the χ 2 /npoints may look reasonable.
Therefore, one is led to include in the set of data samples submitted to the global fit all 3-pion data referred to above, except for the SND φ region data set. The corresponding fit has been performed and is shown in Fig. 5 with c 1 − c 2 = 1.18 ± 0.03; the 3-pion data contribute to the global fit with χ 2 /npoints = 220/179. The result shown at the last line in Fig. 3 shows that the global fit performs, as expected, a good 34 Here also, one may wonder that the top right Figure corresponds to a quite reasonable fit quality. We thus remind once more that, for all figures shown, the effects of correlated uncertainties is not-cannot beshown. In the case of SND, this is larger than 5%. Along the same lines, one should mention that the errors plotted are always the quadratic sum of statistical and uncorrelated systematic uncertainties. 35 In this case, the so-called "diagonal" part of the χ 2 at minimum is larger than 1100. CMD2 ω denotes the fit result of the data from [75], SND ω those from [91], ND + CMD the fit result to the merged data from [93] and [94], CMD2 φ indicates that only the merged data from [88][89][90] have been used in the fit, SND φ corresponds to the fit of the data from [92] and SND + CMD2 φ provides the (simultaneous) fit result of [88][89][90]92]. Finally, the last line shows the result for the selected data consisting of the sample reported in [75, 88-91, 93, 94]. The vertical dotted line serves to show how the fits perform the averaging (fitted) average of c 1 − c 2 . This also indicates that the data sets considered are statistically consistent with each other.

Analysis of the KK annihilation data
As reminded in Subsect. 9.5 above, several data samples are available collected by the CMD-2 and SND Collaborations on VEPP-2M at Novosibirsk. The CMD-2 data are extracted from [88,97,98] and the corresponding SND data from [68]. The quoted systematics are treated as a scale uncertainty and dealt with as explained in Subsect. 9.7.
The published data being cross sections, the fitting function is: for each of the 2-kaon annihilation channels; the kaon form factors F e K (s) are given by (49). Both cross sections are corrected for the intermediate photon dressing. Moreover, for the charged kaon channel, the additional Coulomb factor [47,67] Z(s), reminded in (53), is understood and is not "renormalized" as in [68,69] with Z(m 2 φ ).

Fitting the KK data
In order to perform this analysis, we have done a first series of fits using separately the CMD-2 neutral and charged KK channels and the corresponding data from SND. In order to avoid φ peak information not following from the KK data, we have decided to remove the data from the π 0 γ and ηγ annihilation channels from the fit procedure. However, anticipating on our final results, we have included the three data sets from ALEPH, Belle and CLEO referred to above. Therefore, the additional data sample is composed of all e + e − → π + π + data, all τ ± → π ± π 0 ν data and 18 partial width decays (all V P γ and P γ γ modes and also the three leptonic decays (ρ/ω/φ) → e + e − modes). None among these pieces of information has any direct influence on the description of e + e − → KK, even through the φ mass and width parameters which are, thus, solely determined by the KK data.
The results are shown in Fig. 6, left side for the K 0 K 0 data and right side for K + K − . One observes a good description of the K 0 K 0 data for each of the samples provided by the SND or CMD-2 Collaborations. The picture is quite different for the K + K − data; the CMD-2 data sample is well fitted, while the SND sample is poorly fitted. Additional information for these peculiar fits is displayed in the first two Fig. 4 Simultaneous fit of the e + e − → π + π − π 0 data in the ω and φ regions. Top figures show the case for the merged data from [91,92]. Bottom figures display the fit results for CMD-2 data from [75,89,90] lines of Table 2. One can see there, that the value for χ 2 /N associated with the K 0 K 0 data are the same for both data samples, while they differ significantly for the corresponding K + K − data samples. Fitting simultaneously both CMD-2 and SND K 0 K 0 data samples only, returns the same χ 2 information, illustrating that the corresponding data samples are perfectly consistent with each other. Simultaneous fits of all KK data confirm this property (see third line in Table 2). Interestingly, the χ 2 's at best fit in the third and fourth lines practically coincide with the sum of the corresponding information in the first two lines of the same Table. This illustrates that the socalled additional data set sharply constrain the KK cross sections. Moreover, in view of the fit results for CMD-2 data, one can consider that the constraints are well fulfilled by data, giving a strong support to our modelling. The ratio of cross sections σ (e + e − → K 0 K 0 )/σ (e + e − → K + K − ) is observed to provide a valuable piece of information, as it allows to magnify the effects mentioned just above. This is shown in Fig. 7, where the data for this ratio are plotted normalized to the ratio of cross sections as coming out from our fits. The data ratio plotted in the top Fig. 7 Table 2 Fit quality of the K + K − and K 0 K 0 data. Beside the additional data sample (see text), each line in the first column tells which KK data samples have been included in the fit procedure. χ 2 0 is the χ 2 value for K 0 K 0 data, χ 2 c is the corresponding information for K + K − data. The N 's are the respective numbers of data points. The last data column provides the global fit probability for each case is derived from the information given in [69] and one can estimate its uncertainty to 2.3÷2.4%. The CMD-2 data points normalized to the fit expectations bin per bin is perfectly consistent with 1 over the whole s region covered by the φ resonance. The dotted lines in top Fig. 7 represent the experimental scale uncertainty and do not take into account the uncertainties on the fitting functions. This also illustrates that our modified Breit-Wigner lineshape is very well accepted by the data.
In contrast, the SND data exhibit a behavior reasonably well averaged by the fit function ratio; however, it does not look consistent with flatness-at least as well as for CMD-2 data.
It follows from these considerations that the largest selfconsistent data set for the KK channel is made by merging all CMD-2 data and the K 0 K 0 data provided by SND (see last line in Table 2). As a matter of information, beside getting an appropriate description of both e + e − → K 0 K 0 and e + e − → K + K − cross sections, it is worth remarking that the radiative partial widths included in the fitted data set are also well accounted for. For instance, including also the e + e − → (π 0 /η)γ cross sections in the fitted data set, the remaining set of 10 radiative decays yields a quite remarkable χ 2 /n = 6.5/10, with estimated Γ (π 0 → γ γ ), Γ (η → γ γ ) and Γ (η → γ γ ) at respectively 0.27σ , 1.77σ and 0.23σ from their accepted values [10]. As the corresponding couplings are strongly affected-especially g π 0 γ γ -by A (see (58)), we consider that physics validates our model.

The HLS solution of φ → KK puzzle
The partial width decays φ → KK are defined by: Therefore, one has: where R = 1.528 originates from the ratio of momenta and the Coulomb factor computed at the φ peak is Z(m 2 φ ) = 1.049. The ratio of couplings has been given in (54). Therefore, using A from the last line in Table 2, one gets: The same ratio can be computed from information given by CMD-2 in a recent paper [69] and amounts 36 to 1.47 ± 0.04. The difference between the CMD2 estimate and ours amounts to about 2σ . Our final result, obtained by using the largest possible ensemble of data sets, provides A = (6.34 ± 0.70) × 10 −2 and then the ratio of branching ratios becomes 1.40 ± 0.02.
Therefore, the HLS model, equipped with the (BKY) direct isospin symmetry breaking mechanism, provides a solution to the long-standing puzzle concerning the φ → KK decays as thoroughly analyzed in [47] and more recently 36 The uncertainty might be somewhat overestimated, as one has assumed independent the errors for Br(φ → K + K − ) and Br(φ → K 0 K 0 ). discussed in [68]. In our approach, the mechanism responsible for this is, in fine, the kaon field renormalization which should be performed within the HLS model once isospin symmetry breaking is performed à la BKY-BOC. Indeed, as the neutral and charged kaon field renormalization factors play in opposite directions (see (18)), they pile up in the ratio.
The relatively large value found for A indicates that several sources contributes to the BKY breaking of isospin symmetry. The contribution to A due to the light quark mass difference [58] ( 1%) is certainly not the single source and others-like electromagnetic corrections-are certainly absorbed within the numerical value for A . Moreover, it is also likely that different corrections at the V K + K − and V K 0 K 0 vertices may influence the fit value for A . Being global, the BKY breaking mechanism cannot allow to disentangle the various contributions to A which share a common order of magnitude (each at the percent level). The situation is quite different from the breaking of SU(3) symmetry which is widely dominant numerically and can motivatedly be compared to ChPT expectations [38].

Analysis of the τ decay data
Using F τ π (s), the pion form factor in the decay of the τ lepton (see Subsect. 6.2), the partial width of the two-pion decay is given by (73). On the other hand, the quantity which encompasses the full experimental information in this field is (72): dN(s) ds as, indeed, the lineshape and the absolute magnitude of each experimental spectrum are merged together. The full width Γ τ is taken from the RPP [10]. The last two factors in the middle expression above account for isospin symmetry breaking effects specific of the τ decay: S EW = 1.0235 for short range corrections [103], G EM (s) for long range corrections [104][105][106].
In former studies, it was shown that the lineshape alone was perfectly consistent with annihilation data [24,39]. However, if one also takes into account the absolute magnitude-represented by the branching ratio B ππ in the formula reminded just above-the agreement is poor. In order to reach a satisfactory description of the data, Ref. [24] introduced a mass difference δm 2 and a coupling difference δg between the neutral and charged ρ mesons, which underlays all reported stand-alone fits to τ spectra [13]. However, additional scale factors were needed and their fitted values were found consistent with the reported scale uncertainties [40][41][42].
However, the present study, as reflected by Table 1 above, has clearly demonstrated that isospin breaking of the HLS model does not necessarily result in non-vanishing δm 2 and δg at leading order. 37 As emphasized above, the BKY-BOC breaking scheme instead leads to a difference between the universal vector coupling (g) as it comes in the anomalous sector and in the non-anomalous sector of the HLS Lagrangian (g(1 + Σ V )). We prove, here, that this provides a much better account of all data than only assuming some mass and width differences supplemented with some residual rescaling. Stated otherwise, it is because Direct Isospin Breaking acts differently in the anomalous and non-anomalous sectors that the model yields an almost perfect description for all data, without any need for some additional rescaling. In this mechanism, the single sensible difference between the pion form factor in e + e − annihilations and in τ decays resides in the difference between the transition amplitudes γ -V and W -V . Fig. 8 Global Fit of the dipion spectrum in the decay of the τ lepton. The data points are those from ALEPH [40], Belle [41] and CLEO [42]. The inset magnifies the ρ peak region Figure 8 shows the global fit result for the function H (s) defined just above together with the data points from ALEPH [40], Belle [41] and CLEO [42] Collaborations. 38 The inset magnifies the ρ peak region. One can clearly conclude to a nice agreement between model and data, all along the fitted region-from threshold to 1 GeV. The corresponding pion form factor in e + e − annihilations coming out of the global fit is represented in Fig. 9. These two Figures illustrate that the simultaneous description of e + e − and τ data allowed by the model is, indeed, as successful in both sectors. Figure 10 shows in two different manners the τ residual behavior. Top Fig. 10 displays the usual residuals for the function H (s), while downmost Fig. 10 represents (H model (s) − H data (s))/H model (s). These can be compared with respectively Fig. 3 and Fig. 4 from [24] where the (δm 2 , δg) parametrization of isospin breaking was used. The comparison clearly indicates that the present model better performs for all τ data sets and, especially, for the ALEPH [40] data. 38 When dealing with τ plots, the error bars represent the diagonal errors, i.e. no account of bin-to-bin correlations is attempted.
In order to allow for a deeper comparison with the previous release [24] of the present model, we reproduce in Table 3 (first data column) the fit results reported in [24] together with our new fit results under various conditions.
The second data column in Table 3 is derived excluding the KK data sets in order to be as close as possible to [24]. One observes, for almost all data sets, better fit results than in the former release of our model [24]. There is no effect in introducing the 3-pion data set from SND [91] (covering the ω region) as the χ 2 3π /dof = 1.11 is unchanged. It is also worth noting that the partial width for η → γ γ is found at 0.43σ from its accepted value [10]; the distance is 0.11σ for η → γ γ and 0.47σ only for the newly introduced π 0 → γ γ decay mode.
One may conclude therefrom that the HLS model, equipped with the mixing schemes provided by loops and by the direct isospin breaking procedure, provides a fully satisfactory solution to the e + e − − τ puzzle, both in magnitude and in shape. The relatively poorer fit quality for the BELLE data might be related with the absolute scale issue Fig. 9 Global Fit of the pion form factor squared in e + e − annihilations. The data points are those from CMD-2 [75][76][77] and SND [78]. One has not plotted the so-called "old timelike" data also (mostly) collected at Novosibirsk. The inset magnifies the ρ peak region and the behavior at the ρ-ω interference region Table 3 Comparison of the fit qualities between the fit results of the model as it was in [24] (second data column) and as it is now (third data column). KK data were not submitted to fit in [24]. The '+1' added to the number of data points for τ data stands for the experimentally given r.m.s. affecting the (fitted) global scale. The 3-pion data set information is displayed boldface in order to show the difference in the fit data set: In the second data column, the 3-pion data set from SND [91] has been (newly) introduced and in the last data column only the 3-pion data sets collected below the φ region are considered  The fitted region extends from threshold to 1.0 GeV/c, i.e. over the region where the behavior of the data sets from ALEPH [40], Belle [41]and CLEO [42] reach some agreement Fig. 11 Ratio of the transition amplitudes ρ 0 -γ and ρ ± -W ± , f ργ /f ρW following from the global fit and neglecting loop corrections. This corresponds to the ratio shown in Table 1 and reproduced in Sect. 12. Top figure shows the real part as a function of s, bottom figure the imaginary part. Uncertainties due to fit parameter errors are not given; the uncertainty band for f ργ /f ρW − 1 can estimated to a few percent Fig. 12 Ratio of the couplings g ωππ /g ρππ as a function of √ s, as coming from the global fit (this ratio is explicitly given in Sect. 13). The vertical line locates the PDG mass of the ω meson. The uncertainty band due to fit parameter errors is not shown Fig. 13 A set of recent estimates of the muon anomalous magnetic moment a μ together with the BNL average value [1,2]. These are extracted from [14] (DHMZ10), [16] (JS11), [115] (HLMNT11) and [13] (DHea09). Our own results are figured by A and B for respectively solutions A and B. The statistical significance of the difference between the estimated and measured values of a μ is displayed on the right side of the Figure for each of the reported analyses revealed by the stand-alone fit 39 provided by BELLE [41]. Therefore, one can confirm that: • The main drawback of the breaking model in [24] was a too tight correlation between the universal coupling 39 The fit published by BELLE reveals a very significant improvement if the absolute normalization of their spectrum is left free; instead of returning an absolute scale of 1, the best fit exhibits a significant 2% shift.
in anomalous and in non-anomalous processes. This has been cured by defining the Direct Isospin Breaking mechanism substantiated by a highly significant value for Σ V = (3.74 ± 0.42)%. • The breaking model in [13] may account insufficiently for the difference between the ρ 0 -γ and ρ ± -W ± transition amplitudes.
Therefore, the reported discrepancies between the pion form factor in e + e − annihilations and in τ decays can al-ways be attributed to an incomplete treatment of isospin symmetry breaking. For information, Fig. 11 displays the ratio of the transition amplitudes f ργ and f ρW as coming from the global fit and already given in Table 1: We have found appropriate to provide in the third data column of Table 3 the results of the fit obtained keeping the KK data sets, while excluding all the π + π − π 0 data sets. The fourth data column reports on the fit quality reached using the full data set we considered safe. This means all data sets discussed above, except for two SND data sets: The e + e − → 3π data set collected above 970 MeV [92] and the e + e − → K + K − data set. These have been shown to provide either an unacceptable behavior for the fit solution [92], or a poor χ 2 [68]. In this configuration, one fits 906 data points (including the 10 individual decay modes) corresponding to 881 degrees of freedom. The global fit probability is highly favorable (71%). This configuration will be referred to in the following as "Solution A" or "Configuration A".
In this Solution A, one observes some tension between the KK and π + π − π 0 data groups. Indeed, comparing its content with the second data column, one observes that the π + π − π 0 data group yields a χ 2 increased by 30 units. Instead, comparing Solution A with the third data column in Table 3, one does not observe any significant degradation the fit quality of the KK data group: The χ 2 for the K 0 K 0 data group is improved by 2 units, while the χ 2 for K + K − data group is worsened by 6 units.
As this 30 unit increase of the (π + π − π 0 ) χ 2 may look abnormal, we have tried tracking its origin. This issue is clearly related with having introduced the KK data which influence the model description of the φ region. Therefore, we have redone fits excluding all the π + π − π 0 data sets covering the φ region. One obviously remarks a significant effect; this configuration will be named hereafter "Solution B" or "Configuration B".
In the following, any differential effect between what has been named Solutions A and B is examined carefully and commented.

Structure of the ω → ππ coupling
As noted in Subsect. 6.3, the coupling ω → ππ in the upgraded broken HLS model is given by: This expression exhibits two contributions of different origin. The first part is a constant term generated by the Direct Isospin Breaking procedure defined at the beginning of this paper, the second is generated by the kaon loop mixing procedure already defined in [39,46] and reminded above. This structure resembles that given in [48,49]. It is interesting to examine the behavior of the ratio: as a function of √ s. It is given in Fig. 12, where the vertical line figures the ω mass location. Of course, the effective part of this function is determined by the ω Breit-Wigner distribution and is concentrated within a few tens of MeV's apart from the ω peak position.
From the best fit discussed in the above Section (see the second data column in Table 3), one gets the central values for the fit parameters and their error covariance matrix. These have been used to generate g ωππ by Monte Carlo methods. Computed with using the RPP [10] mass for the ω meson, this gives 40 : The observed useful correlations are δΣ V δ V = −0.056, δΣ V δh V = 0.028 and δh V δ V = 0.232.
In order to stay consistent with [48,49] definitions, one can consider that g I ρππ = ag(1+Σ V )/2 and g I ωππ = ag(1− h V ) V /2 are the couplings of the ideal fields, defined as such before applying the loop mixing. Therefore, the quantity G: should be close to the parameter carrying the same name in [49]. One finds G = (3.47 ± 0.64) × 10 −2 to be compared with the two estimates of the same parameter given in [49]: G = (7.3 ± 3.2) × 10 −2 when relying on the data from [75] and G = (4.4 ± 0.4) × 10 −2 when using, instead, the [76,77] data. Referring to (28), one can conclude 41 that there is much more isospin 0 inside the physical ρ than isospin 1 inside the physical ω. In this case, one also gets for the direct term ag(1 − h V ) V /2 = −0.332 ± 0.024. Comparing this number with Re(g ωππ ), it is clear that ag(1 − h V ) V /2 and 40 The quoted uncertainties for V , h V , a and Σ V are the improved uncertainties returned by the routine MINOS of the MINUIT package [107]. 41 The isospin 0 component inside the physical ρ meson is given by h V V , while the isospin 1 part inside the ω is given by Re(α(m 2 ω )) compensate to a large extent, in such a way that g ωππ is highly dominated by its imaginary part. 42

The π 0 -η-η mixing properties
The mixing of pseudoscalar neutral mesons has been addressed in Sect. 4, especially in Subsects. 4.3 and 4.4. The present Section is devoted to examining how the upgraded breaking scheme developed in this paper performs compared to the results previously derived in this field [38]. In order to perform this study, we let free the pseudoscalar mixing angle θ P , which mostly determines the relationship between the physical η and η fields and their underlying octet and singlet components η 0 R and η 8 R . The parameters and which account for, respectively, the π 0 -η and π 0 -η mixing are also let free.
As shown in [38] and revisited in Subsect. 4.3 above, the ChPT mixing angles [50,51] θ 0 and θ 8 can be expressed in terms of the nonet symmetry breaking parameter λ (or, better, using instead v defined in (22)), z A the SU(3) breaking parameter of the Lagrangian L A and the singlet-octet mixing angle θ P . Therefore, they can be estimated from fitting the data already defined.

and θ P
The mixing angles θ 0 and θ 8 have been recently introduced with the 2-angle description of the η/η mixing [50,51]. The broken HLS model provides expressions for these in terms of the singlet-octet mixing angle θ P and of the breaking parameters z A and λ (see (26) and also [38]).
Therefore, using the fit results (parameter central values and their error covariance matrix) one can reconstruct the values for θ 0 and θ 8 . Having left free θ P , one obtains the results shown in the first data column of Table 4. Therefore, as in former studies, one observes that θ 0 is small and its distance to zero is only 2.8σ ; this should be compared with the estimate θ 0 = −4 • given with no quoted uncertainty in Table 4 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)

General Fit
Constrained Fit  [51]. The value for θ 8 is numerically as expected from other kinds of data [51]. The 'tHooft parameter [35] λ is found of the order 10 %, twice smaller than in [38] where an approximate treatment of nonet symmetry breaking was used. Finally, the singlet-octet mixing angle θ P is still found twice smaller than θ 8 , as in the former study [38].
As the distance to zero of θ 0 is 2.8σ , the non-identically vanishing of θ 0 is on the border of statistical significance. Therefore, imposing the condition θ 0 = 0 is worth being considered; this turns out to algebraically relate θ P to z A and λ by tan θ P = tan B (see (26)). Performing such a fit returns the results shown in the second data column of Table 4 with a quite comparable probability.
It is interesting to observe that the value for θ 8 is nearly unchanged and that the value for λ is affected below the 10 −4 level only. One also observes that the value for θ P generated by the appropriate (26) is found in agreement with its fitted value (when this parameter is left free). We conclude therefrom that assuming θ 0 = 0 does not degrade the fit quality and is consistent with data.
One should also note that the nonet symmetry breaking parameter λ = 8.5% has a statistical significance of 2.4σ . Performing an approximate nonet symmetry breaking [38], the value for λ was overestimated by a factor of 2.
14.2 The π 0 -η and π 0 -η mixing properties These mixing properties are reflected by the parameters named respectively and as displayed in (23). Comparing analogous fits performed by letting free and unconstrained θ P , and , we did not find sensitively different results than those obtained by imposing the constraint on θ P resulting from the condition θ 0 = 0. Therefore, from now on, all presented fit results will refer to this configuration. One should note that the numerical results given in the above Sections have also been derived under these conditions. The global fit returns = (4.89 ± 0.44) × 10 −2 and = (1.68 ± 0.44) × 10 −2 , reflecting that the π 0 -η mixing is certainly much more important than the π 0 -η mixing phenomenon. With the concern of reducing the number of free parameters, we have also assumed [61]: 2 cos θ P −sin θ P √ 2 cos θ P +sin θ P = −2 0 sin θ P √ 2 cos θ P +sin θ P √ 2 cos θ P −sin θ P (85) with θ P still determined by the constraint θ 0 = 0. This reduces the number of free parameters by one more unit. The fit returns 0 = (3.16 ± 0.23) × 10 −2 with an unchanged probability; this corresponds to values for and very close (2σ each) from the corresponding fitted values, while the global fit probability is unchanged. The partial widths for the three decays P → γ γ are all well accounted for: 1.64σ (η), 0.11σ (η ) and 0.06σ (π 0 ). Additional fit detail can be found in Table 3.
The question of whether the present 0 can be identified with the variable carrying the same name in [61] is unclear. 43 Indeed, an important part of isospin symmetry breaking effects are already included in the definition of the renormalized PS fields (see (19) and (21)) which undergo the rotation defined by (23). Therefore, our , and 0 carry only a part of the isospin breaking effects, while another part (governed by A ) has been propagated to all sectors of the effective Lagrangian.

The values of the FKTUY parameters
Our global fit modelling is in position to provide the most accurate information concerning the parameters c 3 , c 4 and c 1 − c 2 defining the scales of the various FKTUY anomalous pieces [28] of the HLS Lagrangian.
From our results, which happen to be the most precise in this field, one may conclude that data only favor a partial fulfilling of the VMD assumptions [23], in the sense that c 3 − c 4 = 0 is in agreement with data at the 2σ level, while c 1 − c 2 + c 4 = 4/3 is badly violated. This can be rephrased as follows: the VMD assumptions [23] are experimentally fulfilled in the triangle anomaly sector and strongly violated in the box anomaly sector. This confirms the previous parent analysis [46] and former studies on the box anomaly in the η/η → π + π − γ decays [37].
In order to go beyond, better data on the annihilation channels involving anomalous couplings ([π 0 /η]γ , π + π − π 0 ) are needed; including new processes like the η/η → π + π − γ decay spectra or information on the l + l − π 0 annihilation channels may also help as their dependence upon c 3 − c 4 or c 1 − c 2 is more important than in the previous channels.
It thus follows from the present analysis that assuming c 3 = c 4 is justified. In this case, one obtains the following results: for Configuration A and: c 4 = c 3 = 0.951 ± 0.016, c 1 − c 2 = 1.169 ± 0.060, g = 5.553 ± 0.012 (89) for Configuration B. In both cases, the correlation coefficient is −0.20. Therefore, the condition c 4 = c 3 drastically reduces the correlation among the surviving FKTUY parameters. Moreover, the fit quality is not significantly changed while assuming c 4 = c 3 . Indeed, configuration A yields χ 2 /dof = 858.08/882 (71.2% probability) instead of χ 2 /dof = 853.98/881 (73.7% probability) and configuration B χ 2 /dof = 728.38/882 (97.0% probability) instead of χ 2 /dof = 722.05/881 (97.9% probability), where the difference mostly affects the set of partial widths which is always well fitted. Therefore, the improvement obtained with the upgraded breaking model is not due to releasing the condition c 4 = c 3 . From these considerings, it is justified to impose c 3 = c 4 for the rest of this study.

Hadronic contributions to g − 2
In [24], one analyzed in full detail the hadronic contribution to g − 2 of most of the data sets used in the present study. The framework was the previous release of the present model studied in detail in [24,46]. Within this framework, only the simultaneous account of both annihilation channels to KK was missing. On the other hand, one might find unsatisfactory that some global rescaling of experimental τ dipion spectra was still playing an important role, even if this rescaling was in accord with expectations. These two issues motivated the present study.
As shown above, the upgraded model allows by itself a satisfactory account of all considered spectra simultaneously. It is therefore worth reexamining within our upgraded framework, how the hadronic contribution to g − 2 is estimated and how this estimate evolves depending on the various kinds of data groups considered.
16.1 The π + π − contribution to g − 2: VMD estimates The most important hadronic contribution to g − 2 is the π + π − channel. Several experiments [75,76,78] and some analyses [13,24,110] give the π + π − contribution to g − 2 provided by the energy region [0.630-0.958] GeV. Therefore, it is worth considering the information provided by this reference region; this allows to substantiate the improvement which can be expected from VMD-like models. Indeed, several kinds of information are worth considering: • While unifying the description of e + e − annihilation and τ decays, one expects an increased precision on the anomalous magnetic moment of the muon a μ (ππ). • While having a framework which encompasses most of the physics up to the φ region, the stability and the robustness of the a μ (ππ) estimates can be examined. The relative statistical consistency of the various data groups is also an issue which can be addressed, relying on their behavior under global fits. Table 5 displays our estimate for the ππ contribution to a μ = (g − 2)/2 provided by the reference energy range under various fit configurations. In each case, the fitted (central) parameter values and their error covariance matrix are used in order to sample several thousand parameter vectors, assuming a n-dimensional Gaussian error distribution. Each vector of sampled parameter values is, then, used to compute a μ (ππ). The corresponding distribution of the a μ (ππ)'s is then fitted to a Gaussian function. The results displayed in Table 5 are the central values and the standard deviations of this distribution which intrinsically takes into account the correlations among the fitted parameters. Unless otherwise stated, the FSR correction is included in all reported contributions of the π + π − channel to a μ .
Beside the experimental spectra, there is always a set of partial width decays submitted to fit. These have been defined in Subsect. 9.6. In the results reported below, one should keep in mind that the accepted values [10] for the (ρ/ω/φ) → (π 0 /η)γ and (ω/φ) → e + e − partial widths are included in the set of partial widths submitted to the fit as long as the experimental spectra for the e + e − → (π 0 /η)γ annihilation channels are not used. As emphasized in [46], this hides some model dependence which might be somewhat conflicting with our own model. This explains why one should prefer any configuration where the e + e − → (π 0 /η)γ data are submitted to the global fit.
Also, when the data for the two annihilation channels e + e − → KK are not considered in the fit, one chooses to fix 0 = A = 0, as we have no real sensitivity to them. Likewise, c 1 − c 2 is absent from fits as long as the e + e − → π + π − π 0 data are not considered. Finally, the parameters fixing the mass and width of the φ meson are left free only when the fitted data allow to constrain them.
In the first line of Table 5, one finds the value for a μ (ππ) derived by submitting to fit the scanned data for the annihilation process e + e − → π + π − -together with the full set of partial width decays. This result compares well with the value derived using the previous release of our broken Table 5 The contribution to 10 10 a μ (ππ) from the invariant mass region 0.630-0.958 GeV/c. The first line provides the fit results using all the e + e − → π + π − annihilation data set group. The next line uses the previous data group and the three τ spectra. By "++" at any given line, we always mean all data sets belonging to the groups referred to in the preceding lines, plus the data set group indicated at this line. FSR corrections are taken into account. An appropriate set of radiative decays is always understood. HLS model, 45 as can be seen by comparing with the relevant piece of information reported in Table 4 of [24]. As there is no longer any mismatch between e + e − and τ data, both in magnitude and in lineshape (see Sect. 12), it is legitimate to merge them. This merging provides the new and important result given in the second line of Table 5. One clearly observes that the merged ππ data give a result perfectly consistent with the e + e − → π + π − data alone with a quite nice probability. The central value for a μ (ππ) is nearly unchanged and the uncertainty slightly improved. This is, of course, the main effect of having upgraded our symmetry breaking procedure of the HLS Lagrangian. In this new framework, there is no need for an auxiliary rescaling [24] of the τ spectra and the net result is a perfect consistency of the e + e − → π + π − data with/without the τ data considered as constraints. This statement can be substantiated by comparing this result with those reported in the entry "NSK + A B C" of [24] (a μ (ππ) = (364.48 ± 1.34) × 10 −10 ) which exhibited a shift of about 5 × 10 −10 produced by the three τ data sets, a 3.6σ effect.
The third line in Table 5, displays the effect of replacing the (ρ/ω/φ) → (π 0 /η)γ and the (ω/φ) → e + e − partial widths by the cross sections for e + e − → (π 0 /η)γ . The central value for a μ (ππ) is practically unchanged, while its standard deviation is increased by 9%. The following line in Table 5 displays the effect of including the full e + e − → π + π − π 0 data group already defined. As in [24], one observes a perfect consistency of the results for a μ (ππ). In total, the standard deviation is slightly reduced (σ (a μ (ππ)) 1.5 × 10 −10 ). At this point one may conclude that the central value is marginally modified by fully including the (ρ/ω/φ) → (π 0 /η)γ and e + e − → π + π − π 0 data groups within the fit procedure. The variations of the uncertainty returned by the fits might rather reveal statistical fluctuations.
The last line in Table 5 displays the effect of including the two e + e − → KK cross sections into the fitted data set. One observes some effect, as a μ (ππ) undergoes a 1.9 × 10 −10 shift upwards while the fit probability remains quite good. This fit configuration-referred to as Solution/Configuration A-encompasses the largest set of data samples considered safe. This turns out to consider that the 30 unit increase of the χ 2 associated with the π + π − π 0 data group, even if large, is not abnormal (see the fourth data column in Table 3).
The result shown in the last line of Table 5, may reveal some tension among the data set groups. In order to explore this issue, one has redone fits excluding the π + π − π 0 data group, and examined the effects of using the selected 45 Even if expected, this proves that the effects produced by having introduced Σ V do not modify the fit description of the e + e − → π + π − data. e + e − → K 0 K 0 and e + e − → K + K − data, either separately or together. The corresponding results are displayed in Table 6. Comparing the statistical information here with those in the last line in Table 5 renders somewhat suspicious the quoted 30 unit increase of χ 2 π + π − π 0 . A final piece of information is provided by performing the fit using the π + π − π 0 data group data amputated from the data points collected in the region above 1 GeV (therefore, excluding the φ region). This fit configuration has already been referred to as Solution/Configuration B. The reason which motivates this removal is that the π + π − π 0 data before introducing the KK data is only constrained in the φ region by the relatively unprecise data on the π 0 γ and ηγ channels. One then obtains: where the result for Solution A is reminded.
These differences indicate that all physics channels covering the φ region are worth to be reconsidered, as already argued from discussing the fit results in Table 3. Indeed, the difference in fit quality between Configurations A and B reveals some tension between the KK data and the π + π − π 0 data collected in the φ region. Fortunately, the physics in the φ region is still accessible at VEPP-2M. It seems also in the realm of the KLOE detector, as this turns out to run DA NE within a ±20 MeV interval apart from the φ mass peak value.

16.2
The π + π − contribution to g − 2: comparison with data An interesting piece of information comes from comparing our (VMD) estimates derived from global fitting with the corresponding estimates provided by the various experimental groups. Table 7 displays the published experimental results concerning the contribution of the 0.630-0.958 GeV/c region to a μ (ππ). We first list the three important results from CMD-2 and SND; as we also use the data sets from OLYA Table 6 The contribution to 10 10 a μ (ππ) from the invariant mass region 0.630-0.958 GeV/c using KK data sets under various conditions. All π + π − π 0 data have been excluded from fit. FSR corrections have been performed  Table 7 The various published estimates of the contribution to 10 10 a μ (ππ) from the invariant mass region 0.630-0.958 GeV/c. The quoted averages always refer to all experimental results displayed in the preceding lines. The line "OLD" information refers to our aver-age performed using the data sets collected before those of CMD-2 and SND (see text). Our fit solutions A and B are derived using the τ spectra from [40][41][42]. KLOE and CMD [80], we also give at the line flagged by "OLD" our average using these data sets together with those from NA7 [111], TOF [112], M2N [113], DM1 [81], all collected before those from [75][76][77][78].
The third data column provides, first, our average derived using the data sets from [75][76][77][78] and, next, also those including the older data sets referred to just above. Our results are directly comparable with these as we do not yet use ISR data.
Both solution A and solution B results favorably compare with the scan (ππ) data averaging as the uncertainty is reduced by a factor close to 2.
The following lines of Table 7 display, for information, the experimental results derived from the data sets collected using the ISR method and the global average of the ISR and scan data.
One should stress that our results for a μ (ππ), derived excluding the ISR data, provide information already comparable in precision to those obtained using them. This motivates to examine the ISR data in view of including them into the fit procedure.
One may also compare our estimates with the weighted average of the τ data [40][41][42] which gives 10 10 a μ (ππ) = 365.21 ± 2.67 exp in the reference region, including FSR corrections; applying the ρ-γ corrections proposed in [16], this becomes 10 10 a μ (ππ) = 361.66 ± 2.67 exp and provides 10 10 a μ (ππ) = 361.15 ± 1.76 exp when averaged with the e + e − data. This indicates that examining the idea proposed in [16] in a wider context is an interesting issue. Indeed, this could lead to another successful VMD-like model and, therefore, may contribute to a motivated evaluation of the model dependence of a μ estimates.
As a summary, one may conclude that our global model provides a good determination of the contribution to a μ (ππ) from the invariant mass region 0.630-0.958 GeV/c. The accuracy of our VMD estimates is found much improved compared to direct averaging of the experimental data and their central values are found consistent within uncertainties. By including ISR data at a later stage, the precision of the result might be further increased.

Hadronic contribution to g − 2
In Table 8, one displays the contribution of each of the examined channels to a μ from their respective thresholds up to 1.05 GeV/c, i.e. slightly above the φ peak.
The first two data columns show the results corresponding to the so-called configurations/solutions A and B. These have been derived by fitting the data sets referred to in the preceding Sections and the motivation to consider both solutions valid can be emphasized from Table 3.
The last two data columns exhibit the averages of experimental data for each of the measured channels submitted to the global fit. These differ by excluding (third data column) or including (fourth data column) in the averaging the ISR data sets collected by KLOE [17,19] and BaBar [18] for the π + π − final state. As we have excluded for now the ISR data from our analysis, the gain due to the global fit can be directly inferred by comparing with the third data column; nevertheless, it is interesting to compare the accuracy of solutions A and B to the averages derived using the high statistics ISR data.
As expected, the improvement generated by the global fit affects all the channels considered and is always a factor of 2 or more (see the π + π − π 0 channel) better than the Table 8 Contributions to 10 10 a μ from thresholds up to 1.05 GeV/c The experimental errors merge the reported statistical and systematic uncertainties in quadrature. FSR effects (3.43 × 10 −10 ) have been in-cluded into the π + π − contribution. The first two data columns display our fit results and the last two data columns report the direct numerical integration of the relevant data average of the same data. The first line even shows that our accuracy is comparable-actually slightly better-than the average derived using the ISR data.
It is interesting to note that the sum of all contributions for solution B is in accordance with the result expected from the standard sum as reported in the third (or fourth) data column. Solution A, instead, gives a smaller sum than the experimental average of the same data; the distance is 2.97 10 −10 , i.e. 1.6σ theor. or 0.7σ exp .
It is interesting to examine the individual channel contributions. Those from the π 0 γ and ηγ channels, as calculated from data, rely on pretty poor statistics and generally cover restricted energy ranges [82][83][84][85][86][87] (see Subsect. 9.2); instead, our model results are estimated (significantly) larger and cover precisely the full energy range from thresholds to 1.05 GeV. This especially concerns the region in between the ω and φ peaks.
Our model estimates for the π + π − π 0 and K + K − channels are found smaller than the experimental averages at the 1 or 2 σ exp levels, while the K 0 K 0 contribution corresponds to the experimental expectation. This confirms the need for a better experimental knowledge of all annihilation channels in the φ region.
The first data line in Table 9 reports the results derived from fits with our global model. The second line ("missing channels") provides the experimental averaged contribution to a μ from the channels unaccounted for within our model (the 4π , 5π , 6π , ηππ and ωπ final states). This has been computed using the trapezoidal integration rule. As the corresponding data are sparse below 1.05 GeV, this estimate might have to be improved.
The line "Total Model" provides the estimate of the full hadronic vacuum polarization (HVP), merging our model results with the additional listed contributions.
In order to illustrate the impact of τ data, we present separately the fit results derived when including or when excluding the τ data sets from the fitted data sets, keeping for the rest the configurations leading to solutions A and B as previously defined.
Including τ data sets results in an increased value of the hadronic VP by 3 × 10 −10 . This will be commented on below. One also remarks that our uncertainties are comparable to the experimental one, even if our estimates are penalized by having-provisionally-discarded the ISR data. Our estimates also compare favorably with the revised estimate excluding all ISR data given by [13]: a μ (e + e − ) = (690.9 ± 5.2 exp+rad ± 0.7 QCD ) × 10 −10 .
16.4 The anomalous magnetic moment of the muon a μ Table 10 displays our final results concerning a μ . We still report on the results derived in the fit configurations A and B, using or not the τ data in the fit procedure. The leadingorder (LO) hadronic VP discussed in the previous Subsection is reminded in the first line. In order to yield our estimate of a μ under the various quoted configurations, one should add the effect of higher-order hadronic loops taken from [16], the light-by-light contribution [5]; we took the latest estimate of the pure QED contribution 46 [3] and the Table 9 Hadronic VP contributions to 10 10 a μ with FSR corrections included. Numbers within brackets refer to respectively statistical and systematic errors. Numbers within square brackets are the total uncertainties electroweak (EW) contribution is taken from [4]. Summing up all these, one obtains the values given as "Total Theor." which should be compared with the average [1] of the different measurements for a μ , recently updated [2]. The difference between our theoretical estimates and the experimental average [2] is finally given together with their respective statistical significance. The significance of this difference varies between 4.07σ (solution B including τ 's) to 4.65σ (solution A excluding τ 's). The difference between including τ 's and excluding them is a 0.4σ effect. [13] provides an estimate excluding the KLOE data [17]-and the more recent ISR data sets not available at that time-reaching a difference with the BNL average [2] of (30.1 ± 8.6) × 10 −10 , a 3.5σ significance. Our least significant estimate (solution B including τ 's) is, instead, 4.07σ .  [14]. The following entry is the estimate given in [16] which combines e + e − and τ data (after correcting for the ρ 0 -γ mixing). The last entry [115] is derived including the ISR data (HLMNT11); this is the latest result using the final KLOE [19] and BaBar [18] data.
We have also displayed the latest result [13] derived excluding ISR data which directly compares to ours. This indicates that the improvement provided by the global fit method corresponds to increase the discrepancy of the BNL measurement [2] with the Standard model prediction by 0.6÷0.8σ . Therefore, the discrepancy starts reaching an interesting significance. 16.5 Influence of data set choices on the estimate for a μ In order to derive our estimates for a μ , we have defined a paradigm, unusual in this field. Indeed, one usually performs the average using all data sets contributing to a given final state in isolation; the prescription used is the S-factor technics of the Particle Data Group. However, this supposes the simultaneous handling of statistical and systematic uncertainties. The most common way of performing this handling is to use as weights the quadratic sum of statistical and systematic uncertainties [9].
In our approach, especially in this paper, the underlying paradigm is different and can be formulated in the following way: • All different channels are correlated by their underlying common physics and an Effective Lagrangian approach is presently the best tool to deal with the non-perturbative QCD regime. • All data sets, covering or not the same physics channel are considered by taking into account the peculiarities of their uncertainties as reported by the experimental groups. There is, in principle, no real difficulty in order to deal with statistical uncertainties. It is commonly assumed that uncorrelated systematics and statistical uncertainties could be added in quadrature and we followed this rule. Other systematics involving bin-to-bin or experiment-to-experiment correlations should be treated as such; the method is standard 47 and has been sketched in Subsect. 9.7. • The Lagrangian model should allow for a good description of a large number of data sets in as many different physics channels as possible. The goodness of the global fit should be accompanied by a good description of each group of data sets-ideally each data set. As tag for this property, we choosed the χ 2 /n points value for each data set group; this tag should not too much exceed 1. Referring to our case, the π + π − , π 0 γ , ηγ physics channel data and the reported partial width decays already represent an acceptably good starting point, allowing a critical examination of the data associated with further additional channels. 47 In the scan experiments we deal with in the present paper, all reported correlated systematics can be considered as global scale uncertainties for which the standard method applies. For ISR experiments [17][18][19], the situation is different as several independent sources of systematics are defined which, additionally, vary all along the spectra. The standard method can be extended to this case [46]; however, it should better be reformulated in a way which avoids introducing as many scale factors to be fitted as sources of different systematics. Indeed, this may produce fit instabilities and, on the other hand, one has to deal with correlations between physics parameters and these scale factors which may be uneasy to handle.
• Including a new data set, or a new group of data sets, should not result in a significant degradation of the already accounted for data sets. This should be observed at the global level and at the local levels (i.e. for each group). Following from the analyses in Sects. 10 and 11, peculiarities of their fit behavior led us to discard from our global fit the K + K − data set and one of the π + π − π 0 data sets provided by SND. This turns out to require that the (large) set of data samples considered be statistically self-consistent: Only 2 data sets out of 45 did not pass this consistency criterium.
At this point, given the (broken) Lagrangian one uses, the selection criteria are only the global fit quality and the "local" (data set specific) fit properties reflected by the various χ 2 /n points values, discarding any possible consequence for the value for a μ . With Solutions A and B, one has also avoided any kind of data set reweighting by discarding the two data sets exhibiting some faulty behavior compared to the rest.
Nevertheless, it is a simple exercise to switch on the two discarded SND data sets within our fitting code. For information, this leads to a μ = (a μ ) exp − (a μ ) th = (34.00 ± 8.21) × 10 −10 , a 4.14σ effect. However, this is associated with an exceptionally poor global fit probability (1.75%) and to χ 2 π + π − π 0 /n points = 331/212 = 1.56 and χ 2 K + K − /n points = 93/62 = 1.50. Interestingly, and somewhat unexpectedly, the χ 2 /n points for the other data sets are practically unchanged compared to Table 3, except for the decay data set account which is sharply degraded: χ 2 decays /n points = 20.5/10 2. This may reflect that our broken HLS model is so sharply constrained that poor data sets are mostly reflected by poor global fit probabilities.
A tag value of χ 2 /n points = 1.3, as yielded for the chosen π + π − π 0 final state data, is on the border of what could look reasonable to us (see third data column in Table 3). Nevertheless, compared with χ 2 /n points = 1.1 (see second data column in Table 3), it looks acceptable; however, this corresponds to an increase by 30 units of the absolute magnitude of χ 2 π + π − π 0 , when introducing the selected kaon data. One may, indeed, consider that this indicates some tension within the φ region data calling for a closer experimental examination which can be performed at the existing facilities covering the φ region.
Awaiting for better data in the φ region, we have been left with two challenging solutions: Solution A which uses all the data sets we have considered as secure, and solution B obtained by removing all π + π − π 0 data sets above the KK threshold.
16.6 The differential effect of the various τ data samples In view of the discussions above, we have chosen to display all our final results for a μ , in the fit configurations corresponding to solutions A and B. On the other hand, as can be read off Table 3, at the fit properties level, one can consider that the so-called e + e − − τ puzzle is over.
However, one still observes a (2÷2.5) × 10 −10 increase of the returned values for a μ produced by the τ data. As stated already above, the τ data are essential in order to return a reasonably precise value for our fit parameter 48 Σ V . Therefore, the shift attributable to the τ data can be considered as a normal consequence when fitting a model with a more constraining set of data samples.
Nevertheless, Table 3 indicates that the χ 2 /n points are sensitively different for ALEPH ( 0.43), CLEO ( 1.26) and BELLE 49 ( 1.77). This difference of fit quality leads us to examine the effects of removing the CLEO data sample and/or the BELLE data sample for our fitted data set.
When keeping only the ALEPH data sample, we get a μ = 38.47 ± 8.22 (a 4.68σ significance) and a μ = 36.81 ± 8.90 (a 4.13σ significance) for respectively solutions A and B. As can be seen from Table 10, these strikingly resemble the corresponding values for a μ derived when keeping only e + e − data in our fit procedure (i.e. excluding all τ data). In these peculiar configurations, the ALEPH data fit quality which was already very good (χ 2 /n points 16/37), becomes impressively better (χ 2 /n points 4/37).
Going a step further, we have examined the effect of considering only ALEPH and CLEO data. In this case, our fit returns a μ = 36.02 ± 8.22 (4.38σ significance) and a μ = 34.74 ± 8.26 (4.21σ significance) for respectively solutions A and B. One can check with Table 10 that these values become closer to their partners when fitting excluding τ samples.
Therefore, using only the τ data samples from ALEPH [40] and/or CLEO [42] returns values for a μ consistent well within errors with those derived using only e + e − data. The slightly different behavior of BELLE data may be related with the normalization issue sketched in footnote 49.

On the significance of the HLS value for a μ
In view of the considerations developed in the two preceding Subsections, one can certainly consider that the most conservative estimates for a μ are those derived while including τ data as they are reported by ALEPH, BELLE and CLEO. This corresponds to the information provided in the first two data columns of Table 10. 48 The numerical accuracy of the scan e + e − data alone does not permit a precise determination of Σ V which is returned by MINUIT with large errors. 49 Leaving free the absolute normalization of their dipion spectrum improves the stand-alone fit of the BELLE Collaboration [41] from 80/52 to 65/51. This corresponds to a best normalization of 1.02±0.01. Such a re-normalization of their absolute scale has some influence on the value for a μ . One should remind that we do not have any longer fitted rescaling factors in our fitting functions. This means that the disagreement between the BNL measurement [2] and the Standard model prediction for a μ lays in between 4.07 and 4.33σ . Moreover, from our analysis of the differential effects of the various available τ data samples, one may consider these bounds as conservative and that the significances in the right part of Table 10 cannot be discarded.
In view of this, in the perspective of taking into account relatively poor data set group, one has rerun our code in order to get the solution when weighting the contributions of 50 : • all π + π − π 0 data in our global sample by 179/232.41, • the BELLE data sample by 19/32.31, • the CLEO data sample by 29/36.48, in the global χ 2 while leaving the other weights (all equal 1) unchanged. This turns out to rescale globally the uncertainties associated with the corresponding data sets by the inverse of these weights, assuming that their relatively poor quality is only due to an overall underestimate of the uncertainties by a factor of respectively 1.14 (π + π − π 0 ), 1.30 (BELLE) and 1.12 (CLEO). This may look as a way to infer some sort of S-factors inside the global fit procedure.
Going a step further, another check may look appropriate. As the contributions of the π + π − π 0 , BELLE and CLEO data to the total χ 2 have been weighted in order to reduce their influence, one can do alike with those groups of data which exhibit too favorable individual χ 2 's. Still referring to fitting with configuration A, this turns out to weight the "Old Timelike" data by 82/56.61, the π 0 γ data group by 86/68.37, the ηγ data group by 182/123.31, the ALEPH data by 37/15.92 while keeping unit weights for the "New Timelike" and both KK data groups. This leads to an hadronic VP of (685.00 ± 4.58) × 10 −10 and to a μ = (36.25 ± 8.21) × 10 −10 corresponding to a 4.41σ discrepancy. This is almost identical to the value found with Solution B, excluding τ 's, as can be seen from Table 10.
Therefore, these exercises enforce our conclusion that the most conservative value for a μ exhibits a discrepancy of 4.07σ and values as large as (4.30÷4.50)σ are not unlikely. 50 The weights used in this Subsection refer to partial χ 2 's obtained by fitting under Configuration A with assuming c 3 = c 4 ; it is the reason why they slightly differ from the corresponding numbers given in Table 3. 51 We have also made a fit leaving free scale factors affecting the covariance matrices of the 3-pion data as a whole, of the BELLE and CLEO data. The hadronic VP we get is (686.73 ± 4.49) × 10 −10 , quite similar to this value.

Conclusion and perspectives
Several aspects should be emphasized. They can be grouped into two items: Low energy hadronic physics description and g − 2 related topics.
Concerning the first item, the present study indicates that the HLS model suitably broken is able to encompass most low energy physics in an energy range extending up to the φ meson mass. More precisely, among the non-baryonic possible final states, one covers 52 most channels with multiplicity n < 4.
More precisely, equipped with the so-called upgraded direct symmetry breaking-in the u, d and s sectors-and including the mixing of neutral vector mesons produced at one-loop, the HLS model accounts quite satisfactorily for all the examined physics pieces of information. This covers the 6 annihilation channels having significant cross sections up to the φ meson mass and a few more spectra like the dipion spectrum in the τ decay and, also, an additional list of partial width decays. Previous studies [37,46] have also shown that the dipion spectra in the η/η → ππγ decays fall inside the scope of the HLS model.
It is an attractive feature of this framework to exhibit a parent character between the long reported issues represented by the e + e − − τ and the φ → KK puzzles: Indeed, it is the same breaking mechanism implemented in the L A and in the L V pieces of the HLS Lagrangian which provides a solution to both. It permits-together with the s-dependent vector meson mixing-to finalize the consistency of the e + e − and τ physics and to reproduce the branching fraction ratio φ → K + K − /φ → K 0 K 0 . This is materialized by a satisfactory simultaneous fit of both e + e − → KK cross sections and of the pion form factor in both e + e − annihilation and τ decay. The upgraded model thus provides a tool allowing a simultaneous treatment of a large number of experimental spectra. It also permits a critical analysis of the fit behavior of any data set in consistency with the others. Then, one is in position to discard motivatedly some data samples which do not behave satisfactorily within a global fit procedure and could then put some shadow on derived numerical results. We have shown that such data samples are only few: 2 out of the 45 considered spectra. It should be stressed that discarded data sets are always identified because of their full redundancy with some other data sets, which are found to behave normally within the global model; stated otherwise, this removal is not expected to produce a bias and, a contrario, any effect resulting of keeping them is suspicious.
The model provides a tool which has the virtue of exhibiting the physics relationship between the various physics channels. Within the global fit procedure involving the data on each channel, the model parameters yield a better accuracy which propagates to all the reconstructed pieces of information, especially the photon hadronic vacuum polarization and, thus, improves significantly g − 2 estimates.
Indeed, we have shown that the various components of the HVP yield central values in accordance with expectations and an uncertainty improved by a factor of 2 quite uniformly within the fit range. This has been shown for the π + π − , π 0 γ , ηγ , π + π − π 0 , K + K − and K 0 K 0 channel contributions up to 1.05 GeV. Up to this energy, these channels represent altogether more than 80% of the hadronic VP and one of the two dominant sources of uncertainty. 53 In order to figure out the gain in terms of statistics, one can make the following statement: considering globally the existing data sets is equivalent to having ×4 more statistics simultaneously in each of the considered channels without any increase of the systematics. Therefore, considering additionally the high statistics ISR data leaves some room for improved estimates of the HVP, provided the dealing with systematics can be reasonably well performed. One should nevertheless stress that the global method we advocate, used with only the standard scan data samples provides already as good results as all scan and ISR data using the standard numerical integration of the experimental cross sections.
One may also try to figure out the improvement expected from including the high statistic ISR data samples [17][18][19] within the fit procedure. Being optimistic, one may think that the uncertainty on the HVP contribution up to 1.05 GeV could be divided by 2, from 2 × 10 −10 (see Table 8) to 1 × 10 −10 . Let us also assume that the ISR data samples will not rise unsolvable bias problems. Taking into account the rest of the HVP, which carry an uncertainty of 4 × 10 −10 (see Table 9), the uncertainty on the full HVP would decrease from 4.60 × 10 −10 (see Table 9) to 4.25 × 10 −10 . Using the information collected in Table 10, the total uncertainty on a μ would decrease from 5.30 × 10 −10 to 5.00 × 10 −10 and the uncertainty on a μ would decrease from 8.20×10 −10 to 8.00×10 −10 . This may look a marginal improvement; the reason for this is the large value for the systematics generated by hadronic HVP in the region 1.05÷3.10 GeV (see Table 9), which thus becomes a prominent issue for future significant improvements. 54 53 The other dominant error comes from the hadronic VP between 1.05 and 2 GeV. 54 Actually, even if the uncertainty on the HVP contribution coming from the energy region up to 1.05 GeV vanishes, this would not en-However, this is not the end of the story. In the course of the paper, and this is well expressed by Tables 9 and 10, we saw that below 1.05 GeV systematics may produce significant shifts of the central values for the HVP and thus for a μ . This was observed, for instance, in the A and B configurations, where the shift for the HVP-and for a μ -amounts to 2.00 × 10 −10 (see also Subsect. 16.6). Because of this, there is still valuable experimental work to do also in the sub-GeV domain to decrease and/or better understand systematic errors. More precisely, a better experimental knowledge of all channels in the φ mass region-0.95÷1.05 GeV-may result in improving quite significantly our estimate on g − 2 and in resolving some of the ambiguities discussed in the main text. As stated above, the information in this mass region has an important influence down to the threshold regions. This is certainly within the scope of existing machines and detectors. 55 What are the prospects for the future? A new muon g − 2 experiment at Fermilab is expected to come into operation in 5 years from now. The accuracy is expected to improve to 0.14 ppm from its current 0.54 ppm. This also requires a factor 4 improvement of the hadronic vacuum polarization. As demonstrated by our analysis, it is possible to improve the low energy part up to and including the φ by a systematic application of effective field theory methods in form of a resonance Lagrangian approach. However, as mentioned above, the main effort will be required in the range above the φ up to about 3 GeV. In this range, major progress is expected from CMD3 and SND at VEPP 2000 at Novosibirsk, from BESIII at Beijing, as well as from exploiting additional yet unanalyzed ISR data from BaBar and Belle. Within the 5 years available until a new experimental result for a μ will be realized, lattice QCD is expected to be able to produce results which are competitive with standard evaluations based on data. This also would provide important cross checks for the present results and, more generally, for the effective Lagrangian approach.
For now, one can conclude that the paradigm represented by a global model which encompasses the largest possible set of data indeed results in a highly significant improvement of the photon HVP uncertainty and of the uncertainty on g − 2. As the global model allows to detect problematic data sets susceptible of generating biases, it must be accompanied by the most accurate possible treatment of the reported experimental systematics. tail a significant improvement of the global uncertainty for a μ ! Stated otherwise, reducing the HVP error in the region from threshold to 1.05 GeV from 4 × 10 −10 to 2 × 10 −10 has much more dramatic effects than reducing it from 2 × 10 −10 to 1 × 10 −10 . This is a pure algebraic effect following from having to perform quadratic sums for final uncertainties. 55 One may remark that scan data for the e + e − → π + π − cross section in the φ region are still not available.
Taking into account the ambiguities generated by a limited number of data sets, the most conservative estimate for the hadronic vacuum polarization leads to a significance for a non-zero a μ of 4.1σ . Solving these ambiguities discussed in the main text may result in a significant increase of this conservative bound.
Open Access This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.

Appendix A: The full HLS non-anomalous Lagrangian before loop mixing
The non-anomalous Lagrangian of the Hidden Local Symmetry Model can be written: in order to split it up into convenient pieces. Removing the pseudoscalar field kinetic energy term, which is canonical, one has: and by: In the ω R 1 P γ sector, one has: 1 GN c g 0 V R 1 P 0 γ , G = − egc 3 4π 2 f π (113) for each V R 1 = ρ R 1 , ω R 1 , R 1 and P 0 = π 0 , η, η . The g 0 V R 1 P 0 γ can be found in Appendix E in (108), (109), (110). The functions H P 0 V i occurring in (66) provide the couplings of the physical vector fields to a photon and a neutral meson. They are given by: These definitions help in writing the cross sections in a way similar to those in [46]. When expanded, the H P 0 V i functions may contain contributions of order greater than 1 in some of the breaking parameters. These higher-order contributions are irrelevant and can be dropped out.

Appendix H: The functions N i (s) in e + e − → π 0 π + π − annihilations
The amplitude for the transition γ * → π 0 π + π − is much simply expressed in terms of the following complex functions: The integration limits can be found in [73]; they are also reminded in [46]. The Kuraev-Siligadze kernel is: