Vacuum polarization contribution to muon $g-2$ as an inverse problem

We analyze the electromagnetic current correlator at an arbitrary photon invariant mass $q^2$ by exploiting its associated dispersion relation. The dispersion relation is turned into an inverse problem, via which the involved vacuum polarization function $\Pi(q^2)$ at low $q^2$ is solved with the perturbative input of $\Pi(q^2)$ at large $q^2$. It is found that the result for $\Pi(q^2)$, including its first derivative $\Pi^\prime(q^2=0)$, agrees with those from lattice QCD, and its imaginary part accommodates the $e^+e^-$ annihilation data. The corresponding hadronic vacuum polarization contribution $a^{\rm HVP}_\mu= (687^{+64}_{-56})\times 10^{-10}$ to the muon anomalous magnetic moment $g-2$, where the uncertainty arises from the variation of the perturbative input, also agrees with those obtained in other phenomenological and theoretical approaches. We point out that our formalism is equivalent to imposing the analyticity constraint to the phenomenological approach solely relying on experimental data, and provides a self-consistent framework for determining $a^{\rm HVP}_\mu$ in the Standard Model with higher precision.


I. INTRODUCTION
How to resolve the discrepancy between the theoretical prediction for the muon anomalous magnetic moment a µ = (g µ − 2)/2 in the Standard Model and its experimental data has been a long standing mission. The major uncertainty in the former arises from the vacuum polarization function Π(q 2 ) defined by an electromagnetic current correlator at a photon invariant mass q 2 , to which various phenomenological and theoretical approaches have been attempted. For instance, the measured cross section for e + e − annihilation into hadrons has been employed to determine the hadronic vacuum polarization (HVP) contribution in a dispersive approach, giving a HVP µ = (693.9 ± 4.0) × 10 −10 [1] (see also a HVP µ = (692.78 ± 2.42) × 10 −10 in [2]). This value, consistent with earlier similar observations [2][3][4][5], corresponds to a 3.3σ deviation between the Standard Model prediction for a µ and the data [6], a exp µ − a SM µ = (26.1 ± 7.9) × 10 −10 . The above phenomenological determinations of a HVP µ , solely relying on experimental data, suffers a difficulty: the discrepancy among individual datasets, in particular between the BABAR and KLOE data in the dominant π + π − channel, leads to additional systematic uncertainty [1]. Therefore, theoretical estimates of the HVP contribution to the muon g − 2 are indispensable, which have been performed mainly in lattice QCD (LQCD) (see [7] for a recent review). Results, such as a HVP µ = (654 ± 32 +21 −23 ) × 10 −10 in [8], are comparable to those from the phenomenological approach. It has been known that the finite volume in LQCD makes it unlikely to compute the vacuum polarization at low momenta with high statistics, for which a parametrization is always required to extrapolate lattice data.
In this paper we will calculate the vacuum polarization function in a novel method proposed recently [9], where a nonperturbative observable is extracted from its associated dispersion relation. Taking the D meson mixing parameters as an example [9], we separated their dispersion relation for D mesons of an arbitrary mass into a low mass piece and a high mass piece, with the former being regarded as an unknown, and the latter being input from reliable perturbation theory. The evaluation of the nonperturbative observable is then turned into an inverse problem: the observable at low mass is solved as a "source distribution", which produces the "potential" at high mass. The resultant Fredholm integral equation allows the existence of multiple solutions as a generic feature. However, it has been demonstrated that nontrivial solutions for the D meson mixing parameters can be identified by specifying the physical charm quark mass, which match the data well. This work implies that nonperturbative properties can be extracted from asymptotic QCD by solving an inverse problem.
Here we will solve for the vacuum polarization function Π(q 2 ) via an inverse problem, and derive the HVP contribution a HVP µ to the muon g−2. The electromagnetic current correlator is decomposed into three pieces according to the quark composition of the ρ, ω, and φ mesons. A dispersion relation is considered for each resonance, and converted into a Fredholm integral equation, which involves the unknown constant Π(q 2 = 0) and the imaginary part ImΠ(q 2 ) corresponding to the e + e − → (ρ, ω, φ) → hadron spectra of nonperturbative origin. We solve the Fredholm equation with the perturbative input of the leading order correlator at large q 2 , and select the solution, which best fits the e + e − annihilation data for the resonance spectra. The determined Π(0), together with the resonance spectra at low q 2 and the perturbative input at high q 2 , then yields Π(q 2 ) from the dispersion relation. It will be shown that our predictions for Π(q 2 ), including its first derivative Π (q 2 = 0), and for a HVP µ agree with those obtained in the literature.
We point out that the above formalism is equivalent to imposing the analyticity constraint to the phenomenological approach solely relying on experimental data: when fitting the data, we search only parameters involved in ImΠ(q 2 ) that satisfies the Feldholm equation, instead of tuning them arbitrarily. To impose the analyticity constraint, one may, for instance, check whether a dataset obeys the Fredholm equation, namely, whether its dispersive integral reproduces the perturbative Π(q 2 ) at large q 2 , so that inconsistent datasets, such as the BABAR and KLOE data mentioned above, can be discriminated, and the precision in the individual datasets can be fully exploited. That is, we provide a self-consistent theoretical framework for deriving a HVP µ in the Standard Model with a potential to reach higher precision.
The rest of the paper is organized as follows. In Sec. II we present our formalism for extracting the nonperturbative vacuum polarization function Π(q 2 ) at low q 2 , and solve the corresponding Fredholm equation. The similar procedure is extended to compute the slope Π (0), that gives the leading contribution in the representation of Π(q 2 ) − Π(0) in terms of Padé approximations [10][11][12], and serves as a key ingredient in the "hybrid" approach proposed in [13]. We evaluate the HVP contribution to the muon anomalous magnetic moment numerically in Sec. III, and compare our prediction a HVP µ = (687 +64 −56 ) × 10 −10 , where the uncertainty comes from the variation of the perturbative input, with those from other phenomenological and LQCD approaches. Section IV is the conclusion.

