BHLS$_2$, a New Breaking of the HLS Model and its Phenomenology
M. Benayoun L. DelBuono F. Jegerlehner









   =0cm
   =0cm
  =15.9cm
 
   =-1cm
  
    =23cm



    ∂↔→←η^'ωγγγ i
        MeV ''[ 0.mm]0.mm6.mm[-3.mm]0.mm7.mm[-3.mm]0.mm9.mm[-1.mm]0.mm8.mm[-3.mm]0.mm3.mm[-3.mm]0.mm2.mm[-3.mm]0.mm4.mm[-3.mm]0.mm12.mm[-3.mm]0.mm9.mm[-1.mm]0.mm5.mm[-5.mm]0.mm12.mm
             e F M
               
            
            
             
            
            
              BHLS_2, a New Breaking of the HLS Model and its Phenomenology
                
            M. Benayoun^a, L. DelBuono^a, F. Jegerlehner^b,c
            ^a LPNHE des Universités Paris VI et Paris VII, IN2P3/CNRS, F–75252 Paris, France 
            ^b Humboldt–Universität zu Berlin, Institut für Physik, Newtonstrasse 15, D–12489 Berlin,
            Germany 
            ^c Deutsches  Elektronen–Synchrotron (DESY), Platanenallee 6, D–15738 Zeuthen, Germany
            
                December 30, 2023
            ===================================================================================================================================================================================================================================================================================================================================
            
            
            Previous studies have shown that the  Hidden Local
            Symmetry (HLS) Model, supplied with appropriate
            symmetry breaking mechanisms, provides an Effective Lagrangian
            (BHLS) able to encompass a large number of processes within a unified
            framework. This allowed to design  a global fit procedure
            which provides a fair simultaneous description of the e^+ e^-
            annihilation  into 6 final states (π^+π^-, π^0γ,
            ηγ, π^+π^-π^0, K^+K^-, K_L K_S),
            the dipion spectrum in the τ decay and some more
            light meson decay partial widths. In this paper, additional breaking schemes are
            defined which improve the BHLS working and extend its scope
            so as to absorb spacelike processes within a new framework (BHLS_2).
            The phenomenology previously explored with BHLS is fully
            revisited in the BHLS_2 context with special emphasis on the ϕ mass region  using
            all  available data samples.  It is  shown that BHLS_2 addresses
            perfectly the close spacelike region covered by NA7 and Fermilab data; it
            is also shown that the recent Lattice QCD (LQCD) information on the pion form factor
            are accurately  predicted by the BHLS_2 fit functions derived from fits to
             only annihilation data. The  contribution
            to the muon anomalous magnetic moment  a_μ^ th of these
            annihilation channels over the range of validity  of BHLS_2
            (up to ≃ 1.05 GeV) is updated within the new BHLS_2 framework
            and shown to strongly reduce the former BHLS systematics. The
            uncertainty on  a_μ^ th(√(s)< 1.05 GeV) is much
            improved compared to  standard approaches relying on direct
            integration  methods of measured spectra. Using the BHLS_2
            results, the leading order HVP contribution to the muon anomalous moment
            is a_μ^ HVP-LO= 686.65 ± 3.01 +(+1.16,-0.75)_ syst in units
            of 10^-10. Using a conservative
            estimate for the light–by–light contribution, our evaluation for the muon anomalous
            magnetic moment is a_μ^ th= [11 659 175.96 ± 4.17
            +(+1.16,-0.75)_ syst] ×
            10^-10. The relationship between the dispersive and LQCD approaches to the
            ρ^0-γ mixing is also discussed which may amount to a shift of
            δ a_μ[ππ]_ργ=+(3.10± 0.31) ×10^-10 at LO+NLO,  
            presently treated as additional systematics. Taking also this shift into account,
            the difference a_μ^ th-a_μ^ BNL exhibits a significance not smaller
            than 3.8 σ.
            
            
            § INTRODUCTION
             The Standard Model is widely recognized as the (gauge) theory which unifies
             the whole  realm  of weak, electromagnetic and strong interactions among
            quarks, leptons and the various gauge bosons (gluons, photons, W^±, Z^0).
            For the physics processes – and quantities – involving strong interactions, QCD is at work
            under two different regimes  tightly connected with the energy involved and  the onset of the perturbative
            regime is  a priori expected  to occur at high energies. However, as clearly shown by
            the data recently collected by  KEDR <cit.> on the ratio
            R(s)=σ(e^+e^-  hadrons)/σ(e^+e^- μ^+μ^-),
             the {u,d,s} sector  of this ratio reaches the
            perturbative regime at energies as low as  ≃ 2.0 GeV (see   <cit.>,
            for instance); above this energy,  the observed departures from perturbative QCD predictions – including, of course,
             the c c and b b  threshold effects therein –
            are only spikes and narrow bumps associated with the charmonium and bottomonium states
            which require an additional  specific treatment in the mass range where they are located.
            
            However, in  the low energy region  where the non-perturbative
            regime of QCD is involved, getting   theoretical predictions able to compete with the accuracy of
            some important experimental measurements  may be challenging. This is, in particular,  the issue met
            with the photon hadronic vacuum polarization (HVP), which contributes  importantly  to the muon
            anomalous magnetic moment a_μ, one of the best measured particle properties. The
            hadronic part of  a_μ is related to the so-called R(s) ratio and is given (at leading order) by :
            
            a_μ^ had =  [ α m_μ/3 π ]^2
            ∫_s_0^∞ds/s^2 K̂(s) R(s) ,
            
            where s_0=m_π^0^2 is the lowest hadronic threshold and K̂(s) is a known smoothly
            varying positive function <cit.>. The 1/s^2 factor strongly enhances the low energy contribution
            of the spectrum,  i.e. just in the validity domain of non-perturbative
            QCD.  Fortunately, as obvious from  the KEDR data <cit.>, one can undoubtedly consider
            that the non-perturbative regime does not extend to energies larger than ≃  2.0 GeV.
            
            So,  the real issue with predictions for objects such as a_μ is to get precise estimates
            of the effects covered by the non-perturbative regime of QCD. Chiral Perturbation Theory
            (ChPT)<cit.>, the low energy limit of QCD,  is of limited help for the present purpose
            as its realm does not extend much beyond the 400 ÷ 500 MeV region and misses the quite
            important  meson resonance mass region.
            
            The most promising approach to the non-perturbative regime of QCD is certainly  Lattice QCD (LQCD)
            which already provides valuable information at low energies <cit.>.
            The  recently derived LQCD evaluations of a_μ^ had<cit.>  have been found in accord with the
            experimental measurement  performed at BNL  <cit.>. However, the magnitude
            of the reported uncertainties is by far too large to fruitfully compare with
            the already existing BNL datum and,  a fortiori, with the measurements expected from
            the  Fermilab experiment <cit.>, already running, or from
            the experiment  planned to start later on at J-PARC <cit.>,
            as both are expected to improve the uncertainties by a factor of 4. So much progress remains
             to be done before getting satisfactory uncertainties from LQCD.
            
            
             
            
            
            This leaves room  for low energy Effective Resonance Lagrangian Approaches,
            which can contribute to improve the knowledge of a_μ^ had by providing
            a good description of the rich amount of experimental data collected in the
            timelike region.
            Among the richest possible Lagrangians, the elegant Hidden Local Symmetry (HLS)
            Model <cit.> is worth to be considered;
            it has been proven to be equivalent to  RχPT <cit.>,
            provided consistency with the QCD asymptotic behavior is incorporated. It thus follows that the
            HLS model is  a motivated and constraining QCD rooted framework, moreover, easy enough to handle
            in phenomenological applications.
            
            The original HLS Model deals with the lowest vector meson nonet and provides a framework
            for hadron production in e^+e^- annihilation, naturally bounded by the ϕ mass region
            – i.e.  up to ≃ 1.05 GeV.
            It thus represents  a  tool giving a handle on
            a mass region contributing for ≃ 83% of the total muon hadronic VP.
            The region extending from just above the ϕ meson mass to 2 GeV only contributes
            for 7%, slightly less than the [2  GeV, ∞]  region (10%).
            
            As such, the non-anomalous HLS Lagrangian  <cit.> sets up a unified framework
            which encompasses the e^+e^- π^+  π^-, e^+e^-  K^+  K^- and
            e^+e^-  K^0 K^0 annihilation  channels;
            the τ^±π^±π^0 ν spectrum also belongs to the same framework, allowing naturally
             a simultaneous treatment of all the mentioned annihilation and decay processes. On the other hand,
            the HLS Model possesses an anomalous sector <cit.> which also brings
             the e^+e^- π^0  γ, e^+e^- ηγ and
            e^+e^- π^+  π^- π^0 annihilation channels inside
            the same framework, together with some radiative
            decay modes and the dipion spectra in the η/π^+  π^- γ decays
            as shown in former studies <cit.>.
            The annihilation channels just listed exhaust almost completely the processes
            contributing to a_μ^ had below  1.05 GeV; indeed the missing
            channels[The e^+e^- γ channel has been considered
            in <cit.> and  found to contributes
            to a_μ^ had at the 10^-12 level only and can thus be safely discarded.]
            (4 π, 2 πη, …) contributes only  ≃ 2 of the full HVP.
            Such a broken HLS (BHLS) Model has already been built up <cit.>
            and shown to provide a pretty good simultaneous  description of almost all data samples
            covering the six channels listed above up to ≃ 1.05 GeV. This work has also proved that
            there was no contradiction between the  e^+e^- π^+  π^- and
            τ^±π^±π^0 ν spectra.
            
            
            The non-anomalous <cit.> and anomalous <cit.>
            sectors of the HLS Model thus open a unified framework able to encompass
            a large corpus of data and  physics processes. However, as such,
            the HLS framework – with only the universal vector coupling g as
             free parameter – cannot pretend to provide a satisfactory simultaneous
             description of the wide ensemble of high statistics  data samples collected by
            several sophisticated experiments in several annihilation channels.
            
            In order to achieve such a program,  the HLS model  must be supplied
            with appropriate symmetry breaking mechanisms not present in its original
            formulation <cit.>. A first successful attempt has been done in
            <cit.> which has set up a model (BHLS)
            based on a breaking mechanism originally proposed in
            <cit.>, hereafter named BKY.
            
            The present work reports on
            a new breaking scheme for the HLS model which aims at improving the behavior
            of the form factors in the dipion threshold region and also in the ϕ region where
            the original BHLS <cit.> meets some difficulty leading to introduce additional
            systematic uncertainties in the evaluation of the muon g-2. The model (BHLS_2)
             presented here will be shown to widen the scope of BHLS, in particular to the spacelike
             domain and improve the BHLS prediction for the muon g-2. Indeed, using
             variants of BHLS_2, together with the original BHLS should allow evaluating the model
             dependence effects in the estimates of derived physics quantities, noticeably the
             photon HVP.
            
            The layout of the paper is as follows. Sections <ref> gives  a brief reminder of
            the original HLS model while Section <ref> briefly reminds the BKY mechanism, still involved
            in BHLS_2. Sections <ref>   and <ref>  introduce a new breaking mechanism at the
            level of the covariant derivative which is the basic object leading to the original HLS model.
            Section <ref>  describes the O(p^4) terms <cit.> of the HLS framework
            which are included in the  BHLS_2 Lagrangian.
            
            At this step, one has at hand the first variant of the BHLS_2 model, named Basic Solution (BS),
            and, relying on its vector meson mass term,  Section <ref> defines its parameter
            properties.  Remarking that there is no fundamental reason why the neutral vector fields involved in
            physics processes should be their ideal combinations (ρ_0^I, ^I, ϕ^I), one allows
            for combinations of these via a mechanism called Primordial Mixing described in Section <ref>.
            This gives rise to a second variant of the BHLS_2 model, named Reference Solution (RS). In Section
            <ref>,  one first reminds the parametrization of the propagators for the  and ϕ
            fields <cit.> in connection with form factor behaviors at s=0;  this part of the work also shows
            that the breaking of nonet symmetry in the vector meson sector must be accompanied by
            breaking of SU(3) of the same intensity.
            
            The dynamical breaking of vector meson  <cit.> is
            revisited in Section <ref> and the role of anomalous loop corrections
            is emphasized, especially in order to connect smoothly the timelike and spacelike
            branches of the pion form factor.
            The derivation of the pion and kaon form factors is the subject of Section <ref>.
            The anomalous sector is fully analyzed in Section <ref> and the cross sections for
            the e^+e^-  (π^0/η) and e^+e^- π^+  π^- π^0 annihilations
             are derived.
            
             Altogether, the above-listed Sections and the Appendix fully report on  the BHLS_2 model properties
             and tools. Section <ref> provides a comprehensive discussion of the available data
             samples falling into the BHLS_2 scope and their peculiarities. Relying on several preliminary studies,
             one comments here on the 3 discarded data samples; one should note that the number of data samples
             found to accommodate with each other within the BHLS_2 framework exceeds now 50 and the total
             number of data points reaches 1237.
            
            
            Sections <ref> and <ref> give a full description of all aspects of the fits
            to the 6 annihilation channels performed within the global BHLS_2 framework. Moreover,
            Section <ref> also provides  a comprehensive analysis
            within our global framework of the existing K K data samples.
            The aim is to substantiate the issue between the BaBar, CMD-3 and CMD-2
            spectra revealed by fits. A summary of the various fit properties is presented in Section
            <ref>.
            
            The description of the spacelike region for the pion and kaon form factors as coming
            from BHLS_2 is the purpose of Section <ref>.
            The existing model-independent data
            for the π^± and K^±  form factors are found to naturally  accommodate our
            global framework, giving support to the low energy behavior predicted by BHLS_2
            for all meson form factors.
            This fair agreement extends to the π^± form factor data provided by several Lattice
            QCD (LQCD) groups.
            One may infer from these data and from other model-dependent data
            that  BHLS_2 should almost certainly apply
            down to ≃ -1 GeV^2. The prediction for the K^0  form factor is shown to cope
            with expectations at s≃ 0. Information on the charge radii of the π^± and K^±
            mesons
            is provided. Section <ref> examines the P-wave ππ phase shift prediction and
            shows its agreement with data and other predictions, including our former BHLS, over
            an energy range extending up to the ϕ mass region.
            
            The values found for some parameters, noticeably the a HLS parameter and some BKY
            breaking parameters, deserve a specific account given in Section <ref>.
            
            The main motivation of the present work is the study of the BHLS_2 predictions
            and their consequences on the muon g-2.  This topic is addressed in Section <ref>.
            The part (≃ 83%) of the muon HVP contribution to a_μ covered by BHLS_2 is 
            estimated
            in several contexts to derive a_μ^ HVP-LO(√(s) < 1.05   GeV) and estimate
            additional systematics possibly due to modeling effects and to observed tension among some data samples.
            Complementing this piece by the non-HLS part of the HVP derived by more usual means, one
            gives our best evaluation of the full HVP and our estimate of the muon anomalous
             magnetic moment.
            We also discuss an issue in the relationship between the LO-HVP
            results derived by dispersive approaches – in particular BHLS/BHLS_2– and the
            LQCD results, where the ρ^0-γ mixing appears as a NLO effect.
            Comparison is done with other published evaluations for a_μ^ HVP-LO.
            Finally, Section <ref> summarizes our results and conclusions.
            
            
            
            
            § OUTLINE OF THE HIDDEN LOCAL SYMMETRY (HLS) MODEL
             Let us briefly outline  the derivation of the HLS
            model[See, for instance, <cit.> for a full derivation.].
             The usual ChPT Lagrangian <cit.> can be written
            in two different manners <cit.> :
             L_ chiral=f_π^2/4 Tr [ _μ U ^μ U^† ]
            =-f_π^2/4 Tr [ _μξ_L  ξ_L^† - _μξ_R  ξ_R^† ]^2 ,
            
            where f_π (=92.42 MeV) is the pion decay constant and :
            [                            U(x)= exp [ 2iP(x)/f_π]  , ξ_R/L(x)=exp [ ± iP(x)/f_π]  ⟹
            U(x)=ξ_L^†(x) ξ_R(x) , ]
            when working in the so-called unitary gauge  which removes a scalar field term  in the definition of
            ξ_R/L(x); P(x) is the  pseudoscalar (PS) field matrix.
            The hidden local symmetry which gives its name
            to the HLS Lagrangian has transformation properties which affect   ξ_R/L(x)   letting U
            unchanged. Ignoring for the moment the weak sector to ease the discussion,
            the HLS Lagrangian is derived by replacing  in  Equation (<ref>) the usual derivative  by the
            covariant derivative :
            
             D_μξ_R/L =_μξ_R/L -ig V_μξ_R/L +ie ξ_R/L A_μ Q ,
            
            where A_μ is the photon field, Q= Diag[2/3,-1/3,-1/3]  the quark charge matrix and
            V_μ is the  vector field matrix; the expressions for P and V are the usual
            ones – fulfilling the U(3) flavor symmetry –
            and can be found in <cit.>. In the expressions for D_μξ_R/L,
            the universal vector  coupling g
            occurs beside the unit electric charge e. The  and ϕ fields occurring in the
            diagonal of the V matrix, also denoted V^I below, are the so-called ideal combinations
            generally denoted _I and ϕ_I.
            
            Substituting the covariant derivative Equation (<ref>)  to the usual derivative in Equation (<ref>),
             L_ chiral becomes the first HLS Lagrangian piece denoted
            L_A while  another piece L_V shows up which vanishes
            in the reversed substitution  D_μ⇒_μ :
            [ L_A=-f_π^2/4 Tr [ D_μξ_L  ξ_L^† - D_μξ_R  ξ_R^† ]^2  ,    L_V=-f_π^2/4 Tr [ D_μξ_L  ξ_L^† + D_μξ_R  ξ_R^† ]^2 ]
            The full HLS Lagrangian is then defined by :
             L_ HLS= L_A +a  L_V ,
            
            where a is a parameter specific of the HLS approach <cit.>. The usual VMD
            framework is obtained by setting a=2<cit.>  while the presently reported phenomenology
            prefers larger values <cit.>. The explicit expression for this unbroken
            L_ HLS  can be found fully developed in <cit.>. Let us note
            that L_ HLS contains a photon mass term; however this Reference also showed that
            the loop dressing (in the HLS context) cancels out the physical photon mass and thus, this mass term
            can be safely ignored in phenomenological studies.
            
            
            So, the unperturbed ( i.e. unbroken) HLS Lagrangian depends on only two parameters
            (g and a)  which can be adjusted using data;
            it is obviously unrealistic to expect this basic  L_ HLS to describe
            precisely the large amount of data covering its scope with so few parameters. To have any
            chance to
            account for real data, additional  input should enter the HLS Lagrangian  to make it more
             flexible;  this is what motivates the introduction  of breaking procedures within the HLS framework.
            
            As already noted, both HLS Lagrangian pieces fulfill a U(N_f) × U(N_f) symmetry and
            not SU(N_f) × SU(N_f). To reduce this symmetry which introduces a ninth  PS
            meson, one  includes <cit.>  the 't Hooft determinant terms <cit.> which
            break
             the axial U(1) symmetry; this turns out to complement the L_A Lagrangian piece
             with :
             L_ 'tHooft=μ^2/2η_0^2+ λ/2_μη_0^μη_0
            The singlet mass term  manifestly breaks nonet symmetry in the PS sector. Nothing analogous
            affects the vector sector.
            
            
            
            § BREAKING THE HLS LAGRANGIAN I : THE BKY MECHANISM
             Up to now, a single breaking mechanism  for the HLS Lagrangian has been
            proposed <cit.>; it is named  here BKY after its proponent names. However,
            undesirable properties of the original  BKY proposal have led to  modifications
            of this breaking scheme  <cit.> originally proposed
             to break  (only) the flavor symmetry of the Lagrangian.
             The "new scheme" variant from  <cit.> has been analyzed with data
            and its results shown to compare fairly well to ChPT expectations <cit.>.
            It has been extended
            to include isospin breaking effects in <cit.> following the lines of
            <cit.>. This (modified) BKY breaking scheme, intended to cover isospin and
             SU(3) breaking effects within the same framework, has been proved fairly
             successful, allowing for a high quality
             global fit of (almost) all available data covering the validity range  of the HLS model
             <cit.> :  e^+e^- annihilations, τ decay spectra and
             radiative decays of light flavor mesons; in particular, the dipion spectra in
             the η/π^+
             π^- γ decays are fairly well predicted <cit.>, especially the
             detailed  form of the drop-off in the ρ- interference region as recently
             measured at BESIII <cit.>. Nevertheless, as already stated in
             <cit.>, the description of the threshold region in
             the dipion spectrum and of the ϕ mass region  in the three pion annihilation
             channel deserves improvements; this issue is addressed by the present paper.
            
             As BKY is one of the breaking mechanisms  used in the present work,
             let us briefly remind  how it is implemented. To lighten writing, we use the notations
             L=D_μξ_L  ξ_L^† and R=D_μξ_R  ξ_R^†.
            The (modified and extended) BKY breaking which also underlies the broken HLS (BHLS) model
              developed in <cit.> turns out to modify Equations (<ref>) as follows :
             [  L_A=-f_π^2/4 Tr [ (L - R) X_A ]^2
              , L_V=-f_π^2/4 Tr [ (L + R) X_V ]^2
               , ]
            where X_A/V=Diag(q_A/V,y_A/V,z_A/V) are constant real matrices. In practice, one prefers setting
            q_A/V=1+(Σ_A/V+ Δ_A/V)/2 and y_A/V=1+(Σ_A/V- Δ_A/V)/2.
            As z_A/V are affecting the ss entries, their departure from 1 can be large
            compared to q_A/V and y_A/V– which refer to resp. the uu and
            dd entries; so,  the Σ's and Δ's are expected small. Within
            the  previous broken HLS framework – hereafter named BHLS –<cit.>, one got
            z_A≃ [f_K/f_π] ≃ 1.5,  z_V≃ 1.2  while  the Σ's
            and Δ's were found at the few percent level.
            
            
            Once these breakings are applied, the  PS kinetic energy term contained in L_A
            is no longer diagonal and a PS field
            redefinition has to be performed in order to restore the kinetic energy term to canonical
            form. The procedure is fully described in Section 4 of <cit.> and is valid unchanged
            in the present work; it is not  repeated here.
            
            In the BHLS model <cit.>, after the X_V breaking, the vector meson mass term is no
            longer diagonal and writes (m^2=a g^2 f_π^2) :
            
              L_V^ mass= m^2/2 [
            (1+Σ_V) (ρ^2_I +^2_I)+ 2 Δ_V ρ_I ·_I +z_V ϕ_I^2+
            2(1+Σ_V) ρ^+ ρ^-
            ] ,
            
            discarding the unessential K^* sector and neglecting terms of degree greater than 1
            in the breaking parameters.
            This mass term can be diagonalized by a 45^∘ rotation, a quite
            unacceptable solution as it does not give a smooth limit in the symmetry limit Δ_V  0.
            A smooth solution
            is obtained by the following transform   to renormalized (R) fields :
               (
            [ ρ_I;  _I ] )
            =
              (
             [ ρ_R;  _R ] )
            - Δ_V   (
             [      h_V_R; (1-h_V)ρ_R ] ) ,
            
            where h_V is a parameter submitted to fit, together with z_V, Σ_V and Δ_V.
            This solution, adopted in  <cit.>, provides a quite satisfactory
            description of the data.  This transformation introduces a Δ_V ρ_R _R term in
            the kinetic energy term of the vector mesons; this issue is known to imply
            the occurrence of wave-function renormalization factors absorbed in the effective
            couplings <cit.>; in the case of BHLS, they can be considered absorbed in the
            breaking parameters.
            
            The full account of the BHLS model also requires the dynamical breaking <cit.>
            generated by PS loop effects
            which also calls for an additional  s-dependent  renormalization step. It is rediscussed below
            within the context of the present work.
            
            
            
            § BREAKING THE HLS LAGRANGIAN II : THE COVARIANT DERIVATIVE (CD) HANDLE
             As clear in Equations (<ref>),   BKY  breaks the symmetry
            of the HLS Model  in the very definition of its Lagrangian. This is
            quite legitimate  and strongly validated by its remarkable
            description of a considerable amount of data <cit.>.
            Nevertheless, it is of interest to explore other possibilities  to improve
            the description of the data in some specific mass regions and, then, reduce
            or cancel out  the additional systematic uncertainties reported in <cit.>.
            
            
            Let us focus on  Equation (<ref>) which defines the original HLS covariant derivative :
             
                D_μξ_R/L =_μξ_R/L -ig V_μ^I ξ_R/L +ie ξ_R/L A_μ Q   .
            
            The superscript  I assigned to the vector field matrix V aims at reminding that the Isospin zero mesons
            are the so-called ideal combinations ^I and ϕ^I. These naturally occur in the U(3) symmetric
            expression of V.
            
            
            
            The vector meson term in D_μξ_R/L
            can be written g V^I= g ∑_a=0,8 V_a T_a where the V_a's are the vector meson
            fields and  T_a,  (a=1,⋯ 8) are the Gell-Mann
            matrices normalized such that Tr[T_a T_b]=δ_ab/2;   T_0 is  the
            unit matrix appropriately normalized : T_0=I/√(6). Thus, the (nonet)  U(3)_V  symmetry
            of HLS  is obtained by :
            
              * 
            Plugging the appropriately constructed V^I into the covariant derivative
            D_μξ_R/L,
            
              * 
            Assuming the universality of the vector coupling g to the external world.
            
            Therefore, we propose a direct breaking of the covariant derivative, a new tool
            independent of the BKY mechanism.
            Such a breaking mechanism is expressed by the modified covariant derivative :
             
            D_μξ_R/L =_μξ_R/L -ig [V_μ^I +δ V_μ ] ξ_R/L +ie ξ_R/L A_μ Q ,
            
            where δ V_μ can be chosen to break the U(3)_V symmetry in a controlled way.
            For instance,  breaking solely the nonet symmetry  of  V^I turns out to allow
            the coupling to the  singlet V^0_μ to differ from those of the octet fields, preserving
            in this way the SU(3)_V symmetry. It will become clearer below why  a more systematic approach
            must be preferred.
            
            Identifying the field combinations associated with each of the canonical T_a matrices,
            one is led to define the following components which can participate
            to δ V_μ separately or together :
            {[      δ V_μ^0 =ξ_0/√(2)[ √(2)_μ^I + ϕ_μ^I/3] Diag [1,1,1]
             ,; δ V_μ^8 =ξ_8/√(2)[ _μ^I - √(2)ϕ_μ^I/3√(2)] Diag [1,1,-2]  ,;               δ V_μ^3 =ξ_3/√(2)[ρ_I^0/√(2)] Diag [1,-1,0] , ] .
            
            in terms of the usual ideal field combinations.
            
            The (free) breaking parameters ξ_0, ξ_8 and  ξ_3 are  only requested to be real
            in order that δ V_μ is hermitian as V_μ^I itself. Clearly, δ V_μ^0 defines
            a breaking of the nonet symmetry down to SU(3)_V × U(1)_V,  δ V_μ^8 rather
            expresses the breaking of the  SU(3)_V symmetry, while δ V_μ^3 is related to
            a direct  breaking of Isospin symmetry in the vector sector.
            
            
            One could also introduce some breaking affecting the  ρ^± entries of the V^I
            matrix[For the sake of simplicity,
            the case for  the K^* sectors is left aside; substantially,  this would introduce two additional
            real parameters ξ_4,5 and ξ_6,7 which may be fixed by the two radiative K^* decays;
            additional constraints may come from the τ decay to Kπν.
            ].
            The  ρ^± entries are given
            by the sum of the Gell-Mann matrix terms  T_1 and T_2; forcing a breaking
            for these entries requires two real parameters which should be equal
             (ξ_1,2) in order to preserve hermiticity.
            
            So,  ξ_3 and this ξ_1,2 could
            summarize the whole isospin breaking effects in the vector meson side;
            however, one can choose ξ_1,2=0 as
             no theoretical value for the (unbroken) universal coupling g is presently available.
            Indeed, all vector couplings in the HLS Lagrangian
            could then be reexpressed in terms of  g^'=g(1+ξ_1,2); in this case,  all physics
            quantities will depend on g^' and ξ_i^' =ξ_i/(1+ξ_1,2) (i=0,8,3) without any
            other occurrence of ξ_1,2 dependency.
            
            Therefore,
            phenomenologically, the ξ_i's and g itself are defined up to  a normalization factor
             presently out of reach. This scaling
            property has obviously no consequence on physics observables like cross-sections or form factors.
            
            So, from now on, one assumes the maximal breaking of U(3)_V experimentally accessible :
            δ V_μ=δ V_μ^0+δ V_μ^8+ δ V_μ^3. Compared with the original V^I entries,
            this turns out to modify
            only the diagonal entries of V^I by the following substitutions[
            Remind that a  1/√(2) term is factored out in the definition of V^I.] :
            {[                                          ρ_I/√(2)⇒                                ρ_I/√(2) [1+ξ_3]  ,;                                           _I/√(2)⇒ _I/√(2) [
             1+2 ξ_0 +ξ_8/3 ]  +  [ξ_0 -ξ_8]/3 ϕ_I ,;                                              ϕ_I ⇒    ϕ_I  [
             1+ξ_0 +2 ξ_8/3 ]  + √(2)[ξ_0 -ξ_8]/3 _I ] .
            
            Then, the U(3)_V breaking of the covariant derivative generates a breaking of
            the vector coupling universality. For this purpose, one should note that a vector
            mixing is generated,  except  if ξ_0=ξ_8.
            For the sake of conciseness, this mechanism is referred to below as CD breaking.
            
            
            § BHLS_2 : A NEW BROKEN HLS LAGRANGIAN
              Here, we define a new version of the broken HLS Lagrangian, merging
            the two breaking schemes just presented – BKY and CD – as
            these two ways of breaking are conceptually  independent.
            Their contributions are thus complementing each other.
            In order to fully take into account the electroweak sector,
            one should modify correspondingly the covariant derivative to :
             {[ D_μξ_L =_μξ_L -ig [V_μ^I +δ V_μ
            ] ξ_L +i  ξ_L L_μ ,; D_μξ_R =_μξ_R -ig [V_μ^I +δ V_μ ] ξ_R +i  ξ_R R_μ , ] .
            
            where, as usual<cit.> :
             [ L_μ =     e Q A_μ
            +g_2/√(2) (W^+_μ T_+ + W^-_μ T_-)  and                                     R_μ =    e Q A_μ   , ]
            discarding the Z^0 boson terms of no concern for the phenomenology we address. The T_± matrices are constructed
            out of the matrix elements  V_ud and V_us  of the Cabibbo-Kobayashi-Maskawa matrix and can be found
            in <cit.>. Finally, the weak coupling g_2 is related to the Fermi constant by
            g_2=2 m_W √(G_F√(2)).
            
            The expression for the non-anomalous HLS Lagrangian pieces given in Section <ref> remains formally valid,
            being understood that the covariant derivatives are modified according
            to Equations (<ref>–<ref>).
            As a result, the L_V  piece substantially differs from its partner in <cit.> while
            the L_A piece in the present scheme is strictly identical to those given in
            <cit.>. The non-anomalous BHLS_2 Lagrangian is given expanded in Appendix <ref>.
            
            
            
            § THE ORDER O(P^4)  TERMS OF THE HLS LAGRANGIAN
             Beside the L_A and L_V pieces which are  O (p^2),
            the HLS approach also possesses terms of order O (p^4)(see Section 4.3 in <cit.>),
            which modify the V-γ/W couplings in a specific way. As the role of such terms has never been
            really examined in phenomenology[See, however,  the discussion in <cit.> which
            is revisited here.],
            it looks worthwhile examining their relevance when dealing with data of high
            accuracy.
            
             Using L_μ and R_μ  just given,
            one first defines :
            {[                           L_μ,ν= _μ L_ν -_ν L_μ
            -i [ L_μ, L_ν ] ,;                           R_μ,ν= _μ R_ν -_ν R_μ
            -i [ R_μ, R_ν ] , ] .
            
            and also :
            [ L_μ,ν=ξ_L  L_μ,νξ_L^†  , R_μ,ν=ξ_R  R_μ,νξ_R^†  , ]
            where ξ_R=ξ_L^†=exp[iP/f_π]. Furthermore defining :
            [ V_μ,ν=1/2 [  R_μ,ν +   L_μ,ν]
               ,   A_μ,ν=1/2 [  R_μ,ν -  L_μ,ν]   , ]
            the O (p^4) Lagrangian writes <cit.> :
             L_z=
            z_1  Tr [  V_μ,ν V^μ,ν ]+
            z_2  Tr [  A_μ,ν A^μ,ν ]+
            z_3  Tr [   V_μ,ν V^μ,ν ] ,
            
            where  one has generically defined :
            
            X_μ , ν=_μ X_ν - _ν X_μ -i [X_μ,X_ν] ,
            
            and where the z_i are constants  not theoretically constrained. 
             The z_3 term, the most involved in the phenomenology one addresses, writes :
             L_z_3= ge/2 z_3  A_μ,ν [ ρ^0_μ,ν+1/3_μ,ν-√(2)/3ϕ_μ,ν ] + gg_2/4 z_3  [
            V_ud W^-_μ,νρ^+_μ,ν +V_ud W^+_μ,νρ^-_μ,ν ] 
            
            at lowest order in the expansion of the ξ_L/R fields; here one has kept the unbroken
            V_μ matrix for clarity.
            By integrating by part and fixing the gauge condition to
            ∂_μ X^μ=0 for all vector fields, this piece becomes :
             L_z_3=e g z_3  s  A_μ [ ρ^0_μ+1/3_μ
            -√(2)/3ϕ_μ ]
            +g_2 g z_3 s/2 [V_ud ρ^+_μ W^-_μ +
            	V_ud ρ^-_μ W^+ _μ ]
            Thus, the L_z_3 piece exhibits  s-dependent parts for the
             V-γ and V-W transition amplitudes which complement their  constant
            parts provided by the usual O (p^2) HLS Lagrangian. The  'broken' version
            of  Equation (<ref>) is given in  Appendix <ref>.
            
            
            § THE BASIC SOLUTION (BS) AND THE VECTOR MASS TERM
             In the context where  both the BKY  and CD breaking mechanisms are at work,
            the vector mass term in L_V becomes :
              L_ mass=  a f_π^2  Tr [ X_V (V_μ^I+δ V_μ)  ]^2
            Ignoring the K^* sector and setting  m^2=a g^2 f_π^2, it writes :
            [                                                                                                                  L_ mass= m^2/2 [(1+ Σ_V + 2 ξ_3)  [ρ^0_I]^2 +
                   (1+ Σ_V+ 4/3ξ_0 + 2/3ξ_8)  [_I]^2
                  + z_V(1+2/3ξ_0 +  4 /3ξ_8 ) [ϕ_I]^2 .;                                                                                                                                                               . 
            + 2Δ_V ρ^0_I·_I +
            2√(2)/3   (1+z_V) (ξ_0 - ξ_8) _I·ϕ_I ]
             + m^2  (1+ Σ_V) ρ^+ ·ρ^- ]
            
            Therefore the CD breaking, via its δ V_3 component,  allows for a
             ρ^0-ρ^± (HK) mass difference, provided fits favor
            non-vanishing ξ_3 values. Moreover, all the components
            of δ V defined above  contribute to generate different HK masses
            for the ρ^0 and  mesons as
            Δ m^2_ρ^0-=m^2[ξ_3-(2 ξ_0+ ξ_8)/3]. This contrasts with the
            previous BHLS model <cit.> where m^2_ρ^0=m^2_ρ^±=m^2_
            is fulfilled.
            Because of the non-vanishing HK mass difference between the ρ^0 and ,
            one can already expect an improved treatment  of the dipion threshold
            and spacelike regions in the BHLS_2 framework compared to BHLS.
            
            Equation (<ref>) is manifestly diagonalized by setting :
            [    Δ_V=0   , ξ_0 = ξ_8  . ]
            This solution – which lets z_V unconstrained– defines  our
            Basic Solution (BS). Within
            the framework of this solution,  the breaking parameters to be derived from
            fits are Σ_V,  z_V,  ξ_3, ξ_0 (=ξ_8). Here one feels the issue with
            assuming solely nonet symmetry breaking; indeed, as this turns out to state  ξ_8=0,
            the CD breaking implies ξ_0 =0 and thus breaking HLS via the covariant
            derivative intrinsically implies that nonet symmetry cannot be solely broken.
            
            
            
            
            At leading order in breaking parameters, the vector meson mass term
            in L_V becomes diagonal and one has :
            {[                                              m_ρ^±^2=  m^2  (1+Σ_V)
             ,;                                   m_ρ^0^2=  m^2  [
            1+Σ_V + 2  ξ_3 ]
             ,; m_^2=  m^2  [
            1+Σ_V + 4/3 ξ_0 +  2 /3 ξ_8 ]
            =m^2  [
            1+Σ_V +2  ξ_0 ]
             ,;   m_ϕ^2  =
             m^2  z_V  [ 1+2/3 ξ_0 + 4/3 ξ_8 ]
            =  m^2  z_V  [ 1+2 ξ_0 ];                                                                        ] .
            
            
            So, BHLS_2 a priori yields different HK masses for all the vector mesons.
            
            
            
            § THE PRIMORDIAL MIXING (PM) AND THE REFERENCE SOLUTION (RS)
            
            Another unused mechanism can be invoked; indeed, independently of the BKY and CD mechanisms
            just defined,  one can always
            consider that the neutral vector fields ρ^0, , ϕ involved in physical processes
            are not
            directly the ideal ones but combinations of these. For this purpose,
            let us define an infinitesimal rotation matrix R(U_3) = 1+ O(ϵ), which associates to the
            ideal field vector V_I=(ρ_I, _I,ϕ_I)  a (first step) renormalized field vector
             V_R=(ρ_R,_R,ϕ_R) via :
              (
            [ ρ_I;  _I; ϕ_I;     ] )
             =  (
            [    1  -ψ_  ψ_ϕ;   ψ_    1  ψ_0; -ψ_ϕ -ψ_0    1 ] )
              (
            [ ρ_R;  _R; ϕ_R;     ] )   ,    also  denoted    V_I=R(U_3) · V_R .
            
            This is also a quite legitimate tool to extend the flexibility of BHLS_2.
            
            The matrix R(U_3)  fulfills  R(U_3)R(U_3)=R(U_3) R(U_3)=1, up to
            terms of order O(ϵ^2) which are discarded in our O(ϵ)
            approach.
            As it is a rotation, this transformation preserves the canonical structure of the
            vector field kinetic term provided by the Yang-Mills Lagrangian
            (up to O(ϵ^2) terms).
            
            However, the (real) ψ (Euler) angles from transformation in Equation (<ref>)
            should be chosen in such a way that L_ mass( V_I) remains canonical
            in the change of fields V_I → V_R,  i.e.
             the crossed terms in L_ mass( V_R) should be canceled out.
             Restarting from the mass term given by Equation (<ref>),  three conditions should
             be fulfilled :
            {[                                               ρ_R-_R =0 ⇒               Δ_V  + [
            4/3{ξ_0+1/2ξ_8 -3/2ξ_3
            }ψ_]  =0  ,;                                             ρ_R-ϕ_R  =0 ⇒  [1-z_V] ψ_ϕ
            +√(2)/3{1+z_V}(ξ_0 - ξ_8)
             +[ ψ_ ψ_0
            ]  =0 ,;                                              _R-ϕ_R  =0 ⇒ 3 [1-z_V] ψ_0 + [1+z_V]  √(2)  (ξ_0 - ξ_8) -3 [ ψ_ ψ_ϕ]=0 ] .
            
            The last bracketed terms in these expressions can be discarded as clearly of second order in
             the breaking parameters {ξ_i, ψ_j}; this already implies that Δ_V=0
             and that ψ_ becomes unconstrained.
            
            As for the parameter z_V– generated by the BKY breaking in its s s entry  –
            the situation deserves  further comments :
            
              *  Eitherz_V is found such
            that[ One may think that a reference magnitude for any generic ϵ
            is O(ϵ) ≃ e=√(4πα_ em )≃ 0.3.]1-z_V≃ O(ϵ), then none of the ψ
            angles is constrained at order O(ϵ)
            and, moreover, the last two Equations (<ref>) then imposes ξ_0=ξ_8.
            
              *  Orz_V is returned by fits  such that 1-z_V is not O(ϵ). Then,
            at first order in breakings, the last two equations imply ψ_ϕ=ψ_0 and
            that they are proportional to ξ_0 - ξ_8.
            
            
            Preliminary fits, performed using this mechanism – named from now on 
            Primordial Mixing (PM) –
            indicate that 1-z_V is  in the range (1 ÷ 2) ϵ and thus,
            in this case also ξ_0 =ξ_8 can be imposed.
            Therefore, besides the Basic Solution (BS), one gets an additional one, we name it
            Reference Solution (RS), which also includes  the Primordial Mixing (PM). Also considering
            the former BHLS model <cit.>, one thus has
            at disposal three different models.
            This allows for a better insight into the systematics and model dependence effects.
            
            
            
            § VECTOR MESON PROPAGATORS AND FORM FACTORS
            
            Let us anticipate  on form factor calculations using the Lagrangian given in Appendix <ref>
            and have a look at the form factor values at s=0– before introducing mixing effects due to
            loop corrections. For this purpose, let us start with a preliminary digression on vector
            meson propagators, especially those for the  and ϕ mesons, as these play an important role
            at the chiral point and in the physics of the (close) spacelike region which is also addressed in the present
            work.
            
            For the ρ meson,  following the pioneering work <cit.>,
            the inverse propagator at one loop can be written (also see <cit.>) :
            
                D_ρ(s)=s-m_ρ^2-Π(s) ,
            
            where the self-energy Π(s) is the sum of the
            loops allowed by the non-anomalous Effective Lagrangian given in Appendix <ref>; in our
            context,
            these are essentially pion and
            kaon loops[ Loops generated by the anomalous HLS Lagrangian pieces <cit.>
            also contribute but can be discarded <cit.> as reminded in the next Section;  however,
            some anomalous loops, suppressed by  a factor of e^2, can have to
            be kept as they may play some role in specific
            expressions as will be
            emphasized below. ].
            Π(s) is a real analytic function
            (Π(s^*)=[Π(s)]^*) which vanishes at s=0 because of current conservation, and is real for
            negative s– in fact, this property holds already below m_π^0^2, the lowest energy hadronic threshold. Hence,
            this implies that D_ρ(0)=-m_ρ^2, where m_ρ^2 is displayed in Equations (<ref>).
            
            In principle, this also applies to the  and ϕ propagators; however, taking into account
            their narrow character, it looks unmotivated to enter into such complications when fitting
            objects like e^+e^- annihilation cross sections. In this case, phenomenology has lengthily
            experienced a broad success using  Breit-Wigner (BW) lineshapes. However, as discussed
            in <cit.>, some physics quantities like the contribution of the e^+e^- π^+π^-
            cross-section to the estimate for the muon HVP, may require some care
            about the behavior of approximations close to the chiral point,
             as this region gives an enhanced contribution  to the muon HVP. For this purpose, <cit.>
             proposed the following modified BW lineshape for the  and ϕ mesons :
             D_V(s) =
             s-m_V^2 -s/m_V^2 [ m_V^2 -m_V^2 -im_V Γ_V
              ] ,
              
            where (m_V,Γ_V) are parameters to be fitted and the m_V^2's
            are the  relevant HK square masses as displayed in Equations (<ref>). Numerically, this
            expression for the D_V(s)'s is
            close to usual BW lineshapes[Setting m_V=m_V in this approximation
            gives D_V(s)=s-m_V^2+i s Γ_V/m_V; more common choices for the imaginary
            part here are m_V Γ_V or √(s)Γ_V, quite analogous
            to our own choice.].
            This BW-like parametrization imposes  the  and ϕ meson inverse propagators to
            fulfill  D_V(0)=-m^2_V  as already does  the ρ meson inverse propagator;
            this property of
            the propagators  at s=0  is essential to recover the expected
            values for the pion and kaon form factors at the origin  within the HLS framework.
            
            On the other hand,  if one wishes to examine the analytic continuation of the  and ϕ propagators
            somewhat inside the spacelike region, it might be desirable to make it real there; this can
            be achieved by simply replacing Γ_V in Equation (<ref>) by
            θ(s) Γ_V– or rather, θ(s-m_π^0^2) Γ_V.
            
            
            
            
             Using the Basic Solution (BS) or the Reference Solution (RS) to redefine
             the physical vector meson fields, one can show, using the Lagrangian given in
             the Appendix  <ref>, that :
            
             [   F_π(0)=1  , F_K^±(0)=1  ,  F_K^0(0)=0 , ]
            up to terms of order O(ϵ^2). One can check that the conditions
            Δ_V=0 and ξ_0=ξ_8 are essential in the derivation of these constraints.
            
             So,  the breaking parameters to be derived from fits are
            Σ_V, ξ_3, ξ_0 (=ξ_8) and z_V when working within the
            BS framework; one should also include the ψ rotation angles
            when extending  BHLS_2 to the Reference Solution.
            
            It is worthwhile reminding that the effects of δ V_μ^0 and δ V_μ^8
             are intimately intricated within our new Lagrangian frameworks.
             Here, indeed, a breaking of solely nonet symmetry ( i.e.ξ_0  0) cannot
             occur if not accompanied  by  a breaking of the SU(3)_V symmetry of
             the same intensity.
            
            In conclusion, the present model (BHLS_2) is not a  trivial variant
            of BHLS <cit.>.  Some important properties of BHLS_2 versus BHLS
            mentioned below, will further substantiate this statement.
            
            
            
            § DYNAMICAL MIXING OF VECTOR MESONS
            
            
             §.§ The Squared Mass Matrix at One Loop …
             As previously noted <cit.> and reminded above,
            all variants of the HLS Model, including BHLS_2 (see Appendix <ref>),  exhibit
            ρ^0//ϕ K K couplings. This implies that, at one loop, the
            squared mass matrix of the  ρ^0//ϕ system receives non-diagonal
            entries,  i.e. the vector fields occurring in the Lagrangian Eq (<ref>)
            are no longer mass eigenstates.  Mass eigenstates are constructed using perturbative
             methods as performed in <cit.>. As the loops are (real analytic)
            functions of  s, the mass eigenstates become also s-dependent. Within the BHLS and 
            BHLS_2 frameworks, the kaon loops produce a s-dependent difference between the 
            ρ^± - W^± and ρ^0 - γ transition amplitudes which has provided
             the first solution <cit.> to the long standing  e^+e^- versus τ puzzle
            <cit.>. Nevertheless, in another Effective Lagrangian context,
            this s-dependent difference can be successfully generated by other means <cit.>.
            
            Basically, the Dynamical (or Loop) Mixing of Vector Mesons has
            been first defined[See also  <cit.>.]
            and motivated in <cit.>. In order to ease the reading,  we remind it
            and emphasize the new features provided by the BHLS_2 context.
            
            The (squared) mass matrix of the ρ^0//ϕ sector can be written :
             M^2(s)=M_0^2(s)+δ M^2(s) ,
            
            where :
             M_0^2(s)= Diag (m_ρ^0^2+ Π_ππ(s),m_^2,m_ϕ^2)   
            The Higgs-Kibble masses m_V are displayed in Equations (<ref>), and Π_ππ(s)
            is the pion loop weighted by the square of the ρ^0 π^+π^- coupling constant.
            Because all the loop functions are real analytic
            function of s, M_0^2(s), δ M^2(s), and hence M^2(s), are hermitian analytic
            matrices (e.g. fulfilling [X(s)]^†=X(s^*)).
            
            In order to construct explicitly δ M^2(s), it is worth reexpressing  some VPP
            coupling constants in a suitable manner. One can write :
            {[ g^ρ_K^+ K^-= G (1-Δ_A/2) g^±_ρ  ,   g^_K^+ K^-= G (1-Δ_A/2) g^±_  ,  g^ϕ_K^+ K^-= G (1-Δ_A/2) g^±_ϕ ,; g^ρ_K^0 K^0= G (1+Δ_A/2) g^0_ρ  ,   g^_K^0 K^0= G (1+Δ_A/2) g^0_  ,  g^ϕ_K^0 K^0= G
            (1+Δ_A/2) g^0_ϕ , ] .
            
            having defined G=ag/(4 z_A). Using Equations (<ref>),  common to both
            the BS and RS variants, one finds :
            {[  g^±_ρ =
            {1 +Σ_V+ξ_3+ ψ_ +√(2) z_V ψ_ϕ} ,; g^±_ =
            {1 +Σ_V+ξ_8 - ψ_ +√(2) z_V ψ_0 } ,;   g^±_ϕ={-√(2) z_V (1+ξ_8)
            + ψ_ϕ+ ψ_0 } , ] .
            
            and :
            {[  g^0_ρ =  {1 +Σ_V+ξ_3- ψ_-√(2) z_V ψ_ϕ} ,; g^0_ =-
            {1 +Σ_V+ξ_8 + ψ_+√(2) z_V ψ_0 } ,;   g^0_ϕ=  {√(2) z_V (1+ξ_8) + ψ_ϕ - ψ_0 } ] .
            
            When dealing specifically with the BS variant, these two sets of equations are used
            by simply dropping out the ψ_α parameters.
            
            Denoting by resp. Π_±(s) and Π_0(s)
             the  amputated charged and neutral kaon loops,
             the V_R^i → V_R^j
            transition amplitudes (i,j=ρ^0, , ϕ) are given by :
            ε_i, j(s)=g^i_K^+ K^-g^j_K^+ K^-Π_±(s)+
             g^i_K^0 K^0g^j_K^0 K^0Π_0(s)
            
            using the notations just defined.
            Then, the elements of the δ M^2(s) matrix are :
            
            δ M^2_i, j(s)=ε_i, j(s)
            
            As in the former BHLS <cit.>, δ M^2(s) is non-diagonal
            for non-zero values of s. Therefore, at one-loop order, the field
            combinations defined in both the BS and RS variants are not mass eigenstates,
            as in the previous BHLS release. This calls for a mass-dependent diagonalization which
            is reminded in the next Subsection.
            
            Let us make a few more remarks about the non-diagonal entries of
             δ M^2.  In the no breaking limit, where
            the neutral and charged kaon loops coincide, the ρ^0-
            and ρ^0-ϕ entries identically vanish; however, the -ϕ
            entry,  proportional to the sum of the neutral and charged kaon loops, does
            not vanish, indicating that the -ϕ mixing is always at work
            at one loop order within the HLS framework.
            
            In the former BHLS, however, if no breaking is applied in the L_A
            sector and if there is no mass breaking among the kaons, there was
            no ρ^0- and ρ^0-ϕ transitions. This remains true  within the
            BS framework defined above. However, within the RS framework, the ψ_α's
            come in such a way that none of the entries of δ M^2(s) vanishes.
            
            
            
             §.§ … And its Diagonalization
             The diagonalization of L_ mass
            at one loop order is performed by means
            of another (s-dependent) rotation matrix to a second step of vector field
            renormalization  R ⇒ R^' :
               (
            [ ρ_R;  _R; ϕ_R;     ] )
             =  (
            [     1 -α(s)  β(s);  α(s)     1  γ(s); -β(s) -γ(s)     1 ] )
              (
            [ ρ_R^';  _R^'; ϕ_R^';       ] )   ,    also  denoted    V_R=R(Loop) · V_R^'
            so that the full renormalization procedure is obtained by :
               V_I ⇒[R(U_3) · R(Loop) ]  · V_R^' .
            
            The R^' fields are the physical vector meson fields. When working within the
            Basic Solution variant,
            V_R≡ V_I and R(U_3) can be dropped out from Equation (<ref>);
            otherwise, within the RS framework,  R(U_3) is given in Equation (<ref>).
            
            The complex "angles" in Equation (<ref>) can be derived from the
            δ M^2(s) matrix elements :
              
            [    α(s)= ε_ρ^0 (s)/λ_ρ(s) -λ_(s)  , β(s)= -ε_ρ^0 ϕ(s)/λ_ρ(s) -λ_ϕ(s)  ,        γ(s)= -ε_ϕ(s)/λ_(s)-λ_ϕ(s) , ]
            
            where the numerators manifestly depend on the kaon loops only
            and the λ's are the eigenvalues of M^2(s) matrix.
            At leading order in breaking parameters, these write <cit.> :
              
            [ λ_ρ(s)=m_ρ^0^2 + Π_ππ(s)+ε_ρ^0 ρ^0(s),                      λ_(s)=m_^2+ε_(s),                  λ_ϕ(s)=m_ϕ^2+ε_ϕϕ(s), ]
            the vector meson masses being those displayed in Equations (<ref>). As the
             three m_V^2's  occurring here are different, the three
             angles defined above vanish at s=0; in contrast,  within
            the former BHLS context <cit.>,  because the HK masses for ρ^0
            and  were equal, α(0) has a non-zero
            limit at  s=0, compromising a proper analytic continuation of the form
            factors downwards into the spacelike region[See, in particular,
            the comments in Section 6.3 of <cit.>.].
            
            In summary the full transform from ideal to fully renormalized ρ_0,
            , ϕ fields is given either by (BS solution) :
             R(I ⇒ R^')  ≡  R(Loop)
            
            or by (RS solution) :
             R(I ⇒ R^')  ≡[R(U_3) · R(Loop) ] =
               (
            [            1   -[ψ_+α(s)]  [ ψ_ϕ+β(s)];    [ψ_+α(s)]            1  [ψ_0 +γ(s)]; -[ ψ_ϕ+β(s)] -[ψ_0 +γ(s)]            1 ] ) ,
            
            where breaking terms of order greater than 1 have been discarded. We will refer from now on to
            Equation (<ref>) for both the BS and RS variants, being understood that in the former case, the
            ψ_i's are zero. When relevant, the entries in the  matrix just above will be denoted
            α(s), β(s) and γ((s), using obvious notations.
            Moreover, in the following, the fully
            renormalized fields can be either  indexed by R^'– if useful – or deprived of
            indexation to lighten expressions.
            For  clarity, in all displayed Lagrangian pieces,
            one always writes the couplings for the  idealρ_0, , ϕ fields.
             To go to physical fields, one has to apply
            the appropriate transformation R(I ⇒ R^') and collects the
            contributions in order to yield
            the  physicalρ_0, , ϕ couplings.
            
            
            
             §.§ The Effects of Anomalous Loops
             When defining the loops contributing to δ M^2(s), we have only considered
            those generated by the non-anomalous HLS Lagrangian. However, the anomalous (FKTUY)
            HLS Lagrangian pieces
            <cit.> allow for (anomalous) couplings also generating  loop
            contributions[The Yang-Mills part of the full Lagrangian
            may also contribute with other loops like  K^* K^*.
            The VPPP Lagrangian  generate 2-loop contributions with coefficients of O(ϵ)
            or O(ϵ^2) to the transitions among the vector mesons; they are presently ignored.
            ]
            to δ M^2(s) as, for instance, K^* K or 3-pion loops.
            As all the two particle channels have thresholds above our fitting region, their loops –
            of order O (g^2)– are real  in our
            fitting range and, thus, can be discarded, assuming  their  effects
            numerically absorbed by the subtraction polynomials of the other
            (pion and kaon) loops involved <cit.>.
            
            However, one should also note that the anomalous (FKTUY) sector of the HLS Lagrangian
            generates couplings of the neutral vector mesons to the π^0 γ, ηγ
            and  γ  final states, even
            when setting c_3 = c_4<cit.>. Then the FKTUY sector generates the corresponding loops, which have
            couplings of order O(g^2 e^2) and develop tiny imaginary parts far inside our fitting region for the first two
            and close to its upper border for the third one.
             These contributions, and their imaginary parts, are higher order  and can generally be discarded;
            nevertheless, their tiny imaginary parts were accounted for in our
            previous studies <cit.>
            by simply adding some fixed i ε to the eigenvalue differences
            λ_i(s)-λ_j(s) which occur in the denominators of the mixing
            "angles" shown in  Equations (<ref>).
            
            So, these (complex) mixing "angles" exhibit a dependence upon differences of the M^2(s)
            eigenvalues. However, these eigenvalues – see  Equations (<ref>) –
            being s-dependent, their differences may exhibit zeros for real values of s
            which, accordingly, generate  real poles for the mixing "angles". Poles occurring at real negative
            s are not a real issue as they can be dropped out by means of customary dispersive
             techniques.
            For the others, the  i ε trick emphasized just above, permits to avoid unwanted
            singularities on the physical region s ≥ m_π^0.
            
            
            In the previous BHLS version<cit.>, the main purpose of this  ad hoci ε  was  to
            force[As clear from the first Equations (<ref>), the numerator of α(s)
            behaves as s close to the origin and, as m_ρ^0=m_ in  the previous BHLS
            <cit.>, λ_ρ(s) -λ_(s) exhibits a similar s-behavior.]α(s=0)=0, which
            is no longer  necessary in the present version of BHLS.
            A more natural way to proceed – giving the same results – is, however,  to take into account
            the π^0 γ loop which introduces a small imaginary part to the eigenvalue differences,
            for s_0 ≥ m_π^0^2
            upwards. For this purpose, we use the loop expression derived in  <cit.> (see  its Appendix C)
             which can be written (s_0=m_π^0^2) :
            ϵ_π^0^V(s)= [g_Vπ^0γ^R]^2/96 π^2 [ G(s)-s_0^2 ] ≡ [g_Vπ^0γ^R]^2 ϵ_π^0(s)
            Along the real  axis in the complex s-plane, one has :
            {[          G(s)=-(s_0-s)^3/sln(s_0-s)/s_0   ,                                      s≤ s_0; G(s)=(s-s_0)^3/s [
            ln(s-s_0)/s_0 -i π ]   ,                                      s≥ s_0 ].
            
            The derivation of the g_Vπ^0γ^R couplings is given in Subsection <ref>.
            The constant term inside the
            bracket in  Equation (<ref>) ensures that ϵ_π^0(0)=0.
            Full consistency would impose to add (using obvious notations) :
            δ_2 M^2_i,j(s)= [ g_V_i π^0γ^R g_V_j π^0γ^R  ] ϵ_π^0(s)
            + [ g_V_i ηγ^R g_V_j ηγ^R  ] ϵ_η(s)
            + [ g_V_i γ^R g_V_j γ^R  ] ϵ_(s)
            
            to δ M^2(s) as expressed above; however, because of the e^2 factor, they are  higher order
            and can be neglected, except in the diagonal where they avoid possible real poles in the physical region. Additionally,
            the ϵ_η(s) and ϵ_(s) contributions displayed for completeness
            will also be discarded.
            
            In order to substantiate the effect of the π^0 γ loops,  let us anticipate on the
            global fit information.
            The main effect of these loops is to prevent naturally the occurrence of poles on the physical
            region.
            This effect is prominent for the  "angle" α(s)  which produces a ρ^0- mixing;
             Figure <ref> displays the behavior of Re[α(s)] for the
             BHLS_2  BS and RS variants
             and also, for comparison, its behavior within  the former BHLS, where a iϵ
             trick prevents a real pole close to the 2-pion threshold <cit.>.
             Even if the overall lineshapes are similar, one observes significant differences between the
             various solutions; for instance, α(s) vanishes at s=0 within the two variants of BHLS2,
             as the denominator is non-zero at the chiral point, allowing for a smooth connection between the
             spacelike and timelike regions.
            
             Over the range  [-3.5 ÷ 1.2] GeV^2, the behavior observed
            for the two other mixing angles  β(s) and γ(s) is much smoother
            (see also <cit.>);
            nevertheless, specific
            parametrizations may  lead to a real pole[One should note that a pole at such
            a location is exhibited by the usual Gounaris-Sakurai formula <cit.>.]
            for β(s) at negative s
            below ≃ -1.5 GeV^2. As noted above, such poles can be eliminated by means of usual dispersive
            methods.
            
            
            § THE MODEL PION AND KAON FORM FACTORS
             Using the non-anomalous sectors of the HLS model, broken
            as pointed out in the previous Sections (see Appendix <ref>),
            one can already derive some of the form
            factors expressions needed for our global fit framework.
            
            
            
             §.§ The Pion Form Factor in the τ Decay
             The τ sector of BHLS_2 is almost identical to that of BHLS
            <cit.>; indeed the main L_τ piece reminded in Appendix <ref>
            is unchanged, but
            the newly introduced  L_z_3  contributes. This yields :
            
            F_π^τ(s) = [
            1-a/2(1+Σ_V) ]- ag/2(1+Σ_V)F_ρ^τ(s)
            1/D_ρ(s) ,
            
            with :
            {[                                 D_ρ(s)=s-m^2_ρ -Π_ρρ^τ(s) ,;                         F_ρ^τ(s) =f_ρ^τ -g z_3 s - Π_W(s) ,; f_ρ^τ =a g  f_π^2 (1+Σ_V)    ,   m^2_ρ=a g^2 f_π^2(1+Σ_V) , ] .
            
            where z_3 is a new HLS parameter <cit.> introduced by the O(p^4) terms of the
            HLS Lagrangian (see Appendix <ref>).  Π_ρρ^τ(s) is the loop correction to
            the ρ^± propagator and Π_W(s) the loop correction to the W^±-ρ^± transition amplitude;
            both  are defined just below.
            
            Using the following short-hand notations :
            [ G_π= ag/2   , G_K=ag/4z_A , ]
            one  first defines the pion and kaon loop contributions to the ρ^± self-energy
             Π_ρρ^τ(s) :
            [ Π_π^τ(s) =
              [ℓ_π^τ(s) + P^τ_π(s) ]   ,   Π_K^τ(s) =
             [ℓ_K^τ(s)  + P^τ_K(s) ] , ]
            where ℓ_π^τ(s) and ℓ_K^τ(s) denote  the amputated loop functions for resp.
            π^±π^0 and K^± K^0, each having absorbed  only a factor of resp. G_π^2 and  G_K^2;
            in this way, only the breaking parameters affecting the V sector of BHLS_2 are factored out
            and displayed.
            P^τ_π(s) and P^τ_K(s) are subtraction polynomials chosen to vanish at s=0 because of
            current conservation.
            Then one has[Here and in the following, one can be led to keep
            breaking parameter dependencies higher than first order to simplify the expression writing.] :
            {[                                Π_ρρ^τ(s)=  [(1+Σ_V)  ]^2 Π_π^τ(s)
            + [ √(2)(1+Σ_V) ]^2Π_K^τ(s)  ,; Π_W(s)= (1+Σ_V)/G_π [ 1-a/2(1+Σ_V) ] Π_π^τ(s) +
            √(2)(1+Σ_V)/ z_A  G_K [z_A-a/2(1+Σ_V) ]
            Π_K^τ(s) ] .
            
            As the loop functions  vanish at s=0, one clearly has F_π^τ(0)=1.
            One should note that the same subtraction polynomials occur in the W-ρ
            transition amplitude  F_ρ^τ(s) and in the inverse propagator D_ρ(s) but
             differently weighted.
            
             The expression for the loop functions ℓ_π^τ(s) and ℓ_K^τ(s)
             along the real s-axis can be found in Subsection F.2 of the Appendix to <cit.>;
             the subtraction polynomial needed to ensure that they vanish at s=0 has coefficients
             given in Equations (122) of this reference. The coefficients for the constant term  and for
             the second degree one given therein are correct, but the coefficient for the term linear in s should
             be corrected to :
             
                c_1= 1/π[ 1 - m_0^2+m_c^2/2 m_0 m_clnm_0-m_c/m_0+m_c]
            
            using the notations[m_0^2 and m_c^2 are the s values
            of the direct and crossed thresholds energies.] of <cit.>.
            
            
            
             §.§ Parametrization of Loops in e^+e^- Annihilations and τ Decays
             The ρ^0 self-energy and the
            ⇒ (ρ^0//ϕ) transition amplitudes involve important
            loop corrections; these already play a central role in the dynamical mixing
            of vector mesons as shown in Section <ref>. Let us define a parametrization
            common to the computation of the "angles" α(s), β(s) and γ(s),
             to the self-energies and to the A ⇒ V transition amplitudes as these are
            involved in all form factors addressed in our global fitting code.
            
            One first defines the following loops  :
            [        Π_π^e(s) =
              [ℓ_π^e(s) + P_π^e(s) ] , Π_K^± ^e(s) =
             [ℓ_K^±^e(s)  + P^e_K^±(s) ] , Π_K^0 ^e(s) =
             [ℓ_K^0^e(s)  + P^e_K^0(s) ] , ]
            where ℓ_π^e(s), ℓ_K^±^e(s) and ℓ_K^0^e(s) are resp. the
            π^+ π^-, K^+ K^- and K^0 K^0 loops. These are defined
            as the amputated loops, including only a factor of G_π^2 for the former loop and
            of G_K^2 for the two others (see Equations (<ref>)). The expression for these loops
            along the real s-axis  can be found in Section F.1 of the Appendix to <cit.>.
            
            The loop parts being defined this way, P_π^e(s), P^e_K^±(s)
            and P^e_K^0(s) are the subtraction polynomials, chosen of the second degree
            and vanishing at s=0.
            
            However, as already observed in Section <ref>, besides the pion loop, it is rather
            the sum and difference of the kaon loops which are relevant :
            {[ Π_K_ sum ^e(s) =
             [
            ℓ_K^±^e(s)+ ℓ_K^0^e(s)+ P^e_S(s)
             ] ≡Π_S ^e(s)  ,;    Π_K_ diff ^e(s) =
             [ℓ_K^±^e(s)-ℓ_K^0^e(s)+P^e_D(s) ]
             ≡Π_D ^e(s) , ] .
            
            where P^e_S(s) and P^e_D(s) are denoting their respective subtraction polynomials.
            As already noted, one can consider Π_K_ diff ^e(s) and P^e_D(s)
            to be  first order in breaking (O(ϵ)). So, it is certainly more appropriate
            to parametrize P^e_S(s) and P^e_D(s) for fits and propagate them into
            the subtraction polynomials of Π_K^± ^e(s) and Π_K^0^e(s) by
            stating :
             [ P^e_K^±(s) =1/2 [  P^e_S(s)+ P^e_D(s)  ]  ,  P^e_K^0(s) =1/2 [  P^e_S(s)- P^e_D(s)  ] , ]
            keeping in mind that
            P^e_D(s)  is O(ϵ), while P^e_S(s) is O(1).
            
            
            
            
            On the other hand, as one certainly has :
             {[      Π_π^τ(s) ≃Π_π ^e(s) +  O(ϵ) ,; Π_K^τ(s) ≃Π_K_ sum ^e(s) +  O(ϵ) , ] .
            
            and it is appropriate to impose :
             {[ P_π^τ(s) =P^e_π(s)+ δ P^τ(s) ,;           P_K^τ(s) =P^e_S(s) , ] .
            
            and submit to fit the coefficients of the δ  P^τ(s) polynomial
            rather than those for P_π^τ(s) directly. This turns out to attribute
            the full breaking in e^+e^- versus τ to the pion loop subtraction
            term; as our fitting region is mostly below the KK thresholds, this  looks
            a safe assumption. The second relation, additionally, is an assumption which allows to reduce
            the number of free parameters in the minimization procedure by two units without
            any degradation of the fit quality.
            
            As a general statement, all our subtraction polynomials  have been
            chosen of second degree and vanishing at s=0. An exception is made, however, by
            choosing the third degree for
            δ P^τ(s) which is found to provide a significantly better description
            of the τ spectra, especially for the Belle dipion spectrum <cit.>.
            
            
            
             §.§ The Pion Form Factor in the e^+e^- Annihilations
            
            Because of the vector meson mixing, deriving the pion form factor in the e^+e^-
            annihilation is slightly more involved than in the τ decay.
             To start, one needs
            to propagate the transformation in  Equation (<ref>) into the Lagrangian in Appendix
            <ref> and collect all contributions to the fully renormalized ρ,
             and ϕ field couplings.
            The derived pion form factor in e^+e^- annihilations F_π^e(s)  writes :
            
            F_π^e(s) =  [ 1-a/2(1+Σ_V )] -
            F_ργ^e(s) g_ρππ/D_ρ(s)
            - F_ωγ^e(s) g_ωππ/D_ω(s)
            - F_ϕγ^e(s) g_ϕππ/D_ϕ(s)
            
            Using the ρ^I π^+ π^- coupling in the Lagrangian  displayed
            in Appendix <ref> and the I ⇒ R transform given
            by Equation (<ref>), one readily gets :
            [  g_ρππ = a g/2(1+Σ_V+ξ_3)  ,   g_ωππ = -a g/2
            [ψ_+α(s)] , g_ϕππ = a g/2  [ψ_ϕ+ β(s)] , ]
            keeping only the leading O(ϵ) terms in breaking. So, BHLS_2 generates
             and ϕ direct couplings to ππ. The inverse ρ^0
            propagator writes :
             
            D_ ρ^0(s) =s-m_ρ_0^2 -Π_ρρ^e(s) ,
            
            where m_ρ_0^2 is defined in Equations (<ref>) and Π_ρρ^e(s)
            is the self-energy of the physical ( i.e. fully renormalized) ρ^0.
            This can be defined in terms of the loop functions constructed in the Subsection just
            above.  Π_ρρ^e(s)  is made up of two pieces :
             {[                                              Π_π^e(s) =   [g^ρ_π ]^2
              [ℓ_π^e(s) + P^e_π(s) ]   ,   g^ρ_π =(1+Σ_V+ξ_3) ,; Π_K^e(s) =
              [ g^ρ_K± ]^2
            ℓ_K^±^e(s) +  [g^ρ_K0 ]^2
            ℓ_K^0^e(s)+
            1/2 [  [g^ρ_K± ]^2+
             [g^ρ_K0 ]^2
             ]  P^e_S(s)+  O(ϵ^2) , ].
            
            where the couplings g^ρ_K± and
            g^ρ_K0 are given in Appendix <ref>.
            A term[  [g^ρ_K± ]^2-
             [g^ρ_K0 ]^2  ] P^e_D(s)/2 has
            been omitted in the expression for Π_K^e(s)  as being O(ϵ^2).
            
             As for the inverse propagators D_ (s) and D_ ϕ(s),
             one uses the expression given by Equation (<ref>)
            with the appropriate  values for HK masses m_V^2, and having
            ([m_V, Γ_V], V=, ϕ) as parameters to be determined
            by fit.
            
            
            
            The loop corrected V-γ transitions amplitudes F_Vγ^e(s) are given by :
             F_Vγ^e(s)= f_Vγ -c_Vγ  z_3  s - Π_Vγ(s)  ,  
            (V=ρ^0_R, _R, ϕ_R)  ,
            
            where the new-comer z_3 s terms occur as in F_ρ^τ(s) above.
            Using together the ideal couplings Equations (<ref>), the renormalization matrix
            Equation (<ref>) and the  RS conditions (see Section <ref>),  one derives :
             {[ f_ργ=a g  f_π^2  [ 1+Σ_V+ξ_3+[ψ_+α(s)]/3
            +z_V√(2)[ψ_ϕ+β(s)]/3 ] ,;  f_γ=a g  f_π^2/3 [ 1+Σ_V+ξ_8-3[ψ_ +α(s)]
            +z_V√(2)[ψ_0+γ(s)]  ] ,; f_ϕγ=a g  f_π^2/3 [-z_V√(2)(1+ξ_8) +
            3[ψ_ϕ+β(s)] +[ψ_0+γ(s)]  ] , ].
            
            which become s-dependent thanks to the dynamical vector mixing. Similarly, using
            the couplings which can be read off  L_z_3  in Equation (<ref>) and
            the transformation matrix given by Equation (<ref>), one also derives :
             {[    c_ργ=g  [ 1+ ξ_3+[ψ_+α(s)]/3
            +√(2)[ψ_ϕ+β(s)]/3 ] ,;     c_γ=g /3 [ 1+ξ_8-3[ψ_ +α(s)]
            +√(2)[ψ_0+γ(s)]  ] ,; c_ϕγ=g/3 [-√(2)(1+ξ_8) +
            3[ψ_ϕ+β(s)] +[ψ_0+γ(s)]  ] , ].
            
            neglecting terms of breaking order greater than 1.
             Finally, the  relevant loop transition terms Π_Vγ(s) are collected in Appendix
            <ref>.
            
            
            
            
             §.§ The Charged and Neutral Kaon Form Factors
              The charged and neutral kaon electromagnetic form-factors can
            also be derived from the BHLS_2 Lagrangian given in Appendix <ref> :
             
            
            [                                                                                     F^e_K_c(s)=                                                             [1 -a/6z_A
            (3+2Σ_V-3Δ_A/2) ][0.5cm]                -F_ργ(s) g_ρ K^+ K^-/D_ρ(s)
            -F_γ(s)g_ K^+ K^-/D_(s)
            -F_ϕγ(s)g_ϕ K^+ K^-/D_ϕ(s) ,;                                                                                     F^e_K_0(s)= -aΣ_V/6z_A
            - F_ργ(s) g_ρ K^0  K^0/D_ρ(s)
            -F_γ(s)g_ K^0  K^0/D_(s)
            -F_ϕγ(s)g_ϕ K^0  K^0/D_ϕ(s) , ]
            
            where the γ - V transition amplitudes F_Vγ and the propagators have
            been already defined and the VKK couplings are displayed in the Appendix
            <ref>.
            
            The experimental spectra  measured in the spacelike region, in particular by the NA7 Collaboration
            <cit.>,  being  squared form factors, the modulus squared of the form factors
            given in Equation (<ref>) and in  Equations (<ref>) apply directly.
            In the timelike region, experimentalists rather publish cross-sections :
            σ(e^+e^-  P P)=8 πα_ em^2/3 s^5/2 q_P^3 |F^e_P(s)|^2 ,
            
            where  q_P=√(s-4 m_P^2)/2 is the meson
            momentum in the center-of-mass system (P=K^+, K^0). In the case of the
            K^+ K^- final state, the cross-section  should take into account the significant
            interaction between the charged
            kaons emerging from the e^+e^- annihilation. Traditionally, this is done
            by multiplying the e^+e^-  K^+ K^- cross-section by the leading term describing
            the Coulomb interaction as formulated by <cit.> :
            
            Z(s)=  [ 1 + πα_ em1+v^2/2 v ]   ,   v=√(s-4 m_K^±^2)/s
            
            As the data have become more and more precise, it looks more convenient
            to use the exponentiated expression derived in <cit.> which includes also
            the Final State Radiation (FSR) correction factor; indeed, the K^+ K^- spectra
            provided by the various experiments are not amputated from FSR effects.
            
            
            
             §.§ Form Factors in BHLS versus BHLS_2
             It looks worth pointing out the difference between the form factor
            expressions in   BHLS and in BHLS_2.   As already noted,
            the L_A Lagrangian is  common to both  frameworks.
            In contrast, the L_V Lagrangian is significantly different.
            Obvious differences for L_V  have
            already been noted.  For instance, within BHLS_2, Δ_V=0 and its sharing parameter h_V
            drop out whereas they are a key ingredient in the diagonalization procedure
            of BHLS.
            
            Three parameters[Additional constraints arising in the diagonalization  procedure of the vector
            mass term, have reduced this number to two as ξ_0=ξ_8 is requested. ]
             arise from the Covariant Derivative breaking (ξ_3, ξ_0, ξ_8).
             Moreover, a key role is played by the ψ rotation angles within the RS solution
             of BHLS_2. None of these parameters occur in the BHLS framework of <cit.>.
            
            The effect of a significant  ξ_3  value is obviously equivalent to
             shifting apart the Higgs-Kibble masses of both ρ's and their couplings
             to ππ in a correlated way. As this kind of exercise performed
             within BHLS concluded to a small effect <cit.>,
             one expects a small value for ξ_3. Also, BHLS_2
             deals with different subtraction polynomials – the δ P^τ(s) function
             already referred to above –
            for the π^±π^0 and π^+ π^- loops, a new (loop) mechanism able alone to
            contribute to physical mass and width differences between the ρ^0 and ρ^±
             mesons.
            
            One should also note another difference : In the present framework,
            as already noted in Subsection <ref>, the pion and kaon loops entering
            the amplitudes F_ρ^τ(s) and F_ρ^e(s) carry fitted subtraction
            polynomials already involved in the charged and neutral ρ inverse propagators
            D_ρ(s). In BHLS, one allowed  the subtraction polynomial in
            F_ρ^e/τ(s)  to carry a piece independent of both
            D_ρ(s)'s. Because of the O(p^4) term still introduced, it
            means that a second-degree term in s has been removed from
            F_ρ^e/τ(s) in the BHLS_2 framework. We should see
            below the effect of this on the ππP-wave phase-shift above the ϕ mass.
            
            
            
            § THE ANOMALOUS SECTOR OF BHLS_2
             As in <cit.>, the anomalous FKTUY
            Lagrangian <cit.>  is used to address
             the e^+e^-  (π^0/η) γ and
            e^+e^- π^0π^+π^- processes, and also the radiative decays  with
            couplings of the form VPγ and P γγ which are parts of
            our global fitting framework. One can display
            its various pieces :
            [                                                 L_ anomalous =
            (1-c_4)  L_AAP+c_3 L_VVP+
             [ 1-3/4(c_1-c_2+c_4)   ] L_APPP;                                                                                +(c_1-c_2-c_3)  L_VPPP+
            (c_4-c_3)  L_AVP , ]
            where, temporarily,  the weight of each piece has been
            extracted out  to exhibit its
            dependence upon the unconstrained  c_i constants of the anomalous HLS Lagrangian
            <cit.>.
            If one imposes to yield the Wess-Zumino-Witten terms <cit.>
            at the chiral point – the so-called triangle and box anomalies –
             the condition c_3=c_4 is known to be
            mandatory <cit.>. Reference <cit.> having also shown that
            this condition is indeed favored by fits to the annihilation data,
            one endorses the c_3=c_4 constraint and, thus, cancels out the last term in
            Equation (<ref>). Hence, the VPγ couplings become  entirely generated
            by the combination of VVP couplings with the  V-γ transitions provided by
             the non-anomalous HLS Lagrangian.
            
            The surviving Lagrangian pieces displayed in Equation (<ref>) are given by :
            {[                                                                              L_VVP=                            - N_c g^2/4 π^2 f_π c_3
             ϵ^μναβ Tr[ ∂_μ V_ν∂_α V_β P]  ,;                                                                              L_AAP=                     - N_c e^2 /4 π^2 f_π (1- c_4)
             ϵ^μναβ∂_μ A_ν∂_α A_β Tr[Q^2 P] ,;                                                                             L_VPPP=           - i N_c g /4 π^2 f_π^3 (c_1-c_2-c_3)
               ϵ^μναβ Tr[V_μ∂_ν P ∂_α P ∂_β P] ,;                                                                             L_APPP= - i N_c e/3 π^2 f_π^3 [1- 3/4(c_1-c_2+c_4)]
              ϵ^μναβ A_μ Tr[ Q ∂_ν P ∂_α P ∂_β P] , ] .
            
            where N_c is the number of colors and the c_i  factors have been  reabsorbed
             in the normalization of the various Lagrangian pieces.
            P and V denote the   bare pseudoscalar  and vector field matrices. In the present
            context, one has, however :
               V_μ= V_μ^I +δ V_μ ,
            
            where V_μ^I is the (usual) U(3) symmetric vector field matrix and
            δ V_μ=δ V_μ^3+δ V_μ^0+δ V_μ^8 which has been
             expressed in term of the ideal field combinations in Equations (<ref>).
            
            For the reader's convenience, Appendix <ref> reminds  briefly the renormalization
            procedure which leads from the bare PS fields to their renormalized partners <cit.>.
            
            
            
             §.§  Renormalization of the VVP Couplings
             As displayed in Equation (<ref>), referring to :
            
                V=(ρ_0, , ϕ)  ,
            
            the relation between the fully renormalized vector fields (denoted R here and from now on)
             and their ideal combination  is given by :
            
                V^I ⇒[R(U_3) · R(Loop) ]  · V^R ,
            
            where the product R(U_3) · R(Loop),   given by Equation (<ref>),  will be denoted from now on :
             R(s) = R(U_3) · R(Loop)
             Let us point out that R(s) is a real analytic 3× 3 matrix function and
              fulfills at  order O(ϵ) :
             R(s+iϵ) · R(s+iϵ)=
             R(s+iϵ)  · R ( s+iϵ)=1
            
            along the right-hand cut on the physical sheet <cit.>. The VVP Lagrangian
            displayed in Appendix <ref> can
            be symbolically written :
            
             L_VVP^ ideal(P)=
            C/2ϵ^μναβ{ P_0  ∂_μ V^I_i ν  G_i,j^I (P_0) ∂_α V_j β^I
             +G_i^I (π^±) [
             π^-  ∂_μρ^+_ν∂_α V_i β^I
             + π^+   ∂_μρ^-_ν∂_α V_i β^I
             ]
             } ,
            
            in terms of ideal field combinations, P_0 being any of the π^0, η or  fields;
            summation over greek and latin indices is understood.
            The 3-vector G^I (π^±) and the 3×3 matrix G^I (P_0) can easily be constructed using
            the relevant pieces of information given in Appendix <ref>. One thus has :
            G^I (π^±)= ( g^±_ρ^0 ρ^±,
            g^±_ρ^±,
            g^±_ϕρ^± ) ,
            
            with the obvious  (ρ, , ϕ) component indexing. On the other hand, one has defined :
            
            G^I (P_0) =
              (
            [    g_ρρ  1/2g_ρ 1/2g_ρϕ;  1/2g_ρ      g_  1/2g_ϕ; 1/2g_ϕρ  1/2g_ϕ    g_ϕϕ ]) ,
            
            where, for each P_0,   g_i j is the coupling of the form P_0  V_i  V_j
            which can be read off the relevant expressions given  in  Appendix <ref>.
            
            Performing  the replacement V^I= R(s)  V^R in  Equation (<ref>),  one can derive the fully
            renormalized VVP Lagrangian :
            
             L_VVP^ ren(P)   =   C/2ϵ^μναβ [
            P_0  ∂_μ V^R_i ν  G_i,j^R (P_0) ∂_α V_j β^R  . [0.2cm]       . + π^-  ∂_μρ^+_νG_i^R (π^±)∂_α V_i β^R
            + π^+   ∂_μρ^-_νG_i^R
            (π^±)∂_α V_i β^R  ] ,
            
            where, the coupling vector G^R (π^±) and the coupling
            matrix G^R(P_0) inherit from R  a s-dependence;
            they are related to their ideal partners by :
            [  G^R (π^±) = R(s) G^I (π^±)   , G^R (P_0)= R(s)G^I
            (P_0) R(s) , ]
            which is a concise way to express the renormalized couplings.
            
            
            
             §.§ The Renormalized AVP Effective Lagrangian
              As already noted,   setting  c_3=c_4 turns out to cancel
            the direct FKTUY AVP Lagrangian piece;  hence, all  AVP couplings inside
            BHLS_2  are
            generated by VVP couplings followed by one  V γ transition.
            Let us consider the V-γ transition term of the non-anomalous BHLS_2 Lagrangian (see
            Equation (<ref>)). Having defined the transposed vector
            f_V γ^I=(f_ργ^I, f_γ^I, f_ϕγ^I),
            one can rewrite it :
             L_AV = -e  [ f_V γ^I· V^I_μ ] A^μ
            =-e  [ f_V γ^R(s) · V^R_μ
            ] A^μ ,
            
            which defines the 3-vector f_V γ^R as :
            
            f_V γ^R=  R(s)f_V γ^I
            
            The effective L_AVP Lagrangian can be derived from Equation (<ref>)
            by replacing one neutral vector meson, say V_i^R, using the rule :
             V^R_i μ⇒ (-ief_γ V_i^R) ( -i/D_V_i(0) ) A_μ
            = ef_γ v_i^R/m_i^2  A_μ≡ e H_i γ A_μ ,     (i=ρ^0, , ϕ)
             
            where  m_i^2  is the renormalized squared mass of the V_i vector
            meson as given in Equations (<ref>); let us note that each
            H_i γ carries a hidden factor of 1/g. Starting with the charged pion part of
            Equation (<ref>), one gets :
            
              L_AVP^ ren(π^±)=
            eC/2 [G^R(π^±) · H_γ V]  ϵ^μναβ∂_μ A_ν[ ∂_αρ^+_βπ^- +
            ∂_αρ^-_βπ^+] ,
            
            where H_γ V is the 3-vector constructed from the H_i γ
            defined just above :
            H_γ V= ( f_γρ_0^R/m_ρ_0^2, f_γ^R/m_^2, f_γϕ^R/m_ϕ^2 )
            
            
            For the P_0 part of the Lagrangian  Equation (<ref>), the algebra is
            slightly more involved and gives :
             L_AVP^ ren(P_0)=
            
            eC  ∑_i=ρ,,ϕ
             [H_γ V G^R(P_0)]_i   ϵ^μναβ∂_μ V_i ν^R ∂_α A_β   P_0   ,   P_0=(π^0, η, ) ,
            
            with :
            
            C=-N_c g^2/4π^2 f_π c_3
            This effective Lagrangian is the tool to parametrize the e^+e^-  (π^0/η) γ
            cross-sections and the VPγ radiative decays.
            
            
            
             §.§ The e^+e^-  (π^0/η) γ Cross-Sections
             Using the AAP and AVP Lagrangians derived above in terms
            of renormalized vector and pseudoscalar fields, one can construct the amplitudes
            for γ^*  (π^0/η) γ transitions. These can be written :
            
            T(γ^*  P_0 γ)=
            -i3 α_ em/π f_π [g^2 c_3 H_γ V G^R(P_0) K_γ V
            -(1-c_3) L_P_0 ]= -i3 α_ em/π f_π F_P_0 γ(s) ,
            
            where the (co)vector H_γ V is defined in Equation (<ref>), the
            G^R(P_0) matrix in Equation (<ref>) (and  Equation (<ref>)). One has also defined :
             K_γ V= ( F_γρ^0^R(s)/D_ρ^0(s), F_γ^R(s)/D_(s), F_γϕ^R(s)/D_ϕ(s) )  ,
            
            the functions F_γ V^R(s) and D_V(s) having been  defined
            in Subsection <ref>. Finally, one has :
            [     L_π^0= g_π^0 γγ/3   , and  
            L_η= g_ηγγ/√(3)   , ]
             with g_P_0 γγ, derived in <cit.>,
             and reminded in Appendix <ref>.
            
             In the case of the e^+e^-  (π^0/η) γ annihilations, the available data
             are always cross-sections which are related to F_P_0 γ(s) just defined by :
             σ(e^+e^-  P_0 γ,s)= 3 α_ em^3/8 π f_π^2 [  s-m_P_0^3/s ]^3 |F_P_0 γ(s)|^2
            
             §.§ The e^+e^- π^0π^+π^- Cross-section
              The amplitude for the  γ^* π^0π^+π^- transition
            involves most of the FKTUY Lagrangian pieces; it can be written :
            
             T(γ^* π^0π^+π^-)=T_APPP+T_VPPP+T_VVP  ,
            
             labeling each term  by the particular piece of the FKTUY Lagrangian
            from which it originates. As already noted, because  c_3=c_4 is assumed,
            there is no T_AVP piece.
            
            
            The T_APPP contribution to the full  T(γ^* π^0π^+π^-)
            amplitude can be derived from the information given  in Appendix <ref>; it
            writes :
            [ T_APPP= C_APPP [1  -G(δ_P)  ]
            ϵ_μναβ  ϵ_μ(γ) p_0^ν
            p_-^α p_+^β , ]
            where one has defined :
            [        G(δ_P)=  Δ_A/2+
            ϵ_1 sinδ_P -ϵ_2 cosδ_P  , C_APPP= -ie/4 π^2 f_π^3 [ 1-3/4 (c_1-c_2+c_4)  ] ]
            
            Three pieces are coming from the VPPP Lagrangian also displayed in Appendix  <ref>;
            they are collected in :
            [                                                                             T_VPPP= C_VPPP [ ∑_V=ρ,,ϕF^e_Vγ(s)/D_V(s)  g_Vπ^R (s)  ]
            ϵ_μναβ  ϵ_μ(γ) p_0^ν p_-^α p_+^β , ]
            where the renormalized vector couplings   g_Vπ^R (s) to 3 pions have to be derived
            from using Equations (<ref>) and Equation (<ref>),  R(s) being given by
            Equation (<ref>). The V-γ amplitudes F^e_Vγ(s) have been constructed
            in Subsection <ref> and the inverse propagators in  Subsection <ref> for
            the ρ^0 meson. The inverse propagators for the   and ϕ mesons have
            been discussed and defined in  Section <ref>. One has also
            defined :
            
            C_VPPP= -3ige/4 π^2 f_π^3 [
            c_1-c_2-c_3  ]
            
            
             The VVP Lagrangian piece in Equation (<ref>) has
            been rewritten in terms of renormalized fields in Equation (<ref>) with
            (renormalized)  couplings derivable by means of Equation (<ref>).
            The simplest way to write down T(γ^* π^0 π^+ π^-)  in a way
            easy to code within our global fit procedure is displayed just below.
            
            One first defines the H_i(s) functions :
             {[                                       H_0(s) =1   , H_1(s) =  1/D_ρ(s_+-) +
            1/D_ρ(s_0+) + 1/D_ρ(s_0-) ,;                               H_2(s)=1/D_ρ(s_+-)  ,    H_3(s)=α (s_+-) 
             [ 1/D_ρ(s_+-) - 1/D_(s_+-) ] , ].
            
            where s is the incoming squared energy and the s_ij's indicate the invariant mass squared of the
            corresponding outgoing (i,j) pairs; the primed mixing angles have been defined in Subsection <ref>. 
            T_VVP depends on the 3 functions (H_i(s),  i=1 ⋯ 3
            with s-dependent coefficients F_i(s) given below.
            
            Collecting all terms, the full amplitude writes :
             
             T(γ^* π^0 π^+ π^-)=
             [ F_0(s) H_0(s)+ C_VVP∑_i=1⋯ 3  F_i(s) H_i(s)  ]
            ϵ_μναβ  ϵ_μ(γ) p_0^ν p_-^α p_+^β ,
            
            with :
              C_VVP= -i 3 e g m^2 /8 π^2 f_π^3 (1+Σ_V) c_3
            In this way,  to write down the full amplitude, the various F_i(s) functions only depend on
            the incoming off-shell photon energy squared s; the dependence upon the various sub-energies s_ij
            is, instead,  only carried by the H_i(s) functions as clear from Equations (<ref>). One has :
            
            {[ F_0(s)= C_APPP [1-G(δ_P) ]
            +C_VPPP [F^R_ργ(s)/D_ρ(s)  g_ρπ^0 (s)
            +F^R_γ(s)/D_(s)  g_π^0 (s)+
            F^R_ϕγ(s)/D_ϕ(s)  g_ϕπ^0 (s)
             ]
             ,;                 F_1(s) = α(s)F^R_ργ(s)/D_ρ(s)
            +  [ 1+2ξ_0+ξ_8/3 ]F^R_γ(s)/D_(s)
            + [√(2)/3 (ξ_0-ξ_8)  +γ(s) ]
            F^R_ϕγ(s)/D_ϕ(s) ,;                                        F_2(s) =  [ ϵ_2 cosδ_P-
            ϵ_1 sinδ_P - Δ_A/2 ] F^R_ργ(s)/D_ρ(s)
            + 3ξ_3/2F^R_γ(s)/D_(s)
            
             ,;                                                                                                       F_3(s) = F^R_ργ(s)/D_ρ(s) ].
            
            where α(s)=α(s)+ψ_ and γ(s)=γ(s)+ψ_0.
            
            One should note, as in BHLS <cit.> and earlier in <cit.>,  that all parameters
            already fitted with the (five) e^+e^-  (ππ/K K/π^0 γ / ηγ)
            processes fully determine all  the (F_i(s),  i=1 ⋯ 3). The only new parameter
            coming with e^+e^- π^0 π^+π^- is c_1-c_2,  only affecting F_0(s); for practical purposes,
            it is handled in a specific manner, as in the previous versions of the (broken) HLS model just quoted.
            
            Indeed, a global fit to all cross-sections
            but e^+e^- π^0 π^+π^- allows to yield the relevant  parameters
            with a good approximation;  thus,  having at hand all ingredients  defining the F_i(s)'s,
            a first minimization run including the e^+e^- π^0 π^+π^- cross-section
            can be performed to also derive a first estimate for c_1-c_2. The output of this  minuit
            minimization run is then used as input for a next minimization step. This initiates an iteration
            procedure involving all cross-sections which is carried on up to convergence –
            when some criterion is met.
            
            This method has been proven to converge in a very few minimization steps <cit.>. What makes such a
             minimization procedure unavoidable is that the e^+e^- π^0 π^+π^- cross-section expression implies
            to integrate over the s_ij's at each step. This is obviously  prohibitively computer time consuming for a negligible gain.
            Hence, at each step, one starts by tabulating  coefficient functions,
            exhibited in the next expression between curly brackets :
            
            [                                                                      σ(e^+e^- π^0 π^+π^-,s)=                                       α_ em s^2/192 π^2× [
            {∫ G dx dy } |F_0(s)|^2 . [0.3cm]                             .
            + C_VVP^2 ∑_i,j=1⋯ 3  F_i(s) F_j^*(s)
            {∫ G H_i H_j^* dx dy } .;                                                                                              . 
            +C_VVP∑_i=1⋯ 3 ( F_0(s) F_i^*(s) {∫ G H_i^* dx dy}
            + F_0^*(s) F_i(s) {∫ G H_i dx dy} )
             ] ]
            and these tables are used all along this step. Equation (<ref>) uses
            the Kuraev-Silagadze parametrization <cit.>
            and its kernel G(x,y) function; these are reminded in Appendix H of <cit.>.
            
            
            
            § THE DATA SAMPLES SUBMITTED TO GLOBAL FITS
             Basically, 48 experimental data samples are presently available which enter
            the scope of the HLS Model.  Relying on the global fits performed in
            <cit.>, one has been led to discard a few of them
            from our fitting procedure; this situation happens again within the present framework.
            Before going on, let us remind  the most important part of the data samples and list
            the newcomers. Newcomers encompass either newly collected samples or
             existing form factor spectra covering spacelike momenta not treated within the former BHLS framework
             <cit.>.
            
             Substantially, the set of available  data samples includes :
              
              *  i/ For the  π^+π^- annihilation  channel :
            
              The earliest
              data samples date back to the beginning of the eighties <cit.>. These have been
              followed by quite precise data samples collected in scan mode by
              CMD-2 <cit.>   and SND <cit.>  on
              the VEPP-2M  collider at Novosibirsk. The published spectra cover the   energy region from 370 MeV
              to 970  MeV; presently,  there is still no published spectrum covering the ϕ(1020) energy region.
            
              These measurements, referred to globally as NSK, have been followed by large statistics data collected
              using the so-called initial state radiation (ISR)  mechanism (see <cit.>, for instance).
              In this way, three data samples <cit.> have been collected with
              the KLOE detector running on   the DAPHNE collider up to  ≃ 1 GeV, just below the ϕ meson mass.
               In the same period of time, another π^+π^-  sample
               has  been collected by BaBar <cit.> covering the energy range up to 3 GeV. By 2015,
               BESIII has also published a  π^+π^- sample <cit.> and, by end
               of 2017, a group from the CLEO-c Collaboration has published the latest to date π^+π^-  sample <cit.>
               collected with the CESR detector.
            
               As a whole, the most precise data samples represent now 8 spectra. One has been led to discard
               from our global framework  KLOE08 <cit.> and BaBar  <cit.>
               for different reasons[As for the BaBar spectrum, a recent comprehensive study <cit.>
               analyzing the ω→π ^+π ^- branching ratio reaches conclusions similar to ours.
               One also has to keep in mind that the BaBar ππ spectrum is in conflict with the τ data corrected 
               for the isospin breakings. This is what we learn from Figure 1 of Ref. <cit.>.]
               discussed in <cit.>. A reanalysis  unifying the three KLOE data
               samples has been recently performed <cit.> and an 85 data-point spectrum has been
               derived which will be commented at the appropriate place. Moreover, the full error covariance
               matrix including the correlations between KLOE08, KLOE10 and KLOE12 is now available; this is used
               here  to account for the (weak) correlations between KLOE10 and KLOE12.
               
              *  ii/ For the  K^+K^- annihilation  channel :
            
             	Up to very recently, the available data samples covering this channel were the three
            	scans collected  by CMD-2   <cit.> and the two scans of SND
            	 <cit.>. They were included within the global fits performed in the (previous) BHLS framework
            	 (see <cit.>).  Quite recently, the CMD-3 Collaboration has found that the two CMD-2
            	 scans  <cit.> should undergo a rescaling by <cit.>1.094 ± 0.04 which represents
            	 an important correction of the spectrum absolute normalization and of its (correlated)
            	 systematics (2.2%→ 4.45%). The influence of this correction on the
            	 physics information previously
            	 derived within BHLS in <cit.> imposes the quite significant update
            	   performed in the present study.
            
            
            	BaBar has published in 2013 the first measurement <cit.> of this cross-section performed
            	in the ISR mode. This high statistics spectrum covers the full energy range from
            	threshold  to 3 GeV. Moreover,
            	very recently, CMD-3 has also produced a new high statistics measurement of this
            	spectrum   in scan mode <cit.>. However, as emphasized  in the CMD-3 publication,
            	this measurement is inconsistent with the BaBar data sample <cit.>. This issue
            	is addressed  below.
                
              *  iii/ For the K_LK_S annihilation  channel :
            
             	To our knowledge, the first reported measurement of the e^+e^- K_LK_S cross section has
            	been provided by CMD-2 in 1995
            	<cit.> and has been followed in 2004 by the 4 scans reported in <cit.>. In the meantime
            	SND also produced 4 data samples, 2 in the charged decay mode of the K_S meson, 2 in
            	its neutral mode <cit.>. These were the data encompassed in our former studies. The present work
            	also includes the CMD-3 data sample <cit.> recently published.
            
                 
              *  iv/ For the (π^0/η)γ annihilation  channels :
            
             		As for the e^+e^-π^0 γ cross section, in our energy range,  the
            	 available data samples are scarce : One sample has been provided by CMD-2  <cit.>
            	 and another one by SND <cit.>– which supersedes the former
            	<cit.> used in our previous studies<cit.>;
            	the older SND π^0 γ spectrum given in <cit.> is also considered.
            
            	Regarding the e^+e^-ηγ cross section, the situation is more favorable as
            	six data samples corresponding to different η decay modes have been collected
            	by CMD-2  <cit.> and SND <cit.>.
            	Nevertheless, no newly collected data sample for this channel  has been reported.
            
             	  
              *  v/ For the  π^+π^-π^0 annihilation  channel :
            
            	As no newly collected  sample covering the energy region up to the ϕ
            	mass has been reported, we are left with only
            	the samples already examined in our previous analyses <cit.>.
            	These  have been collected by CMD-2  <cit.>
            	and SND <cit.> and cover either the  peak region
            	or the ϕ peak region. We also consider the  3-pion spectrum given in <cit.> as it
            	mostly
            	deals with the region in between the  and ϕ peaks, allowing for a valuable constraint
            	on the ρ^0 meson physics background.
            
            	However,  preliminary fits showing that the sample
            	in  <cit.> exhibits an average χ^2 per data point much above 2 – as also yielded by
            	 <cit.> in a different theoretical context – we have been led to discard it from our fits.
            	 Overall, one is left with nine data samples, including an old sample from the former CMD
            	 Collaboration.
            
             	  
              *  vi/ The pion and kaon form factors in the spacelike region :
            
            	One of the motivations to develop a new breaking procedure for the HLS Model is to
            	design it in such a way that the resulting broken model (BHLS_2) could encompass
            	accurately the close spacelike region
            	of the pion and kaon form factors.
            	Indeed, a fair  account of F_π(s) in a region which stretches on both sides of
            	s=0 is the best way to
            	ascertain the low energy behavior of the e^+e^- π^+π^-  cross section
            	 which provides an enhanced contribution to (g-2)_μ.
            
            	The NA7 Collaboration has measured the pion form factor <cit.> in the spacelike range
            	Q^2 ∈ [-0.253,-0.015] GeV^2  with good statistics and a small
            	reported normalization error (0.9%); two years later, NA7 has also published
            	a measurement of the charged kaon  form factor <cit.> over the range
            	Q^2 ∈ [-0.0950,-0.0175] GeV^2 with also a small reported normalization uncertainty
            	of 1%. On the other hand, an experiment at Fermilab had previously published less
            	precise spectra on the pion and kaon form  factors <cit.>
            	over resp. Q^2 ∈ [-0.092, -0.039] GeV^2 and  Q^2 ∈ [-0.1145,-0.0409] GeV^2.
            	Unlike  former experiments[ Reference
            	 <cit.> has collected  the results of several of them.]  at higher
            	virtualities and, more modern experiments[Former references can be found
            	in the quoted papers.]<cit.>,
            	Refs.  <cit.>– and  <cit.> as well – report on
            	results derived in
            	resp. π and K  scattering on atomic electrons and are, thus, model independent.
            
            	As a matter of fact, one should also note that there is no K^0 form factor spectrum
            	measurement in the spacelike region, which would be a real challenge; nevertheless,
            	 measurements  of <r_K^0^2> exist  <cit.> which  report negative values.
            
            
            
            	A common feature of the data samples we have to deal with is that their systematic	
            	errors are quite generally reported as correlated; this implies that the corresponding
            	error covariance matrices W can be of large dimension (especially for the BaBar
            	and KLOE samples). As the input to the χ^2 expressions which feed our fitting code
            	are just their inverses W^-1, the issue has been to ascertain the quality of these. For this
            	purpose, we have chosen using the subroutines of the  NAG library. Numerical checks have 
            	nevertheless been performed on the eigenvalues, the diagonal and non–diagonal elements of the products 
            	W · W^-1 to verify that numerical departures from unit matrix properties 
            	were indeed negligible.
            	
            
            
            
            § COMMENTS ON THE MODEL PARAMETERS
             Before going on, it is worthwhile commenting
            on the fit parameters of the broken HLS Models.
            Basically, the number of parameters of the original (unbroken)  HLS Model
             <cit.>
            to be fixed from data is small. As for its non-anomalous sector, besides f_π, these are
            the universal vector coupling g and the specific  HLS parameter
            a, generally found to slightly depart from 2, the value expected
            from standard VMD approaches. Its anomalous sector <cit.> involves a
            few more parameters already
            displayed in the Lagrangian pieces presented in Section <ref> :
            c_1-c_2, c_3 and c_4. There
            are strong motivations <cit.> to impose c_3=c_4.
            So, when unbroken, the HLS model depends on only four unconstrained parameters.
            
            
            
            
            There are two parameter categories introduced in our approach: the coefficients of
            the subtraction polynomials which
            supplement the corresponding loops and, on the other hand, the breaking parameters inherently
            affecting the present Lagrangian – and those defined in <cit.>.
            
            Traditionally, as illustrated by the Gounaris-Sakurai formula <cit.>, the coefficients of the
            subtraction polynomials are reexpressed in terms of the resonance mass and width, which are determined by
            a fit to the data. In a global approach like ours, the subtraction polynomials affect loops which
            come  simultaneously  in  the various amplitudes associated with the various annihilation processes
             simultaneously under examination. In this case, as there is no reason to favor a specific channel
            to fix the meaning and the value of parameters which come in all channels, the strategy adopted is to let the data,  as a whole,
             determine the subtraction parameter values via a  global fit. Even if the parameter
            content may look less intuitive, this approach – also adopted in BHLS <cit.>–
            looks the most motivated. This does not prevent to extract usual physics quantities,
            relying on the fit results, as done since <cit.>, for instance.
            
            
            In the second category of parameters, besides the parameters of the unbroken HLS model reminded
            just above,
            one finds all those introduced by the various aspects of the breaking procedure at work.
            As for the parameters generated by the BKY mechanism and reminded in Section <ref>,
            six in total, the condition Σ_A=0 is already stated since the original BHLS Model
            <cit.>. It has been shown in Section <ref> that the
            diagonalization procedure of BHLS_2 imposes Δ_V=0. Therefore, four BKY parameters only
            (Σ_V, Δ_A, z_A and z_V) remain free.
            
            In Sections <ref> and <ref>, it was also shown that
            ξ_0=ξ_8 comes as a natural constraint on the parameters of the CD
            breaking defined in Section <ref>. Therefore the CD breaking involves
            two more parameters[Rigorously speaking, as already noted, two more breaking parameters
            can be defined addressing the charged and neutral K^* mesons.]ξ_0(=ξ_8) and ξ_3. Finally, switching on the Primordial
            Mixing  introduces three more parameters (ψ_0, ψ_ and ψ_ϕ) into the fit
            procedure.
            
            
            
            § THE BHLS_2 GLOBAL FITS : GENERAL FEATURES
             Various kinds of global fits have been performed using the data
            samples and channels listed in Section <ref>, with or without the
            spacelike data, with or without the τ data, fixing some of the model
            parameters (the ψ_α's angles and z_V) or leaving all of them free.
            
            When mixing as many data samples carrying important –  sometimes dominant –
            overall normalization uncertainties, the effect of the so-called Peelle's
            Pertinent Puzzle (PPP) <cit.> cannot be ignored to avoid getting  biased
            estimates for physics quantities. The use of χ^2 minimization methods,
             first questioned, has been finally justified  <cit.>; however,
            when dealing with distributions or form factors, the pertinent method turns out
            to invoke the  underlying (theoretical) model function which is just what
            fits are supposed to provide. Following <cit.>,  <cit.>
            has defined an iterative procedure proved to
            cancel out biases (see the Appendix in <cit.>).
            
            
            It is worthwhile briefly sketching the conclusions of[See, in particular,
            its Section 4. ]<cit.>  concerning
            the issue just reminded :  1/ The normalization uncertainty (denoted σ,
            possibly s-dependent) of any spectrum T is absorbed within the total covariance
            matrix contributing to the general χ^2,
             2/ The correction which affects the normalization  of T is
             derived from the data, the theoretical function provided by the fit, σ
            and the  statistical covariance matrix – which may absorb the uncorrelated systematics.
            
            Let us make a few general statements about fit results and data before focusing on specific
             topics, mainly the K K and τ sectors and the spacelike region behavior of the
             pion and kaon form factors. The detailed properties will be emphasized in the next Sections.
            
             
              * 
             As for the behavior of the data samples within our global fits, one can state
             that the issues met with the dipion spectra from BaBar <cit.>
             and KLOE08  <cit.> within the
             BHLS framework <cit.> are confirmed within the BS and RS variants of
             the BHLS_2 framework studied here; they are also discarded in the
             present study.
            
            Whatever the specific conditions of the various fits performed,
            the  π^+π^-  data samples from NSK ( i.e. CMD-2 and SND), KLOE10/12,
             BESIII  <cit.> and  CLEO-c <cit.> are as well described within the
             BHLS_2 framework as they were already within BHLS  <cit.>; more precisely,
             NSK and KLOE10/12 contribute an average χ^2 per data point of ≃ 1, while
             the BESIII and CLEO-c samples yield χ^2/N_ pts≃ 0.6.
             Figure <ref> displays the
            spectra and fits for  the pion form factor in the e^+e^- annihilation and
            in the τ decay.
            
            
              * 
            In order to account for the difficulty to fully address the
            π^+ π^- π^0 channel in the ϕ mass region within
            BHLS, a so-called  B model was defined which simply turns out to discard
            from fit this mass region for (only) this channel <cit.>.
            In contrast, BHLS_2  does not meet any issue[As for
             the  discarded data sample <cit.>, see the above Section <ref>
            and also  <cit.>.]  and achieves
            a fair description of the π^+ π^- π^0
             data up to, and including, the ϕ mass region.
              The spectra are shown in Figure <ref>.
            
              * 
            Finally, the e^+e^-  (π^0/η) spectra are also nicely fitted within BHLS_2– as well as
            in the former BHLS – as illustrated by Figure <ref>. Quite generally,
            one obtains χ^2/N_ pts≃ 90/112 and χ^2/N_ pts≃ 120/182 for resp.
            the π^0 and η channels.
             
            
            
            In order to perform global fits, one obviously should make an assumption on the energy
            calibration of each of the data samples considered. For the dipion spectra, one
            has gathered into the fit data coming from
            the VEPP-2M machine at Novosibirsk (CMD-2 and SND), data from DAPHNE (KLOE10 and KLOE12),
            from BESIII and from CESR.
            A significant mismatch of the relative energy calibrations of these would have produced
             a failure of their common fit because of the
             ρ- drop off region; this is clearly not seen neither in the χ^2
            contributions of the various samples, nor in  Figure <ref> (see the inset in
            its leftmost panel). Moreover, comparing both panels in this Figure excludes a significant
            energy calibration  mismatch between e^+ e^- and τ data.
            
            Moreover, the data collected by CMD-2 and SND around the  and ϕ peaks
            in the 3-pion channel as well as in the (π^0/η) final states indicate a
            good consistency
             between the energy calibrations of the various data samples collected
             in different  runs distributed over several years. This is visible
             in the fit properties as well as in the
             Figures where the fit function and the various data are displayed.
            
            
            § THE BHLS_2 GLOBAL FITS : THE E^+E^-  K K CHANNELS
             The available data covering the  K^+K^-
             channel[In the following,
            as already noted in Section <ref>, the two CMD-2 scans <cit.>
            are corrected as  stated in <cit.>.
            As for the BHLS predictions in <cit.>,
            they are updated correspondingly.]
            come from SND, CMD-2 and BaBar. Recently CMD-3 has published a new K^+K^- measurement
             <cit.>.
            As for the neutral channel, most data come from CMD-2 and SND and have been listed in
            Section <ref>. A new data sample covering our fitting range has also been
            produced by CMD-3 <cit.>; unfortunately, the published K_L K_S spectrum from BaBar
            starts at 1.06 GeV, just above our fitting range <cit.>.
            
            
            
             §.§ Preliminary Remarks
            
            As for the K^+K^- channel :  In its <cit.>, the  CMD-3 Collaboration points
            toward disagreements with both the BaBar <cit.> and the original CMD-2 data samples
            <cit.>.  Once the correction factor (1.094 ± 0.040) for the trigger
            efficiency of CMD-2 is applied, the discrepancy between  the CMD-3 and CMD-2 K^+K^-
            samples looks  much reduced <cit.>– possibly washed out.
            Anyway, a full examination of
            the various e^+e^-  K K data samples[As for the contribution
            to the muon HVP 
            of the energy region from s ∈ [4m_K^2, 1.06 GeV^2] to a_μ(K^+K^-), BaBar
            obtains 18.64 ±0.16_ stat± 0.13_ syst± 0.03_VP while CMD-3 gets
            19.33 ± 0.040 in units of 10^-10, a difference of ≃ 0.7 ×10^-10
            between  the central values, a 3.25 σ effect.] is worth   addressing.
            
            On the other hand, preliminary BHLS_2 analyses based on the CMD-3 and BaBar data,
            made us aware of a visible
            degradation of  both  K K channel data description above
            1.025 ÷ 1.030  GeV;
            therefore, in the present work, one limits the fit region for the  K K data
            to the √(s)∈ [2m_K, 1.025 GeV] range. The explanation for this effect
            is unclear. It could be
            an intrinsic weakness of the BHLS/BHLS_2 models[The standard
            HLS framework <cit.>, which underlies the BHLS/BHLS_2 models, does not account
            for vector particles  not belonging  to the fundamental vector nonet.] which, as it stands,
            does address  a possible  onset of the high mass vector meson contributions; it could
            also reflect an experimental issue. For the sake of consistency, we have chosen to impose
            the same fitting range also to the  former CMD-2
             and SND  data samples in both the charged and neutral modes,
            
            
            
             §.§ A Specific Issue of the K K Data
             As already pointed out, one does not detect any indication for a mismatch
             in the energy calibration of the various data samples covering
             the 3-pion and (π^0/η) channels, all collected with the CMD-2 and SND
             detectors on the same VEPP-2M collider. One also does not  observe any mismatch
             between the energy calibration of the 3-pion and (π^0/η) data and those
             collected by the same collaborations in the K K channels.
             Nevertheless[E. Solodov, Budker Institute, Novosibirsk, private communication.],
             CMD-2 and SND  carry a (correlated) absolute energy uncertainty in the range
            of 30÷ 40 keV, while for
            CMD-3, the  absolute energy uncertainty is ≃ 60 keV, statistically
            independent of that of CMD-2 & SND. Moreover, BaBar also reports
            a 50 keV correlated energy scale uncertainty<cit.>.
            When dealing with an object as narrow as the ϕ meson within a global framework,
            these uncertainties should be
            accounted for otherwise, the fits cannot succeed.
            
            
            As the bulk of data samples for both KK channels has been provided by CMD-2
            and SND, it looks quite natural to choose their common energy
            calibration as a reference; this will
            be denoted  E_ NSK in the following. For this purpose, one
            defines  E_ CMD3=E_ NSK+δ E_ CMD3 and E_
            BaBar=E_ NSK+δ E_ BaBar to
            recalibrate the CMD-3 and BaBar data point energies in a fully correlated way and derive
            their E_ NSK equivalent energy.
            
            Therefore, in order to deal with the  K K data, we are faced with
            a case where two kinds of correlated uncertainties have to be accounted for: Energy
            calibration and absolute normalization of the cross-sections. Hence,  one merges
            in the minimization procedure the method developed in <cit.> and the
            two calibration parameters δ E_ CMD3 and δ E_ BaBar to be also determined
             within the same fits.
            
             For this purpose, one assumes that there is no  functional  relation
             between   the normalization uncertainty λ_i– of an experiment labeled i–
             and the energy scale uncertainty  δ E_i; then,
              as ∂λ_i/∂ [δ E_i] =0 for each experiment,
             substituting, as emphasized in <cit.>, the value  for λ_i derived from
             ∂χ^2 /∂λ_i =0 into the expression for χ^2 is unchanged
             compared to  <cit.>.
            
            
            
             §.§ Fitting the K K Channels : NSK, BaBar and CMD-3 Separately
            
            In the following, the non-kaon sector data and channels used in our fits
            include, unless otherwise stated, the π^+π^-, π^+π^-π^0, π^0, η annihilations,
            the τ decay to π^±π^0 and a few decay partial widths (see <cit.>).
            One starts by discussing the  fits where the CMD-2 & SND (NSK),
            BaBar and CMD-3 kaon data are not influencing each other within the minimization procedure.
            Table <ref> gathers most of our results in this case; these are
            commented right now.
            
            
              * 
            With regard to  NSK :  <χ^2>, the average (partial) χ^2 per NSK data points, is of order 1
            in both the neutral and charged kaon channels; thanks to the fair description of the
            other channels involved in the fit, the global fit probability exceeds the 90% level.
            
            As reminded above, the procedure defined in <cit.> allows us to derive from
            the fit procedure  the normalization (absolute scale) correction to be applied to each
            spectrum[One should stress that such a "renormalization" is the
            natural outcome to the scale uncertainty problem and that the main advantage of
            any kind of global fit is to allow determining such a correction in consistency
            with a (large) set of data samples coming from different sources.].
            For the charged channel, one gets +5.0% for SND (expect. value 7.1%), a 0.7
            σ effect, and  –1.6% for CMD-2 (updated expectation 4.4%),
            a 0.3σ effect; in the neutral channels, one correspondingly
            obtains –1.8%  (expect. value 4.2%) for SND and +0.6% (expect. 1.7%) for CMD-2, resp.
            0.4 σ and 0.3 σ effects.
            
            Therefore the CMD-2 and SND data samples successfully fit  the BHLS_2 framework and,
            moreover, the derived scale corrections meet the reported scale uncertainty expectations.
            
             
              * 
            With regard to  BaBar : One has chosen to supply the NSK neutral data to the fits
            in order to populate each dikaon channel as the BaBar K_L K_S data <cit.>
            do not cover our fitting energy range.
            As the BaBar charged spectrum <cit.> has been
            collected in the ISR mode, it looks appropriate to compare, within the fit,
            each spectrum datum in a bin
            to an average value of the model spectrum over the bin width.
            The second data column in Table <ref> shows
            that the BaBar sample yields <χ^2>=1.7, somewhat large, and a global fit probability also
            above  the 90% level.
            
            
            The minimization procedure returns as energy shift vs NSK δ
            E_ BaBar=-125.8± 19.0
            keV; this looks fair when compared to the independent 30÷ 40 keV (NSK) and 50 keV (BaBar)
            energy calibration uncertainties. Moreover, the spectrum scale correction
            is at the 10^-4 level – so that the absolute scale of the BaBar spectrum has not to
            be corrected. On the other hand, the scale corrections for the
            neutral NSK data included in the fit  are –2.1% (SND) and +0.7% (CMD-2), almost identical to
            what has been obtained in the fit of the NSK data in isolation mode.
            
             Interestingly
            one can also derive the scale correction for the charged data of
            SND (+5.2%) and CMD-2 (–3.7 %) even if they are not submitted to the fit; these predictions
            are clearly in the expected  ballpark. This indicates, beforehand, that NSK and BaBar data
            in the kaon sectors are  consistent with each other.
            
            
              * 
            As for CMD-3 : This Collaboration has provided data samples in both dikaon final states,
            When fitting these data, the correspondence is made – as with NSK –
            between each datum and the theoretical model computed at the nominal energy. As for NSK and
            BaBar,
            the systematic uncertainties are treated as fully correlated and serve to build up the
            systematic  covariance matrices to be added up  to the statistical ones; these are
            used in the fitting of the neutral and charged CMD-3 data within the global BHLS_2 framework.
            The fit results are displayed in Table <ref>.
            
            The  <χ^2> for both CMD-3 data samples are clearly very large and the global fit
            probability poor. The scale corrections found are +1.2% (charged channel)
             and -1.2%  (neutral channel) to be compared to the reported r.m.s., resp. 2.0%
             and 1.8%. The possibility of treating the bin information as in ISR
             experiments[Comparing the datum information to the fit function averaged over
             the bin width, as just sketched  for BaBar.] has been examined – specifically, only in the neutral
             channel – and in  view of the outcome, we gave up.
            
             For the sake of completeness, we have nevertheless also performed the fit of the CMD-3 data in isolation
             by discarding the  non-diagonal part  of the total covariance
             matrix[This turns out to consider that the systematic uncertainties are
             fully uncorrelated; some configurations intermediate between full correlation
             and no correlation can lead to similar results.], as seemingly done in
             <cit.> for the BaBar data; the fit results change dramatically as one reaches
             <χ^2>=14/16 and <χ^2>=17/17 for resp. the charged and neutral modes
             of the CMD-3 data and  with a fit probability exceeding 90%.
             Moreover, in this case (CMD-3 data in isolation) the corrections
             to the original  normalization of the cross-sections become negligible (≃ 10^-4)
             for both the charged and neutral channels[The <χ^2>'s   of the  fit prediction
             to NSK data remain in reasonable correspondence with the NSK fit in isolation (55/49 and 105/92)
             while the distance of this prediction to BaBar data degrades more severely (76/27).].
             As for the energy shift vs NSK it moves from -51.4 ± 13.7 keV (see Table <ref>)
             to  -44.6 ± 14.8 keV, still in the expected ballpark.
               
            
             §.§ Fitting the K K Channels : NSK, BaBar and CMD-3 Combined
             This Subsection is devoted to analyzing the consistency of the
            three available groups of data. Quite generally, the covariance matrices used within the fit
            procedure assume full correlation for the systematics, as in the previous paragraph.
            Nevertheless, one  also briefly reports on additional information about a fit performed assuming
            the systematics uncorrelated for CMD-3 data.
            It is worthwhile proceeding with pairwise combinations.
            
              *  As clear from the first data column in Table <ref>– and as could be
            inferred from the results already reported –  one observes a fair consistency between
            NSK and BaBar samples within the global BHLS_2 framework. Comparing the present NSK and BaBar
            contributions to the (total) χ^2 with their analogs  in their  fits in isolation clearly shows
            that there is no significant tension between them. The only noticeable difference is
            the central value for δ E_ BaBar found closer to 100 keV by 1σ.
            
            Other pieces of information concerning the absolute spectrum normalization[
            For SND and CMD-2, as their corrections are computed run by run and sample by sample,
            here we only give average values, for simplicity.]
            are worth to be mentioned. The fit outcome indicates that
             the BaBar data normalization should be downscaled by ≃ 0.4 %,
             the  (updated) CMD-2 charged spectrum by 1% while the SND charged data normalization
             should be increased by 8.1%. These corrections are always  consistent
            with the respective expectations reminded above.
             
              *  The second data column in Table  <ref> displays results from fits
            combining NSK and CMD-3. Here also one observes a striking resemblance between the individual
            partial χ^2's in the combined fit and in their fits in isolation,  i.e. 
            the NSK
            data samples go on yielding  <χ^2>_ NSK≃ 1
            while  one gets <χ^2>_ CMD3 > 6; the global fit
            probability is also poor.
            
            Besides the energy shift δ E_ CMD3=-46.5± 6.7 keV, the spectrum rescalings derived from
            the fit are interesting pieces of information. In this way,  the fit returns an increase of the
            charged CMD-3 data normalization by only 0.2% while the neutral CMD-3 ones should be downscaled by 2.6% (a 1.5 σ effect).
            One also finds that the SND charged data normalization should undergo an increase by 9.2%,
            while for CMD-2  the increase is only 1.5%.
            As for the neutral kaon normalizations, SND should be downscaled by 2.1%,
            while CMD-2 is increased by 0.5%. So, one does not observe major changes in the NSK normalizations.
            
            
              *  The third data column in Table  <ref> displays results from fits
            combining BaBar and CMD-3 kaon data, excluding their NSK analogs. The global fit probability
            clearly indicates that such a combination is not really favored; one should remark, nevertheless,
            that all properties exhibited by CMD-2 and SND are in fair correspondence with their
            analogs yielded in their fits in isolation; additionally, even the χ^2 distances
            of the CMD-3+BaBar fit function to the NSK data (not under fit) – given within parentheses –
            are analog to their values derived in  NSK alone or NSK+BaBar fits.
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            Quite interestingly, the tension observed between CMD-3 on the one hand, and BaBar and NSK on the other hand,
            is not manifest in plots, as soon as the energy calibration of CMD-3 vs NSK is applied. For instance
            Figure (<ref>), shows the KK cross sections
            derived from the fit
            to NSK + BaBar; the corresponding data are rescaled and appropriately shifted as emphasized
            in the previous item; one has  superimposed the (unfitted) data from CMD-3 shifted by -45.9 keV
            and rescaled as  predicted by the  NSK + BaBar fit : -4.3% and -2.4% for resp. the
            charged and neutral CMD-3 data samples.
            
            On the other hand, one has found worth to also perform the fit sketched
            in the rightmost data column in Table <ref>,  i.e. treating the NSK and BaBar
            uncertainties as fully correlated and both CMD-3 data samples as carrying uncorrelated
            uncertainties. Using  the corrections to the absolute normalization
            of each sample as derived from this fit[
            This amounts to the following rescalings w.r.t. the original normalizations :
            –1.3% (CMD-3 K^+K^-), –0.5% (CMD-3 K_L K_S), +8.1%  and –0.3% (K^+K^- from
            resp. SND & CMD-2), +0.4%  and –0.5% (K_L K_S from resp. SND & CMD-2) and
            ≃ -0.3 % for BaBar K^+K^-.], one gets a figure which
             superimposes exactly to  Figure
            (<ref>), in the smallest detail.
            
            This clearly indicates that  the observed tension  is essentially due to
            the non-diagonal parts of
            the CMD-3 error covariance matrices constructed from experimental information,
            as also done for the other data samples, CMD-2 in particular.
            
            
            § FIT PROPERTIES OF THE VARIOUS BROKEN HLS MODELS
            
            
             All fits referred to in this  Section involve  the collection of data samples
             and channels  already described  and listed in Section <ref>; as for
             the K K channels,
            one limits oneself to using only the  samples from CMD-2 (updated), SND and BaBar
             which do not show any obvious tension. The motivation for this choice is to avoid
             disentangling  model dependence effects, on the one hand, and identified tensions
             among data samples[The issue between BaBar and CMD-3 charged kaon
             data has already been noted in <cit.>.] on the other hand. There is indeed
             no gain in increasing the statistics if this turns out to increase the systematics
             in an uncontrollable way.
            
            
             The BHLS_2 fits presented up to now have been derived using the so-called
             Reference Solution (RS); this solution
             combines the so-called Basic Solution (BS) defined in Section <ref> and
             the Primordial Mixing defined in Section <ref>.
            
             Nevertheless, in order to explore model dependence effects, analyzing also the
             fit results derived using BS is certainly relevant. On the other hand, for the same concern,
             we have rerun the BHLS model <cit.> using the SND together with the  updated
             CMD-2 data; it then looked appropriate to also include the BaBar dikaon sample.
            
             Therefore, with the BS and RS variants of BHLS_2 and with BHLS, one
              actually has at disposal  three approaches
             to the same data; as for BHLS, one concentrates on the so-called model B
             which circumvents the issue met at the ϕ mass in the three pion annihilation channel
             which is one of
             the motivations having led to BHLS_2. The other one was to ascertain the description
             of the dipion threshold region and to improve the accountability of the spacelike region.
             The main fit properties are sketched in Table <ref> and deserve some comments.
            
            
            
             The leftmost data column – dealing with BHLS – shows the  various χ^2 for the
             displayed  channels and the sum of the contributions of all data submitted
             to the fit. In all channels the  average χ^2 is of the order 1 or
             better[One should note the smaller number of data points considered in the fit
             procedure for the π^+π^-π^0 annihilation channel, because of excluding the
             data points from the ϕ region.]  and the (global) fit probability is fairly good.
            
             The rightmost data column displays the corresponding information derived with the RS
             variant of BHLS_2. Here also, the fit quality is good and one yields a fair probability.
             One should note the quite favorable <χ^2>'s for the pion and kaon form factors in
             the spacelike energy region in conjunction with good χ^2's in their timelike regions;
             the spacelike data samples – coming from NA7 and Fermilab – are discussed
             in detail in the next Section. As far as we know, the fits presented here are the
             first ones involving the pion and kaon form factors simultaneously in their
             timelike and spacelike regions with such quality.
            
             
            
            
             The third data column displays the fit properties derived from the BS variant of BHLS_2
             using the same data as those just discussed for the RS variant. The fit probability is  
              lower than usual in the HLS context. Scrutinizing the various items displayed in Table <ref>, 
              one clearly observes a significant 
             tension in  the BS fit  absent from the RS fit results  when  the
             τ data are involved.  Indeed, when going from RS to BS  χ^2_τ increases by 6 units, but
             χ^2_3π increases by 12 and  χ^2_ spaceπ by 8; moreover[The different behavior
             of the NSK and KLOE dipion data under BS global fits  reflects the different  magnitude 
             of the uncertainties
             which makes the NSK dipion data – and likewise for those from  BESIII and CESR – more permissive than their 
             KLOE analogs.],  χ^2_KLOE 2π increases by 17 units whereas
             χ^2_NSK 2π is marginally affected (a 2 unit increase). 
              
            
            
            
            
            
            Instead,  the BS fit performed  when  discarding the τ data leads to a picture
            in good correspondence with the RS fit results already commented; indeed, 
             the partial χ^2's reported in the second data column of Table  <ref>  are similar 
             to  those obtained by fitting with the RS variant of the BHLS_2 model.
             Nevertheless,  some tension survives in the 3-pion sector while the dipion spectra 
             and the spacelike data almost recover their optimum fit quality; so, the BS
            variant of BHLS_2   exhibits a quite acceptable behavior once the τ sector is discarded.
            
             In conclusion, despite their different structures – and taking into account the peculiarities 
             of the BS variant of BHLS_2–  these models/variants open the possibility to
             examine  model dependence effects in the evaluation of  some reconstructed physical quantities; this may 
             allow to compare these with the systematics generated by dealing with contradictory samples.
            
            
            
            § THE PION AND KAON FORM FACTORS IN THE SPACELIKE REGION
            
            By allowing  different HK masses for the  ρ^0 and  mesons,
            BHLS_2 naturally permits a smooth connection
            between the spacelike and timelike sectors of the pion  and kaon
            form factors.  Moreover, applying in this context the vector meson dynamical mixing
            (see Section <ref>), does not break this smooth connection because
            α(s=0)=0 becomes automatic and supplements  the (already) vanishing properties
            of the other angles : β(s=0) = γ(s=0) = 0.
            
            
            
             §.§ The Pion and Kaon Form Factors in the Close Spacelike Region
            
            Figure <ref> displays the pion and kaon form factors in the spacelike region
            coming out from the fit. The spacelike data supplied to the global fit procedure
            are only the (model-independent) data samples  collected by the NA7 <cit.>
            and  Fermilab experiments
             <cit.>. Altogether, these experiments cover photon
              virtualities  down to ≃ -0.25 GeV^2  and  ≃ -0.12 GeV^2
              for resp. the pion and kaon form factors.
            
              The pion and kaon form factors derived from the global fit including
             the spacelike data are the red curves in Figure <ref>.
             Global fits have also been performed using the BS and RS variants
             of BHLS_2,   excluding the spacelike data. The relevant fit results
             obtained using, for instance, the BS  variant
             are displayed in Table <ref>. They correspond to the green curve in
             Figure <ref>; as for the kaon form factor
              the predicted (green) curve
             superimposes exactly to the fitted (red) one.
            
            Table <ref> and Figure <ref> obviously
            prove that  the spacelike data are not a real constraint
            in the BHLS_2 global fit approach; stated otherwise, the BHLS_2
            model supplied
            with  only  timelike data fairly well predicts the spacelike behaviors
            far inside the spacelike region for both the pion and kaon form factors.
            It is worthwhile to  note  that these analytic continuations are performed
            across energy squared distances which can be as large as  1 GeV^2.
            For instance, the kaon form factor displayed here is  the analytic
            continuation of the charged kaon form factor fitted in the ϕ mass region,  i.e.
            the curve shown on the rightmost panel of Figure <ref>. Such a property
            is rare enough to be underlined.
            
            
            
            
             Figure <ref> gives   the opportunity to comment on
             the representation of spectra affected by a normalization (scale) uncertainty;
            the NA7  spacelike data are a quite simple illustration of handling the rescaling effects.
             When a  scale uncertainty affects a spectrum, there are obviously two possible
             ways to display the result : Either rescaling the data, keeping the fit/prediction untouched
            or performing the converse.
            
            So, in Figure <ref> one has chosen to
            display the NA7 data rescaled, keeping unmodified the fit (and prediction) function(s),
            the plotted uncertainties being the statistical
            ones. As for the Fermilab data samples, taking into account their low accuracy,
            no rescaling  has been performed
            and the plotted uncertainties are just the reported total errors.
            
            One can, for once, go into  a few plotting details in simple wording.
            Following the method given in  <cit.>, the scale correction is given by :
             λ= f_i V^-1_ij [m_j-f_j]/ f_k V^-1_kl f_l +1/σ^2 ,
            
            where (NA7) σ=0.9% or 1.0%, f_i=  |F_P(s_i)|^2  (P=π^±, K^±),
            V is the statistical error covariance matrix
            and m_i is the measured datum at s=s_i.  The plotted NA7 data  m^' are
            related with the original ones m by :
            
            
            
            
                m_i  m_i^'=m_i -λ f_i  (≃ [1-λ] m_i  ).
            
            Looking at  Figure <ref>,  the precise data (NA7 pion form factor) are
            indeed observed
            to tightly follow the fit function,  once the normalization correction is applied.
            
            It is relevant to provide the precise values for the λ's derived from the
            fitted parameters. One thus finds :
            
                NA pion form factor :λ_π= (1.18 ± 0.02_syst± 0.03_fit) %
            
            
                NA kaon form factor :λ_K= (1.05 ± 0.03_syst± 0.02_fit) %
            
            which exhibit noticeable accuracies. Subsection <ref> below – and  Figure
            <ref>– illustrates  that the predicted scale correction for the pion form factor
            is supported by the LQCD data which should be intrinsically free of scale uncertainties.
            
            
            
             §.§ The Pion Form Factor : The Model-Dependent Data
             As  one cannot measure
            the pion form factor in a model-independent way for photon virtualities
            Q^2 =-s > 0.3 GeV^2, this physics is performed using the scattering
            process p e  e π^+ n, referred to in the literature as
            ^1H(e,e^'π^+)n. In this approach, two questions should be addressed :
            
             i/ Separating out the longitudinal cross section σ_L,
            
             ii/ Extracting F_π from σ_L.
            
            
            One can refer to <cit.> for  comprehensive
            accounts of the issues. As for item  i/, the method has evolved and
            improved in the course of time  from the early experiments (also reported in
            <cit.>) to the more recent ones <cit.>. Item
             ii/ addresses more deeply the nuclear physics content of the process
            in order to account at best for the π^+ p n vertex and for the off-shell
            character of the intermediate π^+. Presently, the favored approach
            is the Regge Model  developed by Vanderhaeghen, Guidal and Laget <cit.>.
            The data extracted for F_π(s<0) have thus improved in time but remain
            unavoidably model-dependent. These data samples have not been included in our fitting approach,
            but simply  compared with the extrapolation of our fit function  at s<0,
             a pure prediction below  s≃ -0.25 GeV^2.
            
            In Figure  <ref>, besides the fitted data  from NA7 and Fermilab,
            one has displayed the ( unfitted) data from the DESY
            experiments <cit.> and the lowest photon virtuality data
            collected at Jefferson Lab by the F_π Collaboration. The data
            points flagged by JLAB 2008 and DESY have been extracted from the
            latest – to our knowledge – publication of the F_π Collaboration
            <cit.>. One can observe that the updated <cit.>  DESY data points fall
            exactly on the  F_π(s) function derived from a fit using also
            the NA7 <cit.> and Fermilab <cit.> data.
             The three JLab plotted data points are found systematically
            below the fit function;  however,  most of these JLab data points can also be found in
            a previous JLab publication <cit.>. The difference between the
            values in <cit.> and <cit.>
            seems mostly coming from using different extraction methods of σ_L(Q^2)<cit.>. The data points from <cit.> seem in closer agreement with
            the continuation of our  F_π(s) than the estimates given in <cit.>. They
            are also in better agreement with what can be expected from the DESY data points
            as updated in  <cit.>; this updating seems  to deal with the extraction of
            F_π from the longitudinal cross section σ_L.
            The  behavior of (the extrapolated) F_π(s)
            beyond -1  GeV^2 is shown in Figure <ref> together with the other
            reported data.
            
            
            
            
            The content of  Figure <ref>
             gives a hint that the extrapolation of the BHLS_2F_π(s) could well be valid down to
            s ≃ -1 GeV^2. This  also reflects a ρ dominance well inside
            the non-perturbative spacelike region. This dominance is illustrated by
            fit properties : When fitting the pion spacelike data, beside the ρ contribution,
            the fit function contains also  those from the  and ϕ mesons; the fit
            returns <χ^2>/N = 55/59; however, computing the χ^2 distance
            of the fit function  amputated from the   and ϕ contributions to
            the  pion data from NA7 and Fermilab, one obtains  <χ^2>/N = 56/59.
            This indicates that the (tiny) contributions of the  and ϕ mesons
            introduce a (tiny) noise and that nothing else than the ρ contribution
            is actually required.
            
            This ρ dominance phenomenon does not happen for the kaon form factor as,
            truncating the prediction for F_K(s) from its  and ϕ contributions
            returns a  disastrous χ^2 distance; for instance, the full
            χ^2 distance of the NA7 kaon data to F_K(s) which is ≃ 14, becomes ≃ 2300
            when F_K(s) is amputated from its  and ϕ contributions. Such
            property is quite unexpected as the parametrization for the   and ϕ<cit.>, reminded in Section <ref>, is nothing but a Breit-Wigner lineshape
            modified in such a way that the inverse propagators fulfill D_V(0)=-m_V^2, these HK
            masses being those in  Equations (<ref>).
             That the analytic continuation,  over more than ≃ 1 GeV^2,
             of such a simple parametrization, with parameter values fixed
             at the ϕ mass, can provide a good   description of the existing kaon spacelike data
            may  look quite amazing.
            
             
            
             §.§ The Pion Form Factor : The Lattice QCD Data
             Some Collaborations have published information on
            the pion form factor derived from Lattice QCD (LQCD) data  at pion masses close
            to the physical pion mass value <cit.>. We compare these LQCD
            form factors to our fit function F_π(s<0) which, already, fairly well
            describes the  accurate NA7 experimental data as obvious
             from Figure <ref> and Tables <ref> or <ref>.
            
             The HPQCD Collaboration has thus produced  a parametrization[We thank L.
             Lellouch for having drawn our attention to this paper.]
            of the pion form factor <cit.>F_π(Q^2) for virtualities in the range Q^2∈ [-0.1,0.0] GeV^2. In order to perform
            a numerical comparison with our fit F_π(s<0), we have plotted their pole parametrization,
            setting the lattice spacing a at zero and using their pion radius squared
             <r_π^2>=0.403(18)(6) fm^2.
            In Figure  <ref>, our fit function for F_π(s<0) is shown by the red curve; the blue curve
            is the HPQCD parametrization at the central value for <r_π^2>. The red and
            blue curves are clearly very close to each other, if not overlapping.
            The green curves display the HPQCD parametrization
            with <r_π^2> at ± 1 σ from its central value. This indicates that the
            LQCD parametrization <cit.> is in fair agreement with the  BHLS_2
            fit function F_π(s<0) much beyond the expected HPQCD range of validity
            and validated by the DESY data points <cit.>;
            Figure <ref> clearly shows that
            the agreement between HPQCD <cit.> and BHLS_2 extends to the
             prediction, shown in this Figure by the green curve, derived
            when fitting with discarding the spacelike data.
            
            
            
            
            Six spectra for F_π(Q^2) have been published by the ETM Collaboration <cit.>;
            two of them correspond to a pion mass close to its physical value, namely the
            gauge ensembles cA2.09.48 and cA2.09.64.
             These spectra are given by the form factor values F_π(Q^2) and uncertainties
             as functions of Q^2/M^2_π. In order to compare with phenomenology and experimental
             data, one has to restore the values for Q^2   from the provided
             values for Q^2/M^2_π. 
             Ref. <cit.>, which sketches the generation of these ensembles,  indicates that,
             numerically, M_π≃ 134.98 MeV, almost identical to the RPP π^0 mass. 
             
            
            
             Figure  <ref> shows the |F_π(Q^2)|^2 spectra derived from
             the  cA2.09.48 and cA2.09.64 ensembles in  <cit.>. They are displayed
             with their reported statistical errors[It is worth noting that
             the original (unscaled) NA7 stay  below the ETM data in the spectrum range
             where they are the most precise; this should be compared with Figure 6 in  <cit.>,
             especially its right panel, which exhibits the same trend.]. 
             The χ^2 distances[Reference <cit.> does not give special 
            information about systematics and performs its own fits using uncorrelated χ^2's.
             So, following this line, we define the χ^2 distance of the ETM data to our prediction by the sum 
             of squared terms d_i^2   where : d_i= (F_ETM(Q^2_i) - F_BHLS_2(Q^2_i))/σ_ETM(Q^2_i),
            using obvious notations. The quantity ∑ d_i^2 is  not submitted to the minimization procedure.] 
            of cA2.09.48 and cA2.09.64 to
             the BHLS_2  form factor squared |F_π(s)|^2 have been computed, assuming
             the systematics
             negligible compared to the statistical errors   <cit.>; they are displayed
             inside  Figure  <ref>. One can remark that the
             average χ^2 distance per datum is 0.67 for cA2.09.48, quite a good value
             for  data
             outside the fit procedure; it is only 0.06 for the cA2.09.64 form factor.
             This can be considered as an independent confirmation concerning the smallness of the systematics
             versus the statistical errors for the considered ETM spectra.
            
             Stated otherwise,
             the data points of the cA2.09.64 form factor are almost exactly on the fit function.
             Out of curiosity, we have redone the fit in the RS variant of  BHLS_2  discarding the spacelike
             experimental data in order to get, online, the predicted form factor for any s<0 using only
             timelike data; the average χ^2 distances become  0.59 and 0.10 for resp.
             cA2.09.48 and cA2.09.64. So,
             the χ^2 distance of the ETM samples to the BHLS_2 prediction indicates that the form
             factor lineshape is quite consistent with LQCD expectations; the inset in Figure <ref>,
             by magnifying the energy region close to the chiral point, illustrates the tight agreement of the prediction
             and of the ETM data.   
             In other words, the ETM spectra provide strong support to the  BHLS_2 fits; 
              this BHLS_2 model spectrum is also
             intimately related to the dealing with the reported correlated systematics as promoted
             in <cit.>.
            
              For illustration, one has also displayed in Figure  <ref>
            the closest among the two data points from DESY and
             a part of the NA7 spectrum rescaled  and  not rescaled. As expected
            from their  χ^2 distances to the fit function,
            the ETM spectra are observed closer to the rescaled NA7 spectrum than
            to the original one.
            
            
            
             §.§ The Neutral Kaon Form Factor at s<0
             As seen in the previous Subsections,
            BHLS_2 predicts  the charged pion and kaon form factors
            in the close spacelike region with  fair accuracy.  Likewise, BHLS_2  also
            predicts the neutral kaon form factor F_K0(s); it is
            displayed in Figure <ref>.  This is
            the analytic continuation of the form factor given in Subsection
            <ref> and numerically determined within the fit procedure by the
            e^+e^-  K_L K_S annihilation data.
            One clearly observes a behavior
            quite different of those for F_K^±(s), with a negative slope down to
            almost s=0, when coming from the spacelike region.
            
            
            
            
            
            
            
            Actually, as already stated in Section <ref> one has
            F_K0(s=0) =0 up to contributions of second order in breaking
            parameters. What is displayed in Figure <ref> is
            F_K0(s) computed at the fitted central value of the fit
            parameters[This is equivalent to performing only one sampling
            on the multidimensional error covariance matrix of the fitted parameters.
            ].
            To take  into account their uncertainties, a Gaussian Monte-Carlo   has
            been run which returns the averaged central value F_K0(0) ≃ 10^-6
            or less, depending on the fit conditions. This indicates that the
            negative excursion close to zero observed in Figure <ref>
              could well be a numerical  artifact.
            
            
            
            
            As for the charge radius squared of the neutral kaon, both  BS and RS variants
            of BHLS_2 return an average value :
            
            <r^2>_K^0 = 0.00 ± 0.10    fm^2
            
            consistent with all reported measurements, in particular, those of
            the NA48 Collaboration <cit.>
            (<r^2>_K^0 = -0.090 ± 0.021 fm^2),  or those of KTEV <cit.>
            (<r^2>_K^0 = -0.077 ± 0.007_ stat± 0.011_ syst fm^2). It also agrees
            with the ChPT result obtained by Bijnens and Talavera <cit.>
            (<r^2>_K^0 = -0.042 ± 0.012 fm^2) and with the former
            measurements collected in the Review of Particle Properties <cit.>.
            We don't know about predictions for F_K0(s<0) at larger
            virtualities to which one could  compare.
            
            
             §.§ Pion and Kaon Charge Radii
             The fit parameter central values and the parameter error covariance matrix can serve
            to generate several samplings (α=1, ⋯ n) and, for each, the corresponding
            pion and kaon form  factors F_P^α(s).
            These can be used to  evaluate derived physics quantities by computing
            the  average and the r.m.s. values of the various estimates obtained from the set
            { F_P^α(s), α=1, ⋯ n}.
            
            
            Among the relevant physical properties of the pion and kaon form factors,
            their behavior close to s=0 deserves special interest. As clear from Section
             <ref>, the charged pion and kaon vector form factors are 	analytic
             functions of s and, so, can be expanded in Taylor series  around the origin.
            Neglecting terms of degree higher than 2 in s, one can write :
            
            F_P(s) =a_P +b_P s+ c_P s^2      ,  P=π^±,K^± ,
            
            where the a_P's are expected equal to 1 (up to second order terms in breakings) and the
            b_P's are related with the so-called charged radii by  b_P=<r_P^2>/6 for which several
            values have been reported.
            
            In order to have a reasonable lever arm, one has chosen a  s interval
            extending  on both sides of s=0 and
            bounded[ Considering the real part of F_P(s) permits to avoid (minor) issues
            with the π^0 γ threshold effects.]
            by ± 0.05 GeV^2. Using the values of each F_P^α(s) at three s values
            (0, ± 0.05) GeV^2,  one constructs the sets
            { [a_P^α, b_P^α,c_P^α]  , α=1, ⋯ n},   and derives
            their averages and r.m.s.'s.  Working within the RS variant of BHLS_2 and using
            the largest set of experimental data samples, one derives :
             [ b_π = [1.840 ± 0.003]  GeV^-2  , c_π = [4.155 ± 0.054]  GeV^-4  , ]
            and :
             [ b_K = [1.159 ± 0.005]  GeV^-2  ,   c_K = [1.15 ± 0.24]  GeV^-4  . ]
            
            As it  is more customary to express the slope term in units of fm^2, one displays
            our values together with external information  in Table <ref>.  Our results are given
            with two uncertainties; one is derived from the sampling on the error covariance matrix
            and is indexed by "fit". One has also  varied  the content of the data sample set
             submitted to the global fit and  compared the results to those derived with the maximal
             sample set. It is observed that the central values for <r_π^±^2>  and  <r_K^±^2>
             get shifted while the uncertainties are almost unchanged. One has also switched from the RS to the
             BS variant of BHLS_2; in this case, the central values get marginally shifted. Merged
             together, these shifts contribute an uncertainty denoted  "mod" for "model" which, however,
             reflects as much  the (weak)  tensions between the various data samples of the largest sample
             set – defined in  Section <ref>– than the model properties by themselves.
            
            As for <r_π^±^2>, Table <ref> shows an agreement of our result
            with the experimental data from NA7 and Fermilab at the 1 σ level or better.
            The result of G. Colangelo  et al.<cit.> (CHS) is derived
            from a global fit  of the timelike and spacelike  pion form factor data  based on a
            parametrization through the Omnès representation of F_π(s) taking into account
            isospin breaking effects () and prescriptions for the
            asymptotics of the ππP-wave phase shift; this estimate  is in perfect agreement
            with our own fit result. Using an elaborate method mixing the low energy ππ phase information
            and some
            form factor modulus information, B. Ananthanarayan  et al.<cit.>
            recently derived  the value flagged by ACD in Table <ref> which is very close to ours.
            The two flavor result of J. Bijnens  et al.<cit.>  flagged as "ChPT  2 flavors"
            as well as the three flavor evaluation in <cit.> and the
            most recent Lattice QCD derivations of the HPQCD   <cit.> and ETM <cit.>
            Collaborations are also in good correspondence with ours. Former LQCD evaluations for
            <r_π^±^2> can also be found in Table 22 of <cit.>.
            
            As for the curvature term, varying the fit conditions as just indicated for  <r_π^±^2>,
            one can assess the effects observed when switching from the RS to the BS
            variant of BHLS_2 and when modifying the set of data samples submitted to fit.
            Altogether, this leads to our final result :
              c_π = [4.155 ± 0.040_ mod± 0.054_ fit]   GeV^-4
            
            The 2 flavor ChPT estimate based on all data samples available
            at that time (1998), <cit.> reports c_π= 3.85 ± 0.60 GeV^2
            and the 3 flavor evaluation <cit.> (c_π= 4.49 ± 0.28 GeV^2)
            look in fair agreement with the BHLS_2 evaluation given in Equation (<ref>).
            
            
            
            As the  pion and kaon form factors, in the close spacelike
            region, are frequently parametrized by a pole expression :
             
            F_P(s) =1/1 -b_P s   ,  P=π^±,K^± ,
            
            we performed alike without, however, introducing
            additional normalization factors as commonly done as indeed  we believe  in
            the relevance of the absolute scale of the F_P(s)'s returned by the BHLS_2 fits.
            Such a parametrization turns out to impose the constraint c_P=b_P^2 to the
            curvature term. This leads to :
            <r_π^±^2> =0. 424 ± 0.002_ mod± 0.001_ fit  and  
            <r_K^±^2>=0.268 ± 0.003_ mod± 0.001_ fit.
            
            
            
            
             For what concerns  <r_K^±^2>, the numbers reported in Table <ref>,
             show that our evaluation is at 1 σ from the NA7 measurement <cit.>
             and coincides with  those from <cit.>. As for the
             ChPT prediction
             from  Bijnens and Talavera <cit.>,  it differs from the value returned by BHLS_2
             at the 1 σ level.
             There is no reported evaluation of the curvature term of the kaon form factor to which
             one could compare our own evaluation in Equations (<ref>).
            
             Therefore,  concerning the low energy behavior of the pion and kaon
             form factors, the overall picture is  a fair agreement of  BHLS_2 with
             the corresponding experimental measurements  as well as with the existing
             LQCD and ChPT evaluations.
            
            
            
            § THE ISOSPIN 1 P-WAVE ΠΠ PHASE SHIFT
             Another relevant piece of information is the I=1  dipion
            P-wave phase shift δ_11. This physics has not been included within our HLS
            frameworks[A relative energy scale calibration issue between phase shift 
            and annihilation data 
             may have been identified <cit.>, probably amplified by the steep rise of
             the phase in the interesting region, combined with  a 20 MeV binning size
             for what concerns the CERN-Munich data  <cit.>. ]; 
             nevertheless, a comparison of the P-wave ππ phase shift
            extracted from BHLS/BHLS_2 with other sources provides an additional cross check.
             
            
            
            Both panels of Figure <ref> display
            the experimental data collected by the CERN-Munich Collaboration <cit.>,
            the phase shift derived using the Roy Equations <cit.>  and the phase shift
            derived in  <cit.> in the context of a common fit of the e^+e^- π^+ π^-
            cross-section and of the dipion spectrum in the τ decay[Plots
            including other data  <cit.> and the ChPT predictions <cit.> can be found in Figure 5 of
            <cit.>; they are not shown here in order to avoid overburdening Figure <ref>.
            ]. The phase shift coming from BHLS_2  is  displayed in the top panel
            while the  phase shift from the previous BHLS <cit.> is displayed in the bottom
            panel. Both HLS based phase shifts displayed correspond to the phase of the pion form
            factor amplitude amputated from the   and ϕ terms (see Equation (<ref>))
            to meet the usual definition of the I=1P-wave ππ phase shift <cit.>. 
            The full ( i.e. non amputated) BHLS/BHLS_2 pion form factor phases are visually
            identical to the ones displayed in Figure <ref> except for additional
            tiny blips in the   and ϕ energy regions as can also be seen in Figure 5
            of <cit.>. 
            
            
            In the case of BHLS_2 (red curve in the top panel of Figure <ref>), the 
            phase shift superimposes almost exactly on the data points <cit.> and onto the JS11 
            curve up to the ϕ mass region. For this comparison
            with the external theoretical expectations <cit.>, one should cancel out
            the  and ϕ terms from  the pion form factor F_π^e(s) given in 
            Equation (<ref>) as, indeed, these expectations do not carry such signals; one
            should remind that,  within the BHLS/BHLS_2 frameworks, these signals are generated 
            via the mixing angles α(s) and β(s) (see Equations (<ref>)) which
            are proportional to the charged and neutral kaon loop difference and this vanishes in the
            isospin symmetry limit of equal kaon masses.
            
            
            
            This argument might not fully apply when comparing the BHLS_2 phase to the 
            experimental data <cit.>. However, the  CERN–Munich 
             Collaboration reported <cit.> a tiny ρ^0- mixing effect – actually 
             unsignificant (see Table 1 in ths Reference) – in the <cit.> phase shift data.
             Nevertheless, let us quickly consider this aspect of the 
             experimental phase shift for completeness.
            
            For a well-grounded comparison between the BHLS_2 phase and the  CERN–Munich data, 
            one should rather use the BHLS_2 phase smeared out over the experimental (20 MeV) 
             bin width. The purpose of Figure <ref> is to display 
            the smeared BHLS_2 phase together with the CERN-Munich data <cit.>.
            In its leftmost panel the smeared BHLS_2 phase plotted is those of the full pion form factor
            F_π^e(s) just as given in Equation (<ref>); in the rightmost panel, the smeared
            BHLS_2 phase is those of the (/ϕ)-amputated F_π^e(s). In both panels, the
            CERN-Munich data points are also plotted but their symbols generally vanish underneath those for
            the smeared BHLS_2 data. One clearly observes a good agreement between BHLS_2 and the
            CERN-Munich data up to the ϕ mass region. Comparing the phase shift of the full and 
            amputated BHLS_2 pion form factors as displayed in Figure <ref> clearly
            shows that the effects of the  and ϕ signals on  δ_11 is quite 
            marginal as expected from the  <cit.> analysis.
            
            
            
            
            In a large energy interval
            (about ≃ 200 MeV wide) centered on the  mass, the 3 curves for the BHLS/BHLS_2, Roy and JS11 phase
            shifts overlap within the thickness of the curves and, so, it is impossible 
            to distinguish any one from each other. Then, 
            both JS11 and BHLS_2 start to depart from the data increasingly with increasing energies. 
            Put together with the remarks already stated
            about the behavior of the dikaon cross sections slightly above the ϕ mass, this gives a hint
            for an  onset of higher mass vector mesons presently not accounted for in the HLS framework.
            
            Thus, the agreement between the  BHLS/BHLS_2 predictions on the one hand and
             the theoretical expectations <cit.> and the CERN-Munich data on the 
             other hand, looks quite satisfactory[As illustrated in Figure 2 of  <cit.> 
             the solution of the Roy equation requires a normalization of the energy scale at
            about 1 GeV.  In the BHLS approach, one could adjust g_ρππ
            to get a better agreement for the high energy tail of the Roy solution.]. 
            
            One should also note that the bottom panel in Figure <ref> exhibits
            a fair agreement with the datat much above the validity range of the BHLS model (the
            ϕ mass region); this can be traced back to having allowed the pion loop in the
            V-γ transition amplitude F_ρ^e and in the ρ self-mass amplitude  to carry
            different subtraction polynomials (see Subsection <ref> above). This clearly accounts for
            the onset of the high mass vector mesons as the agreement remains valid up to 1.3 GeV and somewhat above.
            
            
            
            
            § LAGRANGIAN PARAMETER VALUES IN BHLS AND BHLS_2
             Table <ref> reports on the Lagrangian parameter values[The  full list of
            parameter values and uncertainties, especially   the subtraction polynomial coefficients, are not given;
            they can be provided upon request.]
            derived when running BHLS  with updated kaon data (first data column). Running the basic variant BS of BHLS_2
            with τ data discarded gives parameters values displayed in the second data column, while the third data column
            reports on the  numerical results from the BHLS_2 Reference Solution (RS) using the largest set of consistent
            data samples.
            
            Clearly, the  Lagrangian parameters common to the various modelings cover a wide range of
             values; this does not prevent  the fit properties (reminded in the last two data lines) to
             always exhibit a fairly good description of the data.
            
            
            
             §.§ The (z_V ÷ a) and (Σ_V ÷ a) Anticorrelations within BHLS_2
            
            The value obtained for z_V from fitting  with the
            RS  variant of BHLS_2 may raise  questions – see Section <ref>– as z_V-1
            is numerically large[ Even if only a 2.4 ϵ effect, referring to the numerical
            estimate for generic ϵ's given in Footnote <ref>.],
            even larger than z_A-1. For
            this purpose, it has been worthwhile to explore the RS variant behavior in fits performed
            discarding the 3 π final state data; this has revealed a numerical anticorrelation between the  HLS
            parameter a and the BKY breaking parameter z_V, undetectable in the parameter error
            correlation matrix  as <δ a  δ z_V> ≃ -0.24. Likewise, but of less concern,
            one observes a similar
            numerical anticorrelation between a and Σ_V as <δ a  δΣ_V>≃ -0.18.
            Several solutions with comparable total χ^2 were found and the  most striking difference among them
            just concerns the values for a, z_V and Σ_V. This is well illustrated by the three following
            parameter sets :
            {[    a = 2.340 ± 0.003 ,     z_V=1.122 ± 0.003,  Σ_V = -0.272 ± 0.001, χ^2/N_ pts= 981.2/1079;    a = 2.083 ± 0.003 ,     z_V=1.266 ± 0.005,  Σ_V = -0.182 ± 0.001, χ^2/N_ pts= 981.5/1079;    a = 1.606 ± 0.001 ,     z_V=1.653 ± 0.002,  Σ_V = +0.068 ± 0.002, χ^2/N_ pts= 983.5/1079 ] .
            
            
             which exhibit  clear anticorrelations between the numerical values for a and z_V on the one hand and
             for  a and Σ_V on the other hand.  Despite the large variations of each of a, z_V and
             Σ_V,  the ranges covered by the central values for a z_V and a (1+Σ_V) 
             are much narrower than those for  a, z_V and Σ_V; indeed, one obtains 
             a z_V = (2.63 ÷ 2.66) and a (1+Σ_V) = (1.70 ÷ 1.72),   i.e.
             both products exhibit spreads at not more than the percent level. 
            
             In order to track back the source of these anticorrelations, it is worth considering the
             ρ^0 π^+ π^-  and ϕ K  K couplings. These are given by (see Appendix <ref>) :
              [ g^ϕ_K  K = ∓[a z_V]  g√(2)/4 z_A(1 ±Δ_A/4) (1+ξ_0)   and                 g^ρ^0_π^+ π^- =[a (1+Σ_V) ] g/2 (1+ξ_3). ] which exhibit the dependences upon [a z_V] and [a (1+Σ_V) ].
            
            On the other hand, given the measured cross sections for e^+ e^- π^+ π^- and 
            e^+ e^-  K  K, the observables  
             Γ(ρ^0   e^+ e^-) ×Γ(ρ^0 π^+ π^-) and 
             Γ(ϕ  e^+ e^-)×Γ(ϕ   K  K) should be weakly model/parameter
             dependent within  successful (global) fits. Substantially, these products are 
             the squares of 
              H_ρ=| F^e_ρ(m_ρ^0^2)  g^ρ^0_π^+ π^-| and
             H_ϕ=|F^e_ϕ(m_ϕ^2)   g^ϕ_K  K|,
              the F^e_V (s)'s  being given by Equation (<ref>); these quantities  stretch
              ranges of the order ± 3%
              only :H_ρ=0.753 ÷ 0.803 and H_ϕ=0.124 ÷ 0.136.
              So, as the fits reported by Equations (<ref>) are successful, 
              they all reproduce successfully the resonance peak values  
            for the e^+ e^- π^+ π^- and  e^+ e^-  K  K cross sections, 
            which implies  that the products  a z_V and a (1+Σ_V) weakly vary within different
            fits, even if each of these parameters scans a wider interval[ The other parameters
            involved in the expressions for g^ϕ_K  K and g^ρ^0_π^+ π^-– i.e.g,Δ_A, z_A and the ξ's  –  undergo variations at the percent level or
            better.].
             
            
             
            One has also run the BHLS_2 RS variant by fixing z_V=1 and found the results shown in
            the last data column of Table <ref>. The 6 units difference between
            the minimum  χ^2 obtained fixing z_V=1 and the best  fit
            leaving z_V free is almost entirely concentrated in χ^2(e^+e^-  3 π)
            which becomes 146.3 (compared
            to 140.4 for 158 data points). Running also  the BHLS_2 RS variant  imposing via  minuit, a lower bound a≥ 2 returns :
            
            [         a = 2.0 ± 0.001  ,      z_V=1.330 ± 0.001   ,       Σ_V=-0.128 ± 0.001 , χ^2/N_ pts= 1133.2/1237  .;                            ]
            
            
            
            
             This fit also   returns  a z_V ≃ 2.66 and  a (1+Σ_V) ≃ 1.75
             and  shows fairly good quality (93.6% probability).  Moreover, all solutions
            return consistent  information about the physics observables as will be exemplified below.
            Therefore, the relatively large
              value for z_V-1 obtained in the best fit with the (unconstrained) RS variant of BHLS_2 is not
              a real issue. One should also note that the multiplicity of solutions
            reflected by   Equations (<ref>) reduces to a single one once the 3-pion data are accounted for.
            
            
             §.§ A Remark on the HLS Parameter a
            
            Correlated with the remark on z_V, one should note the value for the standard HLS model
            parameter a coming out of the best RS fit : 1.590 ± 0.001. This is much smaller
            than any of its previous determinations. One should also note the unexpected value for 
            a naturally returned by 
            the BS variant fit : a=1.970 ±  0.002, very close to 2   (see Table <ref>).
            Therefore, one may conclude that, within the BS variant of BHLS_2,  breaking parameters
            supply  the needed γπ^+π^-, γ K^+K^- and γ K_LK_S couplings 
            which, otherwise, would vanish when the a=2 constraint is imposed.
            
            
            
             §.§ The O(p^4) HLS Lagrangian : Effects of the Parameter z_3
             The present work is the first one where a manifest and
            thorough use of the 
            z_3 term from the O(p^4) HLS Lagrangian is performed; Appendix <ref>
            displays its full expression in the BHLS_2 context in terms of the ideal combinations
            of the vector fields. Actually, the existence of such a term
             was implicitly accounted for in the former BHLS <cit.> by disconnecting
            the subtraction polynomial for the V transitions from the subtraction
            polynomials accompanying the pion and kaon loops occuring elsewhere  – notably
            inside the ρ propagators.
            
            Within the BHLS_2 framework, the loop corrections to the V transition amplitudes
            and to the ρ^0 self-mass are tightly related[Compare Subsection 
            <ref> and Appendix <ref>.] and additional terms
             of the form c_V z_3 s come into play – see Equations (<ref>).
             
             Table <ref> shows that z_3 is small but significant.  As the V
             transition amplitudes, where solely z_3 occurs, affect all the cross sections
             considered,  the most appropriate way to examine its effects is to compare fits 
             where z_3 is let  free and fits performed by fixing z_3=0.
             
             The   overall picture coming from these fits is  reported just below. To be concise,
             one denotes here by BHLS_2(z_3) and BHLS_2(z_3=0) the RS based  fitting framework
             running with resp. z_3 free and z_3=0.
              
              *  When all data and channels considered in this paper  are used (the 3–pion channel 
             is discarded, for simplicity),
             one gets χ^2(BHLS_2(z_3)/N_data =982/1079 (93.1% probability),
             whereas χ^2(BHLS_2(z_3=0)/N_data=1049/1079 (50.6%  probability),
             a significant degradation.
             
             A careful comparison  of the fit results obtained in both cases
              shows that the partial χ^2 for the e^+e^- annihilation channels 
              (π^+ π^-, K^+ K^-, K_L K_S, π^0, η
             and even the τ dipion spectra) are comparable[ BHLS_2(z_3) 
             remains, nevertheless, favored against BHLS_2(z_3=0) in these channels.]. 
             
             
              *  The bulk of the degradation when going from BHLS_2(z_3) to BHLS_2(z_3=0)
             is located in the account by the latter of the spacelike sector where the NA7 data 
             impose accurate constraints, hardly absorbed by BHLS_2(z_3=0).
             
             Indeed, for the spacelike experimental data, the fit performed with BHLS_2(z_3=0)
             returns  χ^2(spacelike π)/N_data =78/59
             and χ^2(spacelike K)/N_data =29/25, whereas the corresponding numbers
             in the  (standard) BHLS_2(z_3) framework
             are resp. 56/59 and 18/25, in perfect agreement
             with the information already reported in Table <ref>.
             
             To fully compare with the information reported in Table  <ref>, a  
             BHLS_2(z_3=0) fit, based on only the timelike data ( i.e. 
            discarding the spacelike data), is useful as it provides  
            the  corresponding BHLS_2(z_3=0) predictions for the spacelike sector.
             
            This BHLS_2(z_3=0)  global fit has been performed and displays an improved 
            fit quality, as χ^2(BHLS_2(z_3=0)/N_data =939/995 
            (72.6% probability), nevertheless lower than those for BHLS_2(z_3). The χ^2 distances 
            of the π and K spacelike data to the corresponding BHLS_2(z_3=0) predictions  
            are χ^2(spacelike π)/N_data =85/59 and
              χ^2(spacelike K)/N_data =30/25.  
              
              These two numbers  should
               be compared with the first  data column in Table <ref>. 
               So, the nice feature of BHLS_2(z_3) which predicts accurately the spacelike data
               is lost with BHLS_2(z_3=0).
            
            One should note that, when z_3=0 is imposed, the effects accounted for by z_3  0
            are absorbed by the vector mesons propagators[For instance g jumps from ≃ 7
            to ≃ 10 and the pion and kaon loop subtraction polynomials are severely modified.  ];
            this is what degrades the continuation
            of the form factors to the spacelike region and, fitting the spacelike data, does not solve 
            the issue.
            
            Stated otherwise, a non-zero z_3 is the best way to treat simultaneously the spacelike
            and timelike data with high accuracy and comfortable fit probabilities. This may be considered
            as an evidence pointing toward the relevance of the O(p^4) terms; a range for z_3
            can be proposed  z_3 = -([0.37÷ 0.42] ± 0.003_fit)× 10^-2 
            for its  first experimental determination.
            
            
            
             §.§ Other Physics Results
            
            As for other model parameters, focusing on BHLS and the RS variant of BHLS_2 :
            
              *  One may wonder that the value for the vector coupling g increases by more than 20%  when going from BHLS to
            BHLS_2, keeping good fit qualities.
            However, the relevant piece of information, closer to observables, being the coupling
            g^ρ^0_ππ =ag(1+Σ_V)(1+ξ_3)/2, it is interesting to compare BHLS to the RS
            variant of BHLS_2 : BHLS returns g^ρ^0_ππ = 6.156 ± 0.005
            while  the full BHLS_2 fits gets g^ρ^0_ππ = 6.06 ± 0.013 in reasonable
            correspondence  with each other.
            
              *  The parameters associated with the anomalous sectors show that[In the present work,
             one
            assumes c_3=c_4.](c_3+c_4)/2 yields comparable values  with BHL and with BHLS_2. However, for
            c_1-c_2, the difference is large; in this case, however, one should remind that the
            so-called "model B" variant of BHLS  excludes from the fit the 3 π data in the ϕ mass
            region while the whole 3 π spectra are considered within BHLS_2 and are well
            fitted.
            
            
              *  The values for z_A are also interesting as they are related with f_K. Indeed, one knows
            that <cit.> :
            
                r=f_K^±/f_π^± = √(z_A) [1+ Δ_A/4]  .
            
            is still valid within the BHLS_2 framework. For this ratio, BHLS returns r=1.257 ± 0.003_fit, 
            while the RS variant of BHLS_2  finds r=1.247 ± 0.003_fit. 
            The systematic uncertainty has been estimated to ± 0.020 by varying the running
            conditions of the BHLS_2 fitting 
            code[We made additional runs  fixing z_V=1 (see Table <ref>), 
            including  the CMD-3 dikaon data and, finally, varying also the running of α_em.
            ].
            
             This result should be compared  with the
            Review of Particle Properties <cit.> which recommends r=1.197 ± 0.002 ± 0.006± 0.001 and
            with the LQCD determination for N_f=2+1+1 flavors <cit.>r=1.195 ± 0.005.
            
            
              *  The additional parameters introduced by the CD breaking and the Primordial Mixing are all found
            at the few percent level. The difference for ξ_0 and ξ_3 between the (full) RS variant 
            and the BS variant should be noted; it may have to be revisited later on, when a
            complementary mechanism may allow one to absorb satisfactorily the τ data within the 
            BS global fit procedure. This may also influence the 3 π sector.
            
            
            § THE MUON LO-HVP EVALUATION
              BHLS_2 (or BHLS) fits provide the contributions of
            the various hadronic channels[All along  this Section, one may use
            interchangeably a_μ^ HVP-LO or, preferably,  a_μ to lighten writing; the context will be
            unambiguous.] it covers
             to  a_μ(√(s)≤ 1.05   GeV). An overview of the main
            results, obtained in different running conditions is given in Table <ref>.
            These have been derived by submitting to fit most of the
            data samples commented in Section <ref> which define a set of
            consistent samples covering the six annihilations considered;
            altogether, these 6 channels and the 10 radiative meson decay partial widths represent 1237 data points.
            
            
            
             §.§ Comparison of BHLS_2 Estimates for  a_μ^ HVP-LO with Others 
              Without going lengthily into detailed comparisons, it looks worthwhile to
            provide some elements. For instance, one can extract from Figure 3 in <cit.>,  values
            for G=a_μ^ππ(0.630 <√(s)<0.958  GeV );  our updated
            estimate for  G is  G( BHLS)= (355.59 ± 0.58)× 10^-10, very close to its former
            evaluation (see Fig. 3 in <cit.>) while the numerical integration of the NSK+KLOE
            data was estimated G( exp)= (356.67 ± 1.69)× 10^-10; BHLS_2
            gives[As the (0.630 <√(s)<0.958) GeV  interval is relatively far
            from both the threshold and the ϕ mass regions, one
            could indeed expect similar results from BHLS and BHLS_2. ]G( BHLS_2)= (356.50 ± 0.75)× 10^-10.
            
            In a quite different theoretical context, the recent analysis <cit.>, provides the estimate
            G=a_μ^ππ(√(s) <0.630  GeV )=(132.5  ± 0.4_ stat± 1.0_ syst)× 10^-10 while the
            BHLS gives G( BHLS)=a_μ^ππ(√(s)<0.630 
            GeV)=(130.80 ± 0.22_ fit)× 10^-10.
            In the same region
            one also gets G( BHLS_2)=a_μ^ππ(√(s)<0.630 
            GeV)=(130.03 ± 0.18_ fit
            +[ ^+0.25_-0.15]_ syst)× 10^-10.
            
            Reference  <cit.> also provides
            H=a_μ^ππ(√(s)<1.0  GeV)= (494.7 ± 1.5_ stat±
            2.0_ syst)× 10^-10
            while the updated BHLS gives H( BHLS)=(490.65± 0.84_ fit)× 10^-10 
            and BHLS_2 leads to
            H( BHLS_2)=(490.75± 0.90_ fit+[ ^+0.31_-0.26]_ syst))× 10^-10. The method in
             <cit.> and ours are different, but more importantly for this energy range, the
              dipion BaBar data included in the study <cit.> and discarded
             in ours, certainly contribute to the difference.
            
            One may also compare the BHLS_2 results with others for the contribution to a_μ^ HVP-LO(π^+ π^-) from
            the energy region (0.6 ≤√(s)≤ 0.9)  GeV.  KNT18 <cit.> gets (369.41 ± 1.32)× 10^-10 
            and  CHS <cit.>(369.8 ±  1.3_fit± 1.3_syst)× 10^-10 using 
            BaBar data at low energies. BHLS2 running
            with KLOE10/12, CMD-2, SND, BESIII and Cleo-c data in the π^+ π^- channel yields
            (367.81 ± 0.77)× 10^-10; replacing the KLOE10/12 data by the KLOE85 combination
            <cit.>, one gets (367.61 ± 0.64)× 10^-10. So the
            observed effect from BaBar in this energy region amounts to ≃ 2× 10^-10. One also
            observes the gain on accuracy by using BHLS_2 : A factor of 2 compared to  <cit.>,
            a factor of 3  compared to  <cit.>.
            
            
            
            
             §.§ Estimates for  a_μ^ HVP-LO(√(s)≤ 1.05  GeV)
              In order to get the full contribution to a_μ^ HVP-LO of the energy
            region up to 1.05 GeV,
            one should complement the  contribution from the channels covered by BHLS/BHLS_2 by those
            of the channels presently outside this framework (4 π, 2 πη, …). The most recent
            evaluation of this quantity, estimated by direct numerical integration of data  is :
            
                10^10× a_μ( non HLS, √(s)≤ 1.05  GeV) = 1.28 ± 0.17
            
            
            In view of model effect studies, it has been found motivated to update and correct
            the data samples submitted to the BHLS fit  in the way stated in Section <ref>. One
            has thus derived  the updated BHLS evaluations given in the first data column of
            Table <ref>; this supersedes the BHLS results formerly given in Table 4 of
            <cit.>. Differences between BHLS and the RS BHLS_2 variant for some channel contributions 
            are observed and,  in fine,
            a_μ( HLS, √(s)≤ 1.05  GeV) only increases by  ≃ 1 × 10^-10 while
            its uncertainty slightly decreases (1.03  0.92, in units of 10^-10).
            An evaluation of  possible additional systematics, specific of the BHLS modeling, has
            been performed and summarizes by <cit.> :
            
            10^10×δ_ syst a_μ( BHLS) =  [ ^+1.3_-0.6]_ϕ +  [ ^+1.4_-0.0]_VNSB
            + [ ^+0.0_-0.9]_τ   .
            
            
            Each piece shown here was found to  act by shifting the central value for a_μ  rather than enlarging its
            uncertainty. The first term refers to the uncertainty in the treatment of
            the ϕ mass region in the 3-pion spectra, the second term takes into account
            the uncertainty on the ππ threshold  behavior within BHLS (see Subsection <ref>
            above and, especially, Figure
            <ref>); finally, the third piece reflects differences observed by running
            the BHLS fit procedure <cit.> with the τ data included
            and excluded.
            
            As already stated, the BHLS_2 modeling has been motivated by the aim to cancel out the
            first two sources of uncertainties reported in Equation (<ref>) by a better
            understanding of the s=0 region and a full account of the 3-pion data
            up to the ϕ region.
            
            For this purpose, the three data columns in Table <ref> referring to BHLS_2 show
            that :
            
              *  The evaluation for the various contributions to a_μ(√(s)≤ 1.05  GeV)
            derived using the BS and both RS variants of   BHLS_2 carries differences at a
            few  10^-12 level.
            
              *  Using  the τ data has some  effect. Indeed, Table <ref>
            shows that running the RS variant of BHLS_2 including the τ data improves
            the uncertainty in the dipion channel by ≃ 10 % while leaving the central 
            value almost unchanged. This illustrates that the  e^+e^- and τ 
            data are not conflicting within the RS variant of   BHLS_2.
            
            Another property exhibited in Table <ref> should be noted :
            
              *  Whereas the sum of the various contributions to a_μ derived using
            several variants of the BHLS_2 fit almost coincides with the estimate derived
            by direct integration of the data, the π^+ π^- contributions differ by 3.72
            in units of 10^-10, a significant difference.  
            
            In the direct integration method, the normalization of each sample is the nominal one 
            and  all  uncertainties participate to the definition of its weight in the combination
            with the other samples involved. Instead, in the global fit method, the absolute
            normalization uncertainty is treated specifically and in 
            consistency for all the data samples involved in all the annihilation channels considered
            (ππ, K K, πππ and (π^0/η) γ). It has beeen
            shown <cit.> that the iterative method which underlies the HLS global fit procedures
            gives unbiased results [This is well reflected, in the present work, by the nice
            matching observed between the predicted pion form factor in the spacelike region and the NA7 and ETMC
            data.]. 
            
            Therefore, that values of definite contributions to a_μ 
            may significantly depend on the method used to derive them is not totally
            unexpected. That the
            sum of them may almost coincide  is, however, accidental[ 
            In this respect, it is worth noting  that a quite 
            similar behavior is observed in the comparison of the KNT18 and DHMZ17 data reported in 
            <cit.>. Indeed, Table 5 herein shows that the estimated dipion contributions (based on
            the same data) differ by 3.40 × 10^-10, whereas the total sums giving the  
            HVP-LO's differ by only
            0.2× 10^-10.].
            
            
            The data update having been performed, BHLS and BHLS_2 have been run on the same
            data and, thus, the differences between their respective evaluations reflect modeling effects.
            One can also compare these differences to
            the  numbers listed in Equation (<ref>).
            
            As shown by the
            various topics examined in Section <ref>, one can legitimately consider
            the  π^+ π^- threshold region accurately treated.
            Therefore, the difference
            494.08-493.73=0.35 between the BHLS_2 and BHLS estimates for
            a_μ(π^+ π^-, √(s)≤ 1.05  GeV)  is a fairly well 
            motivated evaluation
            of the unaccounted for nonet symmetry breaking effects within the BHLS framework;
            it is therefore justified  to replace the  former +1.4 estimate
            in Equation (<ref>) by +0.35.
            
            Table <ref> shows that all variants of BHLS_2 take
            well into account the 3 π annihilation data over the whole energy domain up to the
            ϕ mass region, the fit quality of the RS variant being noticeable.
            The difference
            between BHLS_2 and BHLS for this channel ≃ 0.5 × 10^-10 is in the range
            expected from the estimate <cit.> reported by the first term in Equation (<ref>).
            Nevertheless, the 0.06 × 10^-10  difference for the 3 π entry
            between the RS variant of BHLS_2 runnings
            including and excluding the τ data  might  be considered as
            additional  systematics[There might be a numerical effect
            reflecting the  physics correlation between the π^+ π^- π^0 annihilation
            channel and the τ spectrum due to their common ρ^± exchanges
            taking place in their intermediate states. ].
            
            It is interesting to compare the outcomes from fits to the values derived by the direct integration
            of the same data sample set treated in the BHLS_2 framework. One may note that the sum of all HLS channel
            contributions differs from the fit expectations by only
            ≃ 0.6 × 10^-10; however, this hides the fact that the dipion
            channel is fitted
            3.7 × 10^-10 smaller and the 3-pion channel is fitted 1.70 × 10^-10 larger.
            
            
            In summary :
            
            10^10× a_μ( HLS, √(s)≤ 1.05  GeV) = 571.98 ± 1.20_fit
            is the most conservative  final  BHLS_2 answer, where the uncertainty collects statistical and systematic
            errors as reported by the various experiments. The remarks just stated concerning the
            3-pion channel
            allows to keep an additional uncertainty  of -0.06 × 10^-10,
            playing as a bias.
            
            
            The issue is to examine whether more significant additional sources of systematics
            are not at work. For this purpose, the consequence on a_μ( HLS) of the
            various fits referred to in the previous Section carry a relevant piece of information,
            as also for the fits underlying the results given in Equations (<ref>)
            and (<ref>)  or those
            displayed in the last data column in Table <ref>. These fits, covering a wide
            range of parameter values, are relevant for systematics estimates.
            
             
            
            
            The runs leaving aside the 3-pion data briefly reported in Equations (<ref>), show
            that the  dominant ππ contribution can be as low as
            (493.80± 0.89)× 10^-10,  i.e.0.30 units of 10^-10
            smaller than the corresponding reference entry in Table <ref>. Nevertheless, the middle
            entry line in Equations (<ref>) is very close to Equations (<ref>), where the
            3-pion data have been included and exhibits negligible differences with respect to our
            reference in Table <ref>. Finally, the fit corresponding to the last data column
            in Table <ref> has also been analyzed and returns  a value for
            a_μ(π^+ π^-, √(s)≤ 1.05  GeV) = (494.32 ± 0.90)× 10^-10.
            This may indicate another additional systematic
            able to shift the central value by, at most,  0.22 × 10^-10
            affecting our estimate    for   a_μ^ HVP-LO(√(s)≤ 1.05  GeV) given in Equation (<ref>).
            
            The possible sources of additional systematics just emphasized can be considered 
            as a model uncertainty bounded by (+0.22,-0.36) × 10^-10. On the other hand, running the 
            fitting code using the available parametrizations for α_em(s), returns differences for
            the estimates for a_μ( HLS) of ≃ 0.4 × 10^-10, mostly located in
            in the ϕ mass region.
            
            Finally, the question of using the CMD-3 data can be raised. We have run our code
            using the CMD-3 data <cit.> for both K K decay
            channels, discarding the kaon data from BaBar, CMD-2 and SND and using only the diagonal part
            of their error covariance matrices (see Section <ref> above).
            In total, the change for
            a_μ(√(s)≤ 1.05  GeV) is noticeable (572.52 ± 1.12)× 10^-10,
            0.54× 10^-10 larger than our reference in Equation (<ref>), essentially coming from
             the K^+ K^- and K^0 K^0 channels. On the other hand, combining CMD-3,
             BaBar, (corrected) CMD-2 <cit.> and SND, leads to
            a_μ(√(s)≤ 1.05   GeV)=(572.52 ± 1.04)× 10^-10. This leads us to complete
            our estimate Equation (<ref>) :
            
            10^10× a_μ( BHLS_2, √(s)≤ 1.05  GeV) =
            571.98 ± 1.20_ fit +
            [^+1.16_-0.75]_ syst  ,
            
            where the additional systematics come from model variations, α_em
            parametrizations, and data sample consistency issues. The model uncertainties include the
            marginal effect attributable to  the τ data, as this could reflect some
            (marginal) isospin breaking shortcoming.
            We have chosen conservatively as reference  the RS variant fit excluding the τ data.
            One should also note that the tension between dikaon data samples introduces
            non-negligible systematics which  contribute a bias. This may have to be revisited
            when new data will arise.
            
            
            Taking into account the data upgrade, the corresponding quantity for BHLS is :
            
            10^10× a_μ( BHLS, √(s)≤ 1.05  GeV) = 570.97
            ± 0.92_ fit +
            [^+1.16_-0.75]_ syst+[^+0.0_-0.9]_τ  ,
            
            using  the present findings to go beyond the former systematics evaluation
            reminded in Equation (<ref>). The various additional uncertainties just 
            estimated  represent conservative upper bounds.
             
            
            
             §.§ Dispersive vs LQCD Methods : Additional Systematics
            
            
            Experimental e^+e^-→ hadrons data are dressed by photon
            radiative corrections and disentangling strong interaction from the
            electromagnetic effects is a non-trivial task. What is ideally needed
            as input to the dispersion integrals is the one-particle
            irreducible (1PI) hadronic part. It is given in QCD by the correlator
            between two hadronic currents, and this is what is primarily
            calculated as the LO part of order O(α^2) in lattice QCD. We
            remind that the bare cross-section undressed from photon radiation
            effects is not an observable by itself. It requires some theory input,
            primarily the photon radiative corrections, for its extraction from
            the physical dressed data. Mandatorily, one has to  drop out the photon
            vacuum polarization effects to get the undressed bare
            cross-section. By convention, one also includes final state radiation (FSR)
            in the dipion channel by adopting scalar QED (sQED) for its
            calculation. Similarly, the radiative decay channels π^0γ and ηγ 
            are included in our evaluation. Like FSR, they are related to 
            light-by-light scattering insertions.
            
             
            
            However, there are other important QED effects like the
            ρ^0-γ interference which is substantial in the dipion
            channel because this QED effect is magnified by the resonance
            enhancement in the ρ region. In the standard dispersive approach,
            ρ^0-γ mixing effects – inherent in the e^+e^- data – are
            treated as part of the LO-HVP contribution and are included as well in
            the calculations of the higher-order contributions. An updated
            e^+e^- data-based evaluation of the LO and the NLO parts yields the
            results shown in Table <ref>[For the present
            exercise, the BaBar dipion data have been taken into account].
            
            Here we have to remind that, in contrast to the usual dispersive
            estimations, LQCD calculations are based on the hadronic current
            correlator such that ρ^0-γ mixing is not included in the
            leading order LQCD results (unless it is added by hand from
            phenomenological analyses, if not calculated separately).  
            
            So the
            question of how to disentangle the QCD effects from QED ones, has to be
            reconsidered under this aspect. One possibility is to apply effective
            field theory methods like the resonance Lagrangian approach, the other
            is lattice QCD, the only method which can answer the question from
            first principles.  In fact, to disentangle better QCD from QED effects
            in the e^+e^- data-based dispersive approach, one has to subtract
            the ρ^0-γ interference from the e^+e^-→π^+π^-
            data[The similar ω-γ and ϕ-γ effects
            are much smaller and can be neglected.], expecting to get the effect
            back as a part of the higher-order contributions.
            
            
            Here  is the place to remind that the ρ^0-γ mixing is
            responsible for the discrepancy between the isospin rotated isospin
            breaking corrected τ spectral-function and the e^+e^-
            spectral-function <cit.>;
            this has been shown in a straightforward “VMD-II+sQED” simple
            effective resonance Lagrangian model <cit.> as well as in the
            BHLS approach <cit.> as utilized in the present
            analysis. By convention, to bring τ data closer to e^+e^--data,
            it looked natural to correct the τ-data to conform with the
            commonly adopted dispersive approach based on the e^+e^- data. 
            
            In
            fact, the measured τ-decay dipion spectra are devoid of photon
            interference effects. According to what we just have been arguing,
             one should actually not apply this mixing correction to the
            τ-decay data, as it has been done within our BHLS global fit
            framework or in the much simpler model applied in <cit.>. One
             should better subtract the ρ^0-γ mixing from the e^+e^- data,
            to get an object closer to a pure hadronic “blob” and then
            proceed as usual. As a possible subtlety, this requires to work
            in the “VMD-I” basis where the ρ^0-γ system is
            parametrized in terms of fields where a direct lepton-ρ coupling
            is absent. This is the basis  of the BHLS Lagrangian
            formulation. Being proportional to the hadronic production
            cross-section, the hadronic 1PI amplitude and the corresponding R(s)
            (see Eq. (<ref>))
            do not depend on rotations in the (γ,ρ^0) field space. 
            
            The
            field rotation invariance has been directly confirmed by comparing, on
            the one hand, the “VMD-II + sQED” modeling adopted
            in <cit.>, which uses the mass eigenstate basis where a
            direct ρ-lepton coupling is induced by the necessary field
            rotation <cit.> and,  on the other hand, our BHLS modeling
            where a direct ρ-lepton coupling is absent. One thus can apply the
            standard diagrammatic expansion of Fig. <ref>
            with the modified effective R(s)-function and
            obtains the results presented in
            Table <ref> :
            
            which corresponds to what, in principle, lattice QCD is doing. 
            
            The LQCD
            results displayed in Fig. <ref> are full HVP results, 
            which include ud
            light quark, strange and charm contributions including the connected
            and disconnected parts as well as isospin breaking contributions
            either calculated on the lattice or adopted from phenomenological
            estimates in terms of data. They do not include γ-ρ
            mixing to our knowledge, however.  In lattice QCD this effect is counted 
            as a NLO effect, which has to be evaluated separately[
            Contributions from FSR and radiative decays should have been added 
            in a consistent comparison with data-based dispersive results.  
            ].
            
            According to our Tables, the commonly adopted dispersive set up gets a
            shift[The systematic uncertainties in Tables <ref> and 
            <ref> are 100% correlated, except for the sub-dominant
            diagram c.] of δ a_μ [ππ]_ργ≃ +(3.10± 0.03)
            ×10^-10 at NLO,  i.e.  a 1 σ increase.
            See also the results presented in Table 3
            of <cit.> where a ρ^0-γ mixing shift
            of 2.74×10^-10 on the I=1 component at LO has been reported
            as a difference between LQCD data and the 
            phenomenological dispersive approach.The upgraded result 
            for  δ a_μ [ππ]_ργ
            is based on an actualized data compilation, with mixing correction 
            applied to the full vector form-factor and including NLO corrections. 
            
            More or less surprisingly, 
            the sum of LO+NLO differs in the two settings,  i.e. the mixing effect
            subtracted from the LO result does not show up as a higher-order
            effect of comparable magnitude. It is actually not surprising because the
            ρ^0-γ mixing effect, if included in the LO dispersion
            integral,  is larger than if taken into account in higher-order diagrams
            where the effect is weighted by the higher-order kernel 
            functions. 
            Our result based on correcting the e^+e^- di­pion channel in the range up
            to 1 GeV, where sQED estimates can be trusted, indirectly compares with
            the deviation (9.1 ± 5.0)× 10^-10 (1.8
            σ) between the τ-based hadronic contribution to a_μ and
            the e^+e^--based one reported in <cit.>. Our result does
            not include τ data while the τ based result of <cit.>
            includes τ data up to 1.8 GeV which are combined with e^+e^-
            data which still include ρ^0-γ mixing. That the results are
            increasing in both cases is as expected, but there is no reason why the
            numerical results should agree.
            
            
            Our suggestion to base the standard bookkeeping such as Figure <ref> 
            on a purer QCD γ-ρ undressed hadronic blob is motivated by the
            fact that lattice QCD does not include any QED effects if not taken
            into account as an extra contribution and also by the fact that
            τ-data are free of similar mixing effects, and it looks somewhat
            artificial to add effects to τ-decay spectral data that are
            absent in the corresponding measurements. Nevertheless, as already
            said, surprisingly we subtract a non-negligible contribution, which we
            do not get back in higher orders. It is important to note that the
            shift δ a_μ [ππ]_ργ is real as it derives from
            the manifest discrepancy observed between τ and e^+e^-
            spectral-function data in the corresponding time-like ππ
            production processes. As mentioned before, the γ-ρ
            interference is substantial because of the manifest resonance enhancement 
            in the time-like regime. Lattice QCD is taking advantage
            of the fact that, for the calculation of a_μ^ had, only the
            space-like VP-function is needed. However, there is no ρ-resonance 
            peak  in the Euclidean region which could enhance a
            photon-exchange interference effect; the ρ-resonance is
            completely smeared out and a resonance-enhancement cannot be localized
            there in any obvious way. Because of a model-independent clean separation of
            QCD is only possible in lattice QCD,  one has to wait for the relevant
            LQCD results to get a better understanding of QED-QCD interference effects. 
            
            One may conclude that  all  standard dispersive evaluations should 
            either be upgraded by adding +(3.10± 0.31) × 10^-10 to the 
            LO+NLO results[Concerning the model dependence of
            the prediction for the ρ^0-γ interference, we remind that this
            effect depends on only one extra parameter, the well known  
            leptonic width of the ρ^0 or, equivalently, on the coupling
            g_ρ ee, also well known experimentally. Since the effect is
            small, a few units in 10^-10, we can generously assign a very
            conservative 10% model uncertainty, without substantially increasing
            the overall uncertainties coming with the e^+e^- →π^+ π^-
            data.], or accounted for in the systematic errors by adding there
            the amount of this shift. For the time being, we have adopted this second solution.
            
            
            
             §.§ BHLS_2 Evaluation of  a_μ^ HVP-LO
             With the BHLS_2 and BHLS evaluations of a_μ^
            HVP-LO( HLS, √(s)≤ 1.05 GeV
            displayed in Equations (<ref>) and (<ref>) and the 
            information given in the left
            part of Table <ref>, one can derive the corresponding values
            for  the full  a_μ^ HVP-LO.
            One thus gets :
            {[                                              10^10× a_μ( BHLS_2) = 686.65 ± 3.01_ fit +
            [^+1.16_-0.75]_ syst + [^+3.10_-0.0]_ργ ,;                                               10^10× a_μ( BHLS )   = 685.64 ± 2.91_ fit +
            [^+1.16_-1.65]_ syst + [^+3.10_-0.0]_ργ ] .
            
            In the case of BHLS, due to the inclusion of τ channel within the global fit
            procedure, the corresponding systematics have been added linearly
            to the others. Indeed, as stated several times for BHLS and BHLS_2, the observed systematics rather
            come as biases and so, should be treated separately from the fit error which fully absorbs
            the reported statistical and systematics  from the various experiments encompassed inside
            each of our HLS frameworks. 
            
            The additional systematics absorb uncertainties which can be attributed
            to the model workings  and to the (marginal) tension observed between some of the selected
            experimental data samples. For clarity, the error coming from the ρ^0-γ
            interference effect discussed just above is given separately.
            
            Our evaluations for a_μ^ HVP-LO, given by Equations (<ref>), are shown in
            Figure <ref> together with some other recent results derived in different
            phenomenological contexts <cit.> or via LQCD data samples  with N_f=2+1+1. 
            As for their specific
            HLS modeling parts, the corresponding fit probabilities shown by the down-most entries in
            Table <ref> are always fairly good.
            
            
            
            
            
            
             As for modeling effects, the central values for BHLS and BHLS_2– running on the same set of data samples –
             are shifted apart by ≃ 1.5 × 10^-10; however, if one takes into account the magnitude
             of the additional systematics, Figure <ref> clearly indicates that this shift is not significant
             and so no noticeable modeling effect shows up. Moreover, Figure <ref> clearly shows that
             both HLS based estimates for a_μ^ HVP-LO are consistent with the evaluation derived by numerical
             integration of the data[See the top point in Figure <ref>, derived by
             excluding the BaBar dipion sample for a consistent comparison with the BHLS entries. ].
             
             Actually, the estimates for the muon HVP-LO derived by the various dispersive methods (the five data points 
             displayed in the bottom  part of Figure <ref>) are quite comparable with each other 
             as they use as input to evaluate  a_μ^HVP-LO
             the e^+e^- annihilation cross sections as derived from the data or as reconstructed from fit. Therefore, 
             the results thus provided do not split out the electromagnetic corrections; 
             in particular the γ - ρ^0  mixing effects is merged in the plotted systematics 
            only for the BHLS/BHLS_2 entries, not for the other dispersive results plotted.  
            
            
            
             §.§ BHLS_2 Evaluation of  Muon Anomalous Magnetic Moment
            
            Having at hand our estimate   Equation (<ref>), one can
            derive the BHLS_2 evaluation for a_μ=(g-2)/2, the muon anomalous magnetic
            moment. The left part of Table <ref> collects the pieces to add up
            with a_μ( BHLS_2) to get the full leading  order HVP, its value is shown
            as the top entry in the right-hand part of  Table <ref> where
            the other contributions are included.  The value for
            Δ a_μ=a_μ( exp)-a_μ( th) and its statistical significance are 
            the bottom entries there.
            Taking into account all additional systematics, one thus gets :
            
            10^10× a_μ= 11 659 175.96 ± 4.17_ th 
            +[^+1.16_-0.75]_ syst +[^+3.10_-0.0]_ργ
            and, taking into account the shifts produced by all systematics, one gets :
            
            28.88± 4.17 ≤
            10^10×Δ a_μ (= a_μ^exp.-a_μ^theor.)
            ≤ 33.89  ± 4.17  ,
            
            which means that the statistical significance of Δ a_μ is
            predicted greater than 3.82 σ. The shift just reported is partly due
            to the estimates of the various modeling effects involved, partly due to somewhat conflicting
            aspects of some data samples and parly due to the identified γ-ρ^0 mixing 
            effects which add up linearly. 
            
            
            Our estimate for a_μ given in Equation (<ref>) relies on a variant of the
             traditional estimate  <cit.> for the Light-by-Light (LBL) contribution
            ((10.34 ± 2.88)× 10^-10). Danilkin, Redmer and Vanderhaegen (DRV)
            have quite recently published a comprehensive analysis of all ingredients participating
            to the LBL term and proposed for it the much more precise value <cit.>(8.7 ± 1.3)× 10^-10 which leads to :
            
                10^10× a_μ^ DRV= 11 659 174.78 ± 3.28_ th
                +[^+1.16_-0.75]_ syst  +[^+3.10_-0.0]_ργ
            
            and to Δ a_μ^ DRV=(34.32 ± 7.10)× 10^-10. This makes
            the significance for Δ a_μ jumping from (3.8÷ 4.5) σ to (4.2÷ 4.7) σ.
            
            
            
            § SUMMARY AND CONCLUSIONS
             The present paper has addressed several topics covering modeling,
            phenomenology and the evaluation of the muon anomalous magnetic moment which puts
            forward, as well known, a hint for a physics effect beyond the Standard Model expectations.
            For the sake of clarity,  each of these topics is considered separately.
            
            
             §.§ From BHLS to BHLS_2
             The BHLS framework  <cit.>, although quite successful
            <cit.>, has been shown to exhibit a non-optimum behavior
            in limited and well identified kinematical regions : The threshold region
            and some aspects of the ϕ mass region; the former issue renders delicate the
            continuation of form factors across the chiral point, the latter has led
            to discard the ϕ mass region for 3-pion final state – the so-called Model B   <cit.>.
            
            In order to cure these diseases, one  has  initiated a new breaking
            procedure of the general HLS framework <cit.>.
            For this purpose, besides the BKY mechanism <cit.>, one
            has introduced a breaking mechanism taking place at the level of the covariant
            derivative (CD) itself,  i.e. inside the basic ingredient  of the HLS Model.
            This (seemingly novel) CD breaking mechanism allows naturally the ρ^0,  and
            ρ^± fields to carry different Higgs-Kibble masses. This takes its importance when
            going to loop corrections to the vector meson mass matrix <cit.> which
             generate  a dynamical mixing of the vector mesons  calling for an s-dependent
            redefinition of the vector meson fields to recover mass eigenstates.
            
            Thanks to its nonet symmetry breaking piece, the CD breaking scheme generates
            a non-zero  m_ρ^0- m_ difference which  naturally makes
            the three s-dependent mixing angles to vanish at the chiral point. In this way,
            the timelike to spacelike analytic continuations become smooth without any  ad hoc trick.
            
            Relying on BKY and the CD breaking, the Basic Solution (BS) is defined and represents
            the first variant of BHLS_2.
            Another mechanism can also be invoked which turns out to state that the neutral vector
            fields involved in physical processes are not their  ideal combinations but mixtures of
            these, constructed via a rotation. This Primordial Mixing (PM)
            mechanism[The Primordial Mixing
            of neutral vector mesons resembles the mixing scheme for the
            neutral pseudoscalar mesons developed in <cit.>.] of vector mesons
            assumes a 1+ O(ϵ) rotation. Working at first order in
            breaking parameters, the merging of the BKY and CD breakings  together
            with the PM mechanism define our Reference Solution (RS), the second
            variant of BHLS_2  covering the full HLS scope  more easily than the
            BS variant (as it stands presently).
            
            In both variants of BHLS_2, the L_V piece of the non-anomalous Lagrangian
            is modified  while its L_A piece remains
            unchanged[This might have to
            be revisited in order to improve the working of BS in the τ and π^+ π^- π^0
            sectors of the broken HLS Lagrangian.]  compared to the former BHLS <cit.>.
            
            
            
             §.§ The e^+e^- Annihilation Phenomenology
             The BKY and CD breaking schemes, together with the Primordial Mixing, allow
            to construct the non-anomalous and anomalous pieces of BHLS_2, a new (broken) HLS Lagrangian.
            This allows for  a simultaneous  study of the e^+ e^- annihilation into  6 final states (π^+π^-,
            π^0γ,  ηγ, π^+π^-π^0, K^+K^-, K_L K_S) within a unified
            framework where global fits can be worked out.
            
            The fits have been performed by taking into account the
            fact that most of the data samples carry normalization uncertainties, frequently dominant;
            these
            are appropriately considered in order to avoid biased results
            for physics quantities. The iterative method developed in <cit.>
            is thus intensively used in the global fits reported here; the convergence
            criterion chosen to stop the iteration at order n+1 was
            δχ^2= |χ^2_ tot(n+1)-χ^2_ tot(n)| < (0.3 ÷ 0.5). As
             χ^2_ tot≃ 1000 ÷ 1200 for ≃ 1150 ÷ 1240 data points,
            this convergence criterion is clearly constraining.
            
            The fitting procedure has allowed to detect (and discard) a very few data samples which
            hardly accommodate the global framework fed with almost all the existing
            data samples. Beside the already identified <cit.> two π^+π^-
            samples, one also found that, out of the 7 available π^+ π^- π^0 data samples,
            one  <cit.> hardly fits the global context, as also reported
            by others  <cit.>. Figures <ref>, <ref> and
            <ref> together with Table <ref> prove the fairly good quality of the
            global description for the π^+π^-, (π^0/η)γ and π^+π^-π^0
            channels. As for the 3-pion channel, one should note that BHLS_2 encompasses
            easily  the samples collected in the ϕ mass region and so, solves the difficulty
            encountered there  by the former BHLS.
            
            A thorough study of the available data samples has been performed concerning
            both the  K^+K^- and K_L K_S channels for which some tension is expected,
             as reported
            in <cit.>. The analysis in Section <ref> has shown that
            the SND, BaBar <cit.> and the updated <cit.>  CMD-2 data
            are in fairly good agreement with each other and accommodate easily the HLS context
            where all other kinds of data are encompassed. Tension has indeed been observed
            between the CMD-3 data <cit.> and the other samples.
            This tension has been traced back to the off-diagonal part of the CMD-3 error covariance
            matrices for both  the K^+K^- and K_L K_S spectra;
             this analysis is summarized by
             Table <ref>. In order to include the BaBar (and/or CMD-3 data) within the global framework,
            energy shifts relative to the CMD-2/SND energy scale  (chosen as energy scale
            reference) are
            mandatory and have been fitted at values consistent with the reported
            expectations <cit.>.
            To perform this exercise, when using the CMD-3 data, the error covariance
            matrices of these (only) have been amputated from their off-diagonal part.
            
            One can finally remark that, once the energy scales
            of the CMD-3 and BaBar <cit.> data  relative to CMD-2/SND are
            applied, one does not observe a visible tension among them any longer as illustrated
            by Figure  <ref>, where all data samples are nicely superimposing
            onto the same
            fit function.
            
            
            
             §.§ Meson Form Factors in the Spacelike Region
             The BHLS_2 form factors supplied with the parameter values derived
            from global fits to the  timelike region data only  have been
            shown to nicely predict the
            behavior of the π^± and K^± form factors   in the close spacelike region, as they
            can be expected from   NA7 and Fermilab data.
            Moreover, including these (model-independent) spacelike data within the BHLS_2 fit procedure
            does not exhibit  any gain in the description of the spacelike region as obvious from
            Table <ref> and from Figure <ref>.
            
            Our derived   π^± and K^± form factors,  which perform fairly well
            simultaneously in the spacelike and timelike regions,  strongly supports the validity of
            the BHLS_2 predictions in the s-region located in between. If the fairly good continuation
            of F_π^±(s) is actually performed over a limited gap in s to reach the spacelike region, one may
            wonder  about the so successful continuation of F_K^±(s) over an s-gap greater than
            1 GeV^2. This indicates that our modified Breit-Wigner formulas are performing (unexpectedly) well
            far beyond the ϕ energy region where their parameters are determined.
            
            The LQCD Collaborations HPQCD and ETM have provided F_π^±(s),
            the former as a para­metri­zation, the latter
            by corresponding spectra, covering the same space-like s-region as the NA7 experiment.
            Figure <ref> shows the superposition of the BHLS_2 prediction (without any free parameter)
            together with the
             parametrization provided by  HPQCD  and the two ETM spectra; the accord is very good for  HPQCD
             and simply perfect for both the ETM  spectra. There is, unfortunately, no reported LQCD data
             for the kaon form factors to which one could compare.
            
             The success just emphasized for the BHLS_2π^± and K^± form factors compared
             to spacelike experimental or LQCD data gives confidence in other physics information involving
             the chiral point region. Let us limit oneself, in this Subsection,
             to briefly comment on the neutral kaon
             electromagnetic form factor. As reminded in the main text, the slope for the F_K^0(s<0)
             form factor at s 0 is reported negative. This is qualitatively consistent with the BHLS_2 finding for
             F_K^0(s) shown in Figure <ref>, as it is found to carry a negative slope at negative s and
            a  positive slope at  positive s.
             However,  the effects of the neglected O(ϵ^2) corrections actually prevent to provide
             confidently accurate  slope predictions in the neighborhood of s=0 which is a
             delicate region for our O(ϵ) approximations.
            
            The pion form factor measured by the NA7 and Fermilab Collaborations are model independent and,
            actually, it cannot be measured in a model-independent way above |s| ≃ 0.25 GeV^2.
            To go beyond, one should accept to consider data where the extraction of  F_π^±(s)
            is model-dependent. Figure <ref> already indicates that the two  unfitted desy measurements fall (almost) exactly onto the fitting curve together with
             the model-independent  NA7 and Fermilab measurements. The
            JLAB data points are more or less closer to the BHLS_2 expectation depending on the method
            used to extract F_π^±(s). Nevertheless, our analysis  indicates that
            the BHLS_2 prediction could well be valid down to ≃ -1 GeV^2, possibly slightly beyond.
            
            In summary, our analysis clearly indicates
            that,  inherently, BHLS_2 only supplied with consistent e^+e^- annihilation data,
            leads precisely to the spacelike expectations supported by the precise NA7 data or the available LQCD input
            as well. Nevertheless, LQCD information on F_K^±(s) (and F_K^0(s)) would
            be welcome to  better check the open strangeness sector.
            
            
            
             §.§ Physical Quantities and Model Dependence Effects
             BHLS_2 has derived evaluations of several physics quantities (see Sections
            <ref>, <ref> and <ref>) and some   deserve special emphasis.
            As shown in Table <ref>, one reaches precise values for the charged pion and kaon
            charge radii (fm^2) :
            
                <r_π^±^2>=0.430 ± 0.002_ mod± 0.001_ fit   ,  
                <r_K^±^2>=0.268± 0.004_ mod± 0.001_ fit  ,
            
            our estimate for <r_π^±^2> is in fairly good agreement with the most
            recent evaluations <cit.>. Our evaluation for <r_K^±^2>
            is in accord with the reported values shown in Table <ref> and is more precise.
            The Reference Solution of BHLS_2 also leads to :
            
                f_K^±/f_π^±=1.247 ± 0.020_
                syst± 0.003_ fit ,
            
            where the quoted systematic error reflects different variants of the modeling.
            
            
            
            
            Table <ref> collects the values of the BHLS_2 model parameters obtained
            under several fit conditions and discussed in Section <ref>. An interesting
            numerical correlation is observed between the fundamental HLS parameter a  and the
            BKY breaking parameter z_V within the BHLS_2 framework. This can be approximated
            by  a z_V ≃ 2.6; good fit qualities are thus observed with  a  varying over
            a large interval, namely from 1.6 to 2.6. The rightmost two data columns in
            Table <ref>, which collect the fit values of a large sample of breaking parameters
            in  two extreme cases (a≃ 1.6 and a≃ 2.6), also show correlatedly important
            changes in some of the other breaking parameter values.  Table <ref> also shows that
            the fit qualities reached in all cases are almost identical.
            
            One can herefrom conclude  that
            the various model parameter sets efficiently conspire to provide almost identical
            cross-sections and form factors which are, actually, the real observables feeding BHLS_2.
            This has been substantiated by comparing the running of BHLS_2 fed with our consistent set
            of data samples under the various conditions
            underlying the fits reported in Table <ref>.
            
            A motivated piece of information to detect parametrization effects are
            the integrals a_μ(√(s)≤ 1.05 GeV) given in Table <ref>.
            Running the BS variant (excluding the τ data) provides
            a_μ(BS, √(s)≤ 1.05 GeV)= (571.84 ± 1.06) × 10^-10, while
            running the RS variant running on the same data (with the Primordial Mixing angles free) returns
            a_μ(RS, √(s)≤ 1.05 GeV)= (571.98 ± 1.20) × 10^-10. So the difference
            between the various parametrizations happens indeed to be marginal and quite comparable
            in magnitude to the uncertainties generated by using somewhat conflicting data samples.
            Nevertheless, external information
            constraining or relating  some of the model parameters could be useful in order to reduce
            the parameter freedom in fits and so the parameter correlations.
            
            
            
            
            As for the important muon HVP topic : Complementing the BHLS_2 contribution by 
            the whole hadronic contribution not covered by the (BHLS_2) model and reported in the left-hand side
            part of Table <ref>, one gets :
             
                10^10× a_μ^ HVP-LO = 686.65 ± 3.01
                 +[^+1.16_-0.75]_ syst +[^+3.10_-0.0]_ργ  .
            
             Its (statistical + model) uncertainty  (3.01) is  sharply dominated by those of 
            the whole non-HLS contribution : 2.76 × 10^-10– reported in Table <ref>–
            versus the BHLS_2 uncertainty (1.06 ÷ 1.20)× 10^-10 (see Table <ref>). 
            
            So, a significant
             improvement of the uncertainty for a_μ^ HVP-LO
             can only follow from an improved accuracy of the hadronic contributions of higher energies, especially
             those from the (1.05-2.00) GeV region.
              Therefore, presently, as global fit methods like BHLS_2
               constrain the central value for a_μ^ HVP-LO(HLS),
             they do alike for a_μ^ HVP-LO(HLS+non-HLS). The effect of the ρ^0-γ mixing vs LQCD
             is shown separately for clarity.
            
            Our estimate for the full muon anomalous moment is :
             
                10^10× a_μ= 11 659 175.96 ± 4.17_ th
                +[^+1.16_-0.75]_ syst+[^+3.10_-0.0]_ργ ,
            
            which is, at least, at  3.8 σ from the BNL measurement. The systematics are actually
            upper bounds to a possible bias affecting a_μ.
            
            
            
            
            § APPENDICES
            
            
            § THE FULL HLS NON-ANOMALOUS LAGRANGIAN IN THE NEW SCHEME
             For clarity,  the full non-anomalous HLS Lagrangian is written :
             L_ HLS = L_A+a  L_V+ L_τ+ L_p^4
            and  the various pieces will be given just below.
            
            
            
             §.§ The L_A+a  L_V Lagrangian Piece
            
            First, one displays the part of L_ VMD= L_A+ a  L_V relevant for the physics
            we address,  especially e^+ e^- annihilations. This is :
            
            [                                                                                             L_ VMD                                        =  ∂π^+ ·∂π^-
            + 1/2∂π^0 ·∂π^0
            + ∂ K^+ ·∂ K^- +  ∂ K^0 ·∂K^0;                                                                                                                                                            + ie [ 1 -a/2 (1+ Σ_V +Δ_V/3) ] A ·π^- π^+;                                                                                                                                                  + i e a/6z_A (1+Δ_A/2)(z_V-1 - Σ_V +Δ_V) A · K^0 K^0;                                                                                                                                  + ie [ 1 -a/6 z_A{ 2+z_V + 2  Σ_V +2  Δ_V
            }
            (1-Δ_A/2)]
             A · K^-  K^+;                                                                                                    +1/2 [
            m_ρ^0^2  (ρ^0_I)^2 + m_^2 ^2_I + m_ϕ^2 ϕ^2_I ]
            -e
             [f_ργ^I ρ^0_I+ f_γ^I _I + f_ϕγ^Iϕ_I] · A;                                                                                                                                                 + i a g/2 (1+ Σ_V){[ 1+ξ_3 ]ρ^0_I  + Δ_V _I
            }·π^- π^+;                                                                                                                                         +i a g/4 z_A [1-Δ_A/2]
            [ 1+ Σ_V +Δ_V +ξ_3 ] ρ^0_I ·  K^-  K^+;                                                                                                                 +i a g/4 z_A [1-Δ_A/2]
            [
            1+ Σ_V +Δ_V +2/3 (1-z_V)ξ_0 +1/3 (1+2z_V)ξ_8
            ]_I ·  K^-  K^+;                                                                                                                                                +i a g/4z_A[1+Δ_A/2] [1+ Σ_V -Δ_V+ξ_3]
            ρ^0_I · K^0 K^0;                                                                                                              -i a g/4z_A [1+Δ_A/2]
            [   1+ Σ_V -Δ_V+
            2/3 (1-z_V)  ξ_0+1/3 (1+2z_V)ξ_8
            ]  _I ·  K^0 K^0;                                                                                                                -i a g√(2)/4z_A  (1-Δ_A/2)
            [z_V  -1/3 (1-z_V) ξ_0 +1/3 (1+2z_V) ξ_8
            ]  ϕ_I ·  K^-  K^+;                                                                                                                  +i a g√(2)/4z_A  (1 + Δ_A/2)[
            z_V  -1/3 (1-z_V) ξ_0 +1/3 (1+2z_V)ξ_8
            ] ϕ_I ·
            K^0 K^0 ]
            The pseudoscalar fields occurring there are the renormalized ones defined exactly as
            in <cit.>.  For simplicity, we use  ideal vector fields.
            The  V-γ couplings are :
             {[                                                         f^I_γρ                              = m^2/g
            [1+ Σ_V  +Δ_V/3+  ξ_3]
             ,;                                                          f^I_γ = m^2/3g[1 + Σ_V +3Δ_V  +
            2/3(1-z_V)  ξ_0 + (1+2z_V)/3 ξ_8 ] ,;                                                         f^I_γϕ             =  -√(2) m^2/3g[z_V-ξ_0/3(1-z_V) + (1+2z_V)/3ξ_8
            ] ] .
            
            One has discarded the (irrelevant) photon mass term. The vector meson masses
            are given in Eq.
            (<ref>). Δ_V and Σ_V have been defined
            in  Section <ref>.
            
            In practical use, one should perform the change of field given in Eq. (<ref>)
            and collect the terms corresponding to each (neutral) vector meson coupling.
            
            
            
             §.§ The L_τ Lagrangian Piece
            
            The new  L_τ is given below (after performing an integration by parts
            to remove the terms with  ∂ W^±(=0) :
            [                                                                                                                   L_τ                         =  a g^2  f_π^2 (1+ Σ_V)ρ^+ ·ρ^-
            - i V_ud g_2/2√(2) W^+ · [1 - a/2 z_A   (1+ Σ_V ) ] K^0  K^-;                                                                                                                                       -1/2 a g g_2 f_π^2(1+ Σ_V)   [V_ud W^+ ·ρ^-
            +V_udW^- ·ρ^+ ] -i a g /2√(2)z_A (1+ Σ_V)
            ρ^+ ·  K^0  K^-;                                                                                                                                        + i V_udg_2/2 W^+ ·[ 1 -a/2 (1+ Σ_V ) ]
             [
            π^0  - ( ϵ_1 -Δ_A/2sinδ_P )  η
            -( ϵ_2 -Δ_A/2cosδ_P ) ]π^-;                                                                                                                       - i Δ_A V_udg_2/2 W^+ 
            
            
            
             [  sinδ_P η - cosδ_P ]  π^-  -  f_π g_2/4  W^+
            
             [ V_ud∂π^- +  V_us√(z_A) (1+Δ_A/4) ∂ K^-];                                                                                                                                                - i a g /2ρ^+ · [(1+ Σ_V) π^- π^0
            - ( ϵ_1 -Δ_A/2sinδ_P ) π^- η
            -( ϵ_2 -Δ_A/2cosδ_P ) π^- ] , ] where all fields are renormalized. δ_P is related with the PS mixing in 
            the one angle mixing scheme <cit.>; its definition is reminded in Appendix 
            <ref> below together with those for the ϵ_i. 
            
            The appearance of 
            W^± (π^0/η/η^')π^± couplings should be noted;
            one should also note, about the W^± (η/η^') π^± couplings,  that
            the contributions from L_A and L_V  are not proportional
            in contrast with the W^±π^0 π^± case.  Finally, it is also worth noting
            that the ρ^± (η/η^') π^± couplings in Equation (<ref>)
            are needed in the transitions amplitudes for η/η^'π^+ π^- γ 
            in order to yield the right behavior at the chiral point.
            
            
            
             §.§ The L_p^4 Lagrangian Piece L_z_3
            
            As for the L_z_3 Lagrangian, it writes in terms of ideal fields :
            [                                                                         L_z_3= e g z_3  s  A_μ [ (1+ξ_3)
            ρ^0 μ_I +1/3(1+ξ_8) ^μ_I
             -√(2)/3(1+ξ_8)
               ϕ^μ_I  ]                        +g_2 g z_3 s/2 [V_ud ρ^+_μ W^-_μ +
            	V_ud ρ^-_μ W^+ _μ ] ]
            It clearly allows for different Vγ and WV couplings, especially for the ρ meson.
            
            
            
             §.§ The Renormalized Couplings of Vector Mesons to K K
            
            Let us first remind the coupling of the fully renormalized vector mesons to
            π^+π^-; having defined these as g^V_π^±=ag g_π±/2,
            one has :
            [ g_π±^ρ=[1+Σ_V+ξ_3]  ,   g_π±^=-[ψ_+α(s)]  ,     g_π±^ϕ=[ψ_ϕ+β(s)] ]
            
            Using the Lagrangian in Section <ref> and Eq. (<ref>) one can derive the
            renormalized coupling of vector mesons to K K pairs. Defining:
            {[                          g_K±^ρ = (1-Δ_A/2 )
            {
            1 +Σ_V+ξ_3+ [ψ_+α(s)] +√(2) z_V [ψ_ϕ+β(s)]
            } ,; g_K±^ = (1-Δ_A/2 )
            {
            1 +Σ_V+2/3 (1-z_V)ξ_0 + 1+2z_V/3ξ_8
            - [ψ_+α(s)] +√(2) z_V [ψ_0+γ(s)]
            } ,;    g_K±^ϕ= (1-Δ_A/2 )
            {-√(2)[
            z_V-1/3 (1-z_V)ξ_0 +1+2z_V/3ξ_8
            ]
            + [ψ_ϕ+β(s)] + [ψ_0+γ(s)]
            } , ] .
            
            the couplings involved in form factors and cross-sections
            are given by G^V_K^+ K^-=ag /[4 z_A] g_K±^V. One also has :
            {[                           g_K0^ρ =  (1+Δ_A/2 )
            {
            1 +Σ_V+ξ_3- [ψ_+α(s)] -√(2) z_V [ψ_ϕ+β(s)]
            } ,; g_K0^ =-  (1+Δ_A/2 )
            {
            1 +Σ_V+2/3 (1-z_V)ξ_0 + 1+2z_V/3ξ_8
            + [ψ_+α(s)] +√(2) z_V [ψ_0+γ(s)]
            } ,;      g_K0^ϕ=  (1+Δ_A/2 )
            {√(2)[
            z_V-1/3 (1-z_V)ξ_0 +1+2z_V/3ξ_8
            ]
            + [ψ_ϕ+β(s)] - [ψ_0+γ(s)]
            } , ] .
            
            and the full couplings are given by
             G^V_K^0 K^0=ag /[4 z_A] g_K0^V. One  cancels out
             here and in the following the dependence upon Δ_V.
            One should also note that all these couplings depend on s, the square energy flowing through
            the relevant vector line.
            
            
             §.§ The V-γ Transition Loop Terms
            
            As seen in the main text, the V-γ transitions amplitudes write :
            
                F_Vγ^e(s)= f_Vγ - c_Vγ z_3 s - Π_Vγ(s)
            
            for each of V=ρ^0_R, _R, ϕ_R. In this Section, one displays
            the expression for the loop terms Π_Vγ(s) in terms of the basic loop
            functions denoted Π_π^e(s), Π_K±^e(s) and Π_K^0^e(s)
            and defined in Subsection <ref>. For the sake of conciseness, it is convenient
            to define :
            [                  g_π±^=
            [1-a/2[1+Σ_V  ] ]  ,; g_K±^=
            [1-a/6 z_A[2+z_V+2Σ_V -3/2Δ_A ] ]   ,                 g_K0^=
            -a/6 z_A(1-z_V+Σ_V) , ]
            neglecting the O(ϵ^2) terms. These couplings do not depend
            on the "angles" α(s), β(s) and γ(s).
            
            Using these definitions and those displayed in the Subsection just above, on finds :
            {[ Π_ρ=
            g_π±^g_π±^ρ/G_πΠ^e_π(s)
            +g_K±^g_K±^ρ/G_KΠ^e_K±(s)
            +g_K0^g_K0^ρ/G_KΠ^e_K0(s)  ,;    Π_=
            g_π±^g_π±^/G_πΠ^e_π(s)
            +g_K±^g_K±^/G_KΠ^e_K±(s)
            + g_K0^g_K0^/G_KΠ^e_K0(s)
             ,;   Π_ϕ=
            g_π±^g_π±^ϕ/G_πΠ^e_π(s)
            +g_K±^g_K±^ϕ/G_KΠ^e_K±(s)
            + g_K0^g_K0^ϕ/G_KΠ^e_K0(s) ].
            
            Expanding Equations (<ref>) in order to rather parametrize in terms of
            Π^e_K  sum and Π^e_K  diff
            (see Subsection <ref>) allows to exhibit the terms of order O(ϵ^2)
            in breakings which can be dropped out.
            
            
            
            § PSEUDOSCALAR FIELD RENORMALIZATION : A BRIEF REMINDER
            
            As already stated, one goes on using unchanged the breaking procedure in the L_A sector
            of the HLS model as defined in <cit.>. The  broken L_A Lagrangian
            displayed in Eq. (<ref>) implies a first step field redefinition in order to make the pseudoscalar
            kinetic energy term canonical. Additionally, as BHLS and BHLS_2 also include the 't Hooft determinant
            terms <cit.>, the  pseudoscalar kinetic energy receives an additional piece :
             L_ 'tHooft=1/2λ_μη_0^μη_0 + ⋯ ,
            
            where η_0 is the singlet PS field and
            λ is a parameter to be fixed. This term
            imposes a second step renormalization <cit.> of the (π^0, η,) sector.
            
            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 <cit.> :
             [                P^R_1=X_A^1/2 P X_A^1/2  , where   
            X_A= Diag(1 +Δ_A/2,1 -Δ_A/2,z_A) ]
            
            Some pseudoscalar (PS) fields happen to be fully renormalized after this first step; they are
             related to their  bare  partners by :
            {[                π^±_ bare                  = π^± ,;                K^±_ bare =1/√(z_A)(1-Δ_A/4)
            K^± ,;                K^0_ bare   =1/√(z_A)(1+Δ_A/4) K^0 ].
            
            These renormalized fields have already been used to derive the Lagrangians
            given in Appendix <ref>. As for the other PS fields, the first step
            renormalization leads to the R1 renormalized PS fields :
            {[                        π_0^bare =  π^R_1_0
            -Δ_A/2√(3)η_8^R_1
            -Δ_A/√(6)η_0^R_1 ,; η_0^bare =
             -Δ_A/√(6)π^R_1_0
            +√(2)/3z_A -1/z_Aη_8^R_1
            +1/32 z_A+1 /z_Aη_0^R_1 ,;  η_8^bare =
             -Δ_A/2√(3)π^R_1_0
            +1/3z_A  +2/z_Aη_8^R_1
            +√(2)/3z_A -1 /z_Aη_0^R_1 ].
            
            The second step renormalization required by the 't Hooft term
            gives the R renormalized fields <cit.> :
            {[                                    π^R_1_0=π^R_0 ,;   η^R_1_8=1+vcos^2β/1+vη^R_8
            -vsinβcosβ/1+vη^R_0 ,; η^R_1_0=-vsinβcosβ/1+vη^R_8
            + 1+vsin^2β/1+vη^R_0 , ].
            
            where :
            {[ cosβ=2z_A+1/√(3(2 z_A^2+1))  ,  sinβ=√(2)(z_A-1)/√(3(2 z_A^2+1));         v=√(1+λ(2 z_A^2+1)/3z_A^2) -1
            ≃λ/2(2 z_A^2+1)/3z_A^2   . ].
            
            The parameters affecting the PS fields submitted to fit are thus Δ_A, z_A
            and λ (or alternatively v).
            
            The π^0-η-η^' mixing is supplemented by defining
            the physically observable PS fields in terms of their R renormalized
            partners by <cit.> :
             {[                          π_R^0= π^0-ϵ_1  η-ϵ_2  η^' ,;   η_R^8= cosθ_P (η+ϵ_1  π^0)+
            sinθ_P (η^'+ϵ_2  π^0) ,; η_R^0 = -sinθ_P (η+ϵ_1  π^0)+
            cosθ_P (η^'+ϵ_2  π^0) , ].
            
            which exhibits three more parameters to be determined by fit.
            The singlet-octet  mixing angle θ_P occurring here is
             expected in the range (-10^∘÷ -15^∘) and is
            also algebraically related to the usual θ_8<cit.>.
            In order to express the forthcoming VVP and AVP coupling constants, the mixing angle δ_P
            defined by  <cit.> :
             {[ sinθ_P=1/√(3)
            (cosδ_P+√(2)sinδ_P) ,; cosθ_P=1/√(3)
            (√(2)cosδ_P-sinδ_P) , ].
            
            happens to be more appropriate than θ_P. On the other hand, it is also
            useful to define the combinations <cit.> :
            {[ x=1-3 z_A^2/2 z^2_A+1 v ,; x^'=1-3 z_A/2 z^2_A+1 v ,;      x^=1-3/2 z^2_A+1 v , ].
            
            where v is the nonet symmetry breaking parameter  defined in Equations(<ref>).
            x, x^' and x^  reflect the various ways by which
             nonet symmetry breaking occurs in the PS sector of broken HLS Lagrangians.
            
            
            
            
            § THE VVP LAGRANGIAN IN TERMS OF IDEAL VECTOR FIELDS
             The π^0, η and  pieces of the
            L_VVP  Lagrangian are given separately just below.
            If the PS fields occurring in this Appendix are renormalized by the procedure
            reminded in Appendix <ref>, it is simpler to express the various
            VVP Lagrangian pieces in terms of the ideal vector fields.
            The derivation of the couplings involving the renormalized fields is given in the main text.
            
            We split up the VVP Lagrangian into the pieces involving the renormalized
            π^0, π^±, η and  pieces, below.  As for the π^0, we have :
            
            
             L_VVP(π^0)=   C/2ϵ^μναβ{ [
              ( 1+2 ξ_0/3 + ξ_8/3 )
            ∂_μ_ν^I
            +√(2)/3 (ξ_0-ξ_8)
            ∂_μϕ^I_ν
             ] ×
                 [ ∂_αρ^+_βπ^- +∂_αρ^-_βπ^+
             + ( 1+ ξ_3 )∂_αρ^I_βπ^0 ] 
                  
            +  [g_π^0 ∂_μ^I_ν∂_α^I_β +
            g_ρπ^0 ∂_μρ^I_ν∂_αρ^I_β+
             g_ϕπ^0  ∂_μϕ^I_ν∂_αϕ^I_β +
             g_ρ^±π^0  ∂_μρ^+_ν∂_αρ^-_β] π^0  } , 
            
            where :
            {[                                     C=-N_c g^2 c_3/4 π^2 f_π ,; g^0_ρ^±π^0 =
            2  g^0_ρπ^0 =
            ϵ_2  cosδ_P
            -ϵ_1  sinδ_P
            -  Δ_A/2 ,;            g^0_π^0 =1/2 [ ϵ_2  cosδ_P
            -ϵ_1  sinδ_P
            -  Δ_A/2] ,;               g^0_ϕπ^0=-1/z_A√(2) [ ϵ_1  cosδ_P
            +ϵ_2  sinδ_P ] ].
            
            
            The VVη Lagrangian is given by :
            
            [                                                                                L_VVP(η)= C/2ϵ^μναβ{g^0_ϕη ∂_μ_ν^I
            ∂_αϕ^I_β
             +
             g^0_ρη ∂_μ_ν^I∂_αρ^I_β
            + g_ρ^±η ∂_μρ^+_ν∂_αρ^-_β .;                                                                                                         .
            +g_η ∂_μ^I_ν∂_α^I_β +
            g_ρη ∂_μρ^I_ν∂_αρ^I_β+
             g_ϕη ∂_μϕ^I_ν∂_αϕ^I_β}η , ]
            
            with :
            
            {[                   g^0_ϕη=1/3z_A [2 cosδ_P+ z_A √(2)sinδ_P ] (ξ_8-ξ_0)
             ,;                                       g^0_ρη =
             [
            Δ_A/2sinδ_P -ϵ_1  ]
             ,;                 g^0_ρ^±η = 1/3 [√(2) (1-x^')  cosδ_P- (1+2 x)sinδ_P ]
             ,; g^0_ρη= 1/6 [ √(2)  (1-x^')    cosδ_P
            -(1+2 x)sinδ_P  -6 ξ_3 sinδ_P ]
             ,;       g^0_η=1/6 [ √(2)  (1-x^')    cosδ_P
            -(1+2 x+ 4ξ_0+2ξ_8)sinδ_P ]
             ,;   g^0_ϕη=-1/6z_A [ √(2)  (2+x^ +2 ξ_0+4ξ_8)  cosδ_P
            - 2 (1-x^') sinδ_P ];                                                                          ].
            
            
            The VV Lagrangian is given by :
            
            [                                                                              L_VVP()= C/2ϵ^μναβ{g^0_ϕ ∂_μ_ν^I
            ∂_αϕ^I_β
             +
             g^0_ρ ∂_μ_ν^I∂_αρ^I_β
            + g_ρ^± ∂_μρ^+_ν∂_αρ^-_β .;                                                                                                     .
            +  g_ ∂_μ^I_ν∂_α^I_β +
            g_ρ ∂_μρ^I_ν∂_αρ^I_β+
             g_ϕ ∂_μϕ^I_ν∂_αϕ^I_β} , ]
            with :
            {[                   g^0_ϕ= 1/3z_A [2 sinδ_P- z_A √(2)cosδ_P ] (ξ_8-ξ_0)
             ,;                                       g^0_ρ =-
             [
            Δ_A/2cosδ_P +ϵ_2  ]
             ,;                g^0_ρ^± =  1/3 [√(2) (1-x^')  sinδ_P + (1+2 x)cosδ_P ]
             ,; g^0_ρ= 1/6 [ √(2)  (1-x^')    sinδ_P
            + (1+2 x)cosδ_P + 6 ξ_3 cosδ_P ]
             ,; g^0_= 1/6 [ √(2)  (1-x^')    sinδ_P
            + (1+2 x +4 ξ_0+ 2 ξ_8)  cosδ_P ]
             ,;   g^0_ϕ=-1/6z_A [ √(2)  (2+x^ +2 ξ_0+4 ξ_8)  sinδ_P
            + 2 (1-x^') cosδ_P ] ].
            
            Finally, it is worth extracting out from Eq. (<ref>) the part
            involving a charged ρ meson for illustrative purposes. This write :
            
             [                                                               L_VVP(ρ^±,π^0)= C/2ϵ^μναβ [
            g^±_ρ^±∂_μ_ν^I
            +g^±_ϕρ^±∂_μϕ^I_ν ]
            [ ∂_αρ^+_βπ^- +
            ∂_αρ^-_βπ^+] , ]
            with :
            
            
            [ g^±_ρ^±= ( 1+2 ξ_0/3 +ξ_8/3 )
              ,        g^±_ϕρ^±= √(2)/3 (ξ_0-ξ_8) ]
            
            
            § THE AAP AND APPP LAGRANGIANS IN THE HLS FRAMEWORK
             For the reader's convenience, we have found appropriate
            to remind the expressions for the AAP Lagrangian and for the part
            of the APPP Lagrangian relevant for the annihilation process
            e^+ e^- π^0π^+π^-. These have been derived in <cit.>.
            One first have :
             L_AAP=-3 α_ em/π f_π (1-c_4)ϵ^αβμν∂_α A_β∂_μ A_ν[
            g_π^0 γγπ^0/6 + g_ηγγη/2√(3)
             + g_η^'γγη^'/2√(3)] ,
            
            with :
            
            {[                                                           g_π^0 γγ=      1-5 Δ_A/6 +
             ϵ_1/√(3){5z_A-2/3z_Acosθ_P-√(2)5z_A+1/3z_Asinθ_P}              + ϵ_2/√(3){5z_A-2/3z_Asinθ_P+
            √(2)5z_A+1/3z_Acosθ_P} ,;                                                             g_ηγγ =                  cosθ_P/3{5z_A-2/3z_A(1+v)+v 1+2z_A/1+2z_A^2-Δ_A/2} -√(2)sinθ_P/3{5z_A+1/3z_A(1+v)+v 1-z_A/1+2z_A^2 -Δ_A/2}
            -ϵ_1/√(3) ,;                                                           g_η^'γγ =                 sinθ_P/3{5z_A-2/3z_A(1+v)+v 1+2z_A/1+2z_A^2
            -Δ_A/2}  +√(2)cosθ_P/3{5z_A+1/3z_A(1+v)+v 1-z_A/1+2z_A^2-Δ_A/2}-
            ϵ_2/√(3) , ].
            
            using the renormalized PS fields defined in Appendix <ref> and
            their specific parameters.
            
            The pion sector of the APPP Lagrangian is <cit.> :
            
            [ L_APPP=  -iE  g_π^0 π^+π^-ϵ^μναβ A_μ∂_νπ^0
             ∂_απ^- ∂_βπ^+  ,                                                              ]
            with
            
            E=-e /π^2 f_π^3[ 1- 3/4(c_1-c_2+c_4)
            ]
            
            and :
            
             g_π^0 π^+π^- =1/4[ 1-Δ_A/2 + cosθ_P/√(3){ϵ_1+√(2)ϵ_2}
             -sinθ_P/√(3){√(2)ϵ_1-ϵ_2}]
            One may prefer  reexpressing  this formula by :
            
            g_π^0 π^+π^- =1/4 [
            1-Δ_A/2 +ϵ_2 cosδ_P-ϵ_1 sinδ_P ] ,
            
            using the angle δ_P defined by Equations (<ref>).
            
            
            
            
            § THE VPPP LAGRANGIAN AND ITS RENORMALIZATION
             The VPPP anomalous Lagrangian in terms of  ideal
             vector fields can be written :
             {[                     L_VPPP=iDϵ^μναβ{ [ g_ρπ^0 ∂_νπ^0+g_ρη^0 ∂_νη+
            +g_ρη^'^0 ∂_νη^' ]  ρ^I_μ. ,; + [ g_π^0 ∂_νπ^0+g_η^0 ∂_νη+
            +g_η^'^0 ∂_νη^' ] ω_μ^I
            + g_ϕπ^0∂_νπ^0  ϕ_μ^I} ∂_απ^-  ∂_βπ^+ . ,;                                                     with   D=- 3 g (c_1-c_2-c_3)/4 π^2 f_π^3 , ].
            
            where one has limited oneself to display the P_0π^+π^- sector. The leading terms
            of these couplings are :
            
             {[ g^0_ρπ^0 =1/4 [
             Δ_A/2 +ϵ_1  sinδ_P
            -ϵ_2  cosδ_P ]
             ,;                 g^0_π^0 = -3/4 [ 1+ 2ξ_0 + ξ_8/3 ]
             ,;                       g^0_ϕπ^0=-√(2)/4 [ξ_0 -ξ_8  ] , ].
            {[ g^0_ρη=-1/12 [ (1-x^') √(2)cosδ_P-
            { (1+ 2 x)+ 3 ξ_3 } sinδ_P ]   ,;                                g^0_η= -3/4 [Δ_A /2 sinδ_P-ϵ_1
             ]
             ,;                                                        g^0_ϕη= 0  , ].
            {[ g^0_ρ=-1/12 [  (1-x^') √(2)sinδ_P
            +{ (1+ 2 x)+ 3 ξ_3 } cosδ_P ]
             ,;                                  g^0_= 3/4 [Δ_A /2 cosδ_P+ϵ_2
             ] ,;                                                           g^0_ϕ= 0 ].
            
            One may  note, once again  the clarity reached by using δ_P
            instead of θ_P.
            
            The couplings to the renormalized fields can be derived  from  those to the ideal fields
            by means of the R(s) matrix defined in the body of the text :
             
             g^R_P =  R(s) g^0_P_0
            the elements of the g^0_P_0 vectors being defined in  Equations (<ref>),
            Equations (<ref>) or Equations (<ref>) for resp. the P_0= π^0 or η,
            mesons.
            
            For the couplings of the ϕ meson, at leading  order in breaking, one gets :
              g^R_ϕπ^0(s) = -3/4{  [ψ_0 +γ(s)]
            -√(2)/3 [ξ_0-ξ_8] }
            One clearly sees that, at the ϕ mass,  one obtains
            a non-vanishing direct coupling of the ϕ to 3 pions, additionally
            s-dependent.
            
                              h-physrev