II. THE FORMALISM
Start with the correlator The leading order expression for the HVP contribution to the muon anomalous magnetic moment is written, in terms of the vacuum polarization function Π EM (q 2 ), as [14,15] a HVP µ = 4α 2 with the electromagnetic fine structure constant α EM and the muon mass m µ . The first term can be set to Π EM (0) = 0 [16] in the on-shell scheme for the QED renormalization, but is kept for generality, because it also receives nonperturbative QCD contribution. The behavior of Π EM (−s) in the region with a large invariant mass squared s has been known in perturbation theory. We will derive Π EM (−s) in the low s region, where the nonperturbative contributions from the ρ, ω, and φ resonances dominate. The vacuum polarization function obeys the dispersion relation λ being a threshold. The function Π EM (s) for large s can be expressed as m f being a light quark mass. The real parts of the functions Π(s, m f ) at leading order are read off [17] up to an overall normalization, with s > 0 in the on-shell, MS and MS schemes for the QED renormalization, respectively. The imaginary part is given by [17] ImΠ(s, m f ) = 0, s < 4m 2 It is seen that the real parts Π EM (−s) in the above schemes differ by the s-independent terms, which can be always absorbed into the redefinition of the unknown constant Π EM (0) in Eq. (3). It is also clear that our result for a HVP µ will not depend on the choice of a specific renormalization scheme, because the scheme dependence cancels between the two terms in Eq. (2). Hence, we will stick to the on-shell scheme, and omit the subscript OS in the formulation below. We decompose Eq. (3) into three separate dispersion relations labelled by r = ρ, ω, φ, and rewrite them as where the thresholds are set to λ ρ = 4m 2 π + , λ ω = (2m π + + m π 0 ) 2 , and λ φ = 4m 2 K + , with the pion (kaon) mass m π (m K ). The separation scale Λ r will be determined later, which is expected to be large enough to justify the perturbative calculation of the imaginary part ImΠ r (s) in Eq. (8). Equation (7) is then treated as an inverse problem, i.e., a Fredholm integral equation, where Ω r (s) defined by Eq. (8) for s > Λ r is an input, and ImΠ r (s) in the range s < Λ r is solved with the boundary condition ImΠ r (λ r ) = 0 and the continuity of ImΠ r (s) at s = Λ r . That is, the "source distribution" ImΠ r (s) will be inferred from the "potential" Ω r (s) observed outside the distribution. Equation (7) can be regarded as a realization of the global quark-hadron duality postulated in QCD sum rules [18].
Both the real and imaginary parts of the input functions Π r (s) in Ω r (s) are related to Π(s, m f ) via with the charge factors and Ω ρ (s) in Eq. (8) It is easy to find that any two adjacent rows of the matrix A approach to each other as the grid becomes infinitely fine. Namely, A tends to be singular, and has no inverse. We stress that this singularity, implying no unique solution, should be appreciated actually. If A is not singular, the solution to Eq. (7) will be unique, which must be the perturbative results in Eqs. (5) and (6). It is the existence of multiple solutions that allows possibility to account for the nonperturbative ImΠ r (s) in the resonance region. After solving for Π r (0) together with ImΠ r (s) in the whole range of s, we derive Π r (−s) from the three dispersion relations, and Π EM (−s) from their sum to be inserted into Eq. (2).
Knowing the difficulty to solve an inverse problem and the qualitative behavior of a resonance spectrum, we propose the parametrizations according to [19,20], where d r = m r Γ r is the product of the meson mass m r and the width Γ r . The first factors in the above expressions guarantee the vanishing of the resonance spectra at the thresholds. We have adopted the same threshold for the K + K − , K S K L and π + π − π 0 final states of φ decays for simplicity. The parameter κ characterizes the ρ-ω mixing effect, and b r 0 (c r 0 ) describes the strength of the resonant (nonresonant) contribution. The parameters z 1 and z 2 , motivated by the Gounaris-Sakurai mode [21], lead to the effective width and mass of a ρ meson. This can be understood by completing the square of the denominator of the resonance term, with the quartic term being left aside first. The z 1 term then shifts the ρ meson mass and width into m 2 valid for |s| m 2 ρ will be assumed in the numerical analysis. We have confirmed that the quartic term is much less important than the quadratic term in the denominator even for s ∼ m 2 ρ and z 1 and z 2 determined later, so the approximation indeed holds well.
We have examined that the variations of the meson masses m r and widths Γ r and the ρ-ω mixing parameter κ change our results at 0.1% level, so m r and Γ r are set to their values in [6], and the mixing parameter is set to κ = 2.16 × 10 −3 [1]. The free parameters z 1 , z 2 , b r 0 , c r 0 , Λ r and ImΠ r (0) are then tuned to best fit the input Ω r (s) under the continuity requirement from ImΠ r (s = Λ r ). The separation scale Λ r introduces an end-point singularity into Ω r (s) in Eq. (8) as s → Λ r . To reduce the effect caused by this artificial singularity, we consider Ω r (s) from the range 15 GeV 2 < s < 250 GeV 2 , in which 200 points s i are selected. We then search for the set of parameters, that minimizes the residual sum of square (RSS) Such a set of parameters corresponds to a solution of the Fredholm equation in Eq. (7) in terms of the parametrizations in Eq. (10).

III. NUMERICAL ANALYSIS
The scanning over all the free parameters reveals the minimum distributions of the RSS defined in Eq. (11). Typical distributions on the Λ r -Π r (0) plane for other given parameters are displayed in Fig. 2. The minima along the curve, having RSS about 10 −12 -10 −13 relative to 10 −8 from outside the curve, hint the existence of multiple solutions. A value of Λ r represents the scale, at which the nonperturbative resonance solution starts to deviate from the perturbative input. This explains the dependence on Λ r of a solution. It is observed that the solutions for Π r (0), including its sign and magnitude, fall in the same ballpark as LQCD results [8]. In order to select a solution from the multiple solutions, a guidance from experimental data is needed. For the ρ resonance spectrum, we are guided by the SND data in [22], which are consistent with those from all other collaborations as indicated by Fig. 5 in [1]. We do not consider the BABAR and KLOE data due to the discrepancy between them, as mentioned in the Introduction. It means that we are making a conservative prediction for the HVP contribution to the muon anomalous magnetic moment. We are guided by the data for the process e + e − → π + π − π 0 through the ω resonance in [23]. For the φ resonance, the SND data [24] are also adopted, which include the e + e − → K + K − , K S K L , and π + π − π 0 channels. We search for the parameters along the minimum distributions to select the solutions, for which the above e + e −annihilation data are best accommodated, and find  The SND data [22][23][24] are also exhibited for comparison. The data for the three modes e + e − → φ → π + π − π 0 , KSKL and K + K − have been combined with their uncertainties being added in quadrature. for the ρ, ω and φ resonances. The values of Λ r from the best fit are large enough for justifying the perturbative evaluation of the input Ω r (s). Note that the above parameters follow the correlation demanded by the perturbative input via the Fredholm equation, and are not completely free. This correlation, originating from the analyticity of the vacuum polarization, renders our approach distinct from the phenomenological one [1][2][3][4][5], in which the free parameters are solely determined by data fitting. We emphasize that a sensible resonance spectrum should be a solution of the Fredholm equation, i.e., respect the analyticity of the vacuum polarization. Therefore, one may check whether a dataset obeys the Fredholm equation, i.e., whether its dispersive integral reproduces the perturbative vacuum polarization function at large s, before it is employed in the phenomenological approach. This check will help discriminating inconsistent datasets, such as the BABAR and KLOE data mentioned before, and enhancing the precision of the obtained hadronic contribution to the muon g − 2.
The predicted cross sections corresponding to the sets of parameters in Eq. (12) are shown in Fig. 3, which agree with the measured ω and φ resonance spectra well, but deviate from the ρ spectrum slightly. The agreement is nontrivial, viewing the correlation imposed by the analyticity constraint on the parameters. A parametrization more sophisticated than Eq. (10), e.g., the one proposed in [1] below the threshold of the inelastic scattering may improve the agreement in the ρ channel. However, we will not attempt this, but investigate whether the theoretical uncertainty in the present analysis can explain the deviation. The inclusion of higher order QCD corrections to the perturbative input is expected to cause about α s /π ∼ 20% variation at the scale of Λ r around few GeV 2 . We thus increase and decrease the perturbative inputs by 20%, repeat the above procedure to fix the parameters, and use the selected best solutions to estimate the bounds of our predictions, which are represented by the bands in Fig. 3. The bands associated with the ω and φ spectra are too thin to be seen. It is found that most data for the ρ spectrum are covered, except the tail part at low s, which contributes little to a HVP µ anyway. It implies that the estimate of the theoretical uncertainty through the variation of the perturbative input is relevant. Certainly, different choices of the parametrizations for the resonance spectra may also cause theoretical uncertainty. Because our predictions have matched the data satisfactorily, we do not take into account this source of uncertainty here.
Once the imaginary part ImΠ r (s) at low s is derived, its behavior in the whole s range is known (with the perturbative input at high s), and the real part Π r (−s) can be calculated from Eq. (3). The results of the vacuum polarization functions in both the space-like s < 0 and time-like s > 0 regions are presented in Fig. 4. The oscillations of the curves ought to appear, when the photon invariant mass crosses physical resonance masses. The predicted vacuum polarization function from the u and d quark currents, i.e., the ρ and ω meson contributions is exhibited in Fig. 5. In order to compare our result with Π(Q 2 ) ud in LQCD [8], where a photon invariant mass is defined in the Euclidean momentum space, we have converted Eq.
It is obvious that our prediction for Π(Q 2 ) ud agrees with the LQCD one corresponding to the pion mass m π = 185 MeV within the theoretical uncertainty. The LQCD results show the tendency of decreasing with the pion mass, so a better agreement is expected, if a further lower pion mass could be attained.
With the vacuum polarization functions Π r (s) being ready in the whole s range and the relation Π EM (s) = to the muon anomalous magnetic moment, where the uncertainty comes from the variation of the pertubative inputs by 20%, and mainly from the ρ channel. The decomposition of the central value into the three pieces of resonance contributions gives a HVP,ρ µ = 585 × 10 −10 , a HVP,ω µ = 53 × 10 −10 , and a HVP,φ µ = 49 × 10 −10 . All the above results, consistent with those in the literature [8], imply the success of our formalism: nonperturbative properties can be extracted from asymptotic QCD by solving an inverse problem.
A hybrid method has been proposed in [13], which combines the data fitting and the LQCD input for the first derivative of the vacuum polarization function Π EM (0). The final expression for the light-quark HVP contribution to the muon anomalous magnetic moment is written as where the first error largely stems from the data of the e + e − annihilation cross section. The first derivative in the second term is given by the sum Π EM (0) = r=ρ,ω,φ Π r (0) with each piece where the determined parameters in Eq. (12) are taken for the first integral, and the perturbative input is inserted into the second integral. Equation (15) then yields the first derivatives at the origin Π ρ (0) = 0.0870, Π ω (0) = 0.0076, Π φ (0) = 0.0067, which are scheme-independent, though the on-shell scheme has been adopted here. Substituting Eq. (16) into Eq. (14), we have This value, turning out to be very close to that in [1], further supports our formalism for evaluating the vacuum polarization. The accuracy of a calculation in the hybrid approach can be improved by including higher derivatives of the vacuum polarization function [13], which are not yet available in LQCD, but can be derived using our formalism.
At last, we present an alternative expression for the vacuum polarization function, which may be considered for a hybrid approach. Starting with Eq. (3) and following the idea of [13,25,26], we write where the first and third terms can be computed in perturbation theory for a large enough scale Λ, and the second term, receiving the low mass contribution, can take the data input.

IV. CONCLUSION
In this paper we have extended a new formalism for extracting nonperturbative observables to the study of the HVP contribution a HVP µ to the muon anomalous magnetic moment g−2. The dispersion relation for the vacuum polarization function Π(q 2 ) was turned into an inverse problem, via which Π(q 2 ) at low q 2 was solved with the perturbative input of Π(q 2 ) at high q 2 . Though multiple solutions exist, the best ones can be selected, which accommodate the data of the e + e − annihilation cross section. Because the involved parameters are correlated under the analyticity requirement of the vacuum polarization, and not completely free, the satisfactory agreement of our solutions with the data is nontrivial. Our work suggests that this analyticity constraint can be included into the conventional phenomenological approach solely relying on data fitting to form a more self-consistent framework for the determination of the HVP contribution with higher precision. It has been found that our prediction for Π(q 2 ), including its first derivative Π (0), is close to those from LQCD, and contributes a HVP µ = (687 +64 −56 ) × 10 −10 to the muon g − 2 in consistency with the observations obtained in the other phenomenological, LQCD and hybrid approaches.
We state again that the purpose of this work is not to fit the e + e − annihilation data exactly, but to demonstrate how our formalism is implemented, and that reasonable results can be produced even with a simple setup like the leading order perturbative input and the naive parametrizations in Eq. (10). It has been shown that the slight deviation of our prediction for the ρ resonance spectrum from the data could be resolved by considering subleading contributions to the perturbative input. This subject will be investigated systematically in a forthcoming publication, and the corresponding theoretical uncertainty is expected to be reduced. Other sources of uncertainties need to be examined, such as the one from different parametrizations for the resonance spectra. The success achieved in this paper also stimulates further applications of our formalism to the hadronic contributions to the muon g − 2 from heavy quarks and from the light-by-light scattering [16,[27][28][29], for which a lack of experimental information persists, and a theoretical estimation is crucial.