Four-Flavour Leading-Order Hadronic Contribution To The Muon Anomalous Magnetic Moment

We present a four-flavour lattice calculation of the leading-order hadronic vacuum polarisation contribution to the anomalous magnetic moment of the muon, $a_\mathrm{\mu}^{\rm hvp}$, arising from quark-connected Feynman graphs. It is based on ensembles featuring $N_f=2+1+1$ dynamical twisted mass fermions generated by the European Twisted Mass Collaboration (ETMC). Several light quark masses are used in order to yield a controlled extrapolation to the physical pion mass. We employ three lattice spacings to examine lattice artefacts and several different volumes to check for finite-size effects. Incorporating the complete first two generations of quarks allows for a direct comparison with phenomenological determinations of $a_\mathrm{\mu}^{\rm hvp}$. Our final result including an estimate of the systematic uncertainty $$a_{\mathrm{\mu}}^{\rm hvp} = 6.74(21)(18) \cdot 10^{-8}$$ shows a good overall agreement with these computations.


Introduction
The detection of a scalar boson at the LHC leaves us, in case that it turns out to be the standard model Higgs boson, with a somewhat puzzling situation. On the one hand, the standard model seems then to be complete. On the other hand, e.g. the lack of our understanding of dark matter, the amount of CP violation, and the generation of baryon asymmetry in the universe strongly suggest that new physics must exist. It is therefore a most important task, to find indeed signs of this new physics beyond the standard model by looking at promising physical observables to detect its manifestations.
The anomalous magnetic moment of the muon, a µ , can be considered as such an important benchmark quantity to test the standard model of particle interactions, or, alternatively to detect new physics beyond the standard model. It can be measured experimentally very precisely [1,2] and can likewise be computed precisely in the standard model, see e.g. the review [3]. A comparison between the experimental result for a µ and a standard model calculation reveals a discrepancy of about three σ which is persistent over many years now and has been confirmed by computations of a number of groups. The interesting question is then, whether this discrepancy originates from some undetected effect in the experimental or theoretical determination of a µ , or, somewhat more excitingly, whether it points to physics beyond the standard model.
A key ingredient in the calculation of a µ is the leading order hadronic vacuum polarisation contribution, a hvp µ , which presently is the largest source of uncertainty of the theoretical computation of a µ . Since a hvp µ is of intrinsically non-perturbative nature, evidently a lattice QCD computation for this quantity is highly desirable. In fact, in a recent study, using only two flavours of mass degenerate quarks, it could be demonstrated that by employing an improved method to compute a hvp µ on the lattice, an unprecedented precision for a lattice calculation of a hvp µ could be reached [4,5]. In this paper, we want to report on a calculation of a hvp µ with four flavours of quarks, using besides mass degenerate up and down quarks also the strange and the charm quarks. Having such a setup is important for two reasons. First, since the contribution of the bottom quark is so small that neither experimental nor theoretical determinations are currently sensitive to it, the standard model calculation is sensitive to four flavours and disentangling the flavour contributions, for e.g. only the light flavours, is afflicted with ambiguities. Using four flavours is therefore cleaner and allows for a direct and unambiguous comparison to phenomenological determinations.
The second reason is that the charm quark contribution computed in perturbation theory [6] is a hvp µ,c = 1.44(1)×10 −9 . Hence, the charm quark contribution has approximately the same size as the light-by-light contribution [7] and the electroweak contributions [3] which are larger than the current experimental and theoretical uncertainties. Including the charm quark contribution in a lattice calculation is therefore important if values of a hvp µ are to be evaluated with a precision aiming to match the experimental errors. For our computation of a hvp µ we will follow closely the strategy of Refs. [4,5] applying an improved lattice definition of a hvp µ . As we will see below, also for the four flavour calculation discussed here, this method allows for a controlled extrapolation to the physical point leading to an accurate determination of a hvp µ . For our calculations we employ configurations generated for four flavours by the European Twisted Mass Collaboration (ETMC) [8,9]. These sets of configurations are obtained at different values of the lattice spacing and several volumes, thus allowing us to estimate discretisation and finite size effects as systematic uncertainties in our lattice calculation. In addition, at each value of the lattice spacing configurations exist at several values of the pion mass, ranging from 230 MeV m π 450 MeV thus allowing an extrapolation to the physical point.
In this paper, we will give a detailed discussion about the extraction of a hvp µ , providing our fitting strategy for the vacuum polarisation functions, and also address Padé approximants to analyse the vacuum polarisation function as suggested in [10].

Basic definitions
The key quantity for the determination of the leading-order hadronic contribution to the muon g − 2 is the hadronic vacuum polarisation tensor Π em µν (Q) with Euclidean momentum Q. It can be obtained from the correlator of two electromagnetic vector currents in the following way The form of J em µ given above anticipates that our simulations involve four dynamical quark flavours, up(u), down(d), strange(s), and charm(c). Euclidean symmetry and the Ward-Takahashi identities require the vacuum polarisation tensor to be transverse, i.e.
is the hadronic vacuum polarisation function, for which the label "em" will be left out in the following to ease notation. Its renormalised variant appears in the expression for the leading hadronic contribution to the anomalous magnetic moment of the muon in Euclidean space-time [11,12] a hvp 5) in which α denotes the fine-structure constant and m µ the muon mass. Since the weight function w Q 2 /m 2 µ is known from leading-order QED perturbation theory to have the form w(r) = 64 the main task of the lattice calculation is the determination of Π R (Q 2 ) from the vector current correlator.

Setup
As mentioned before, the calculations are performed employing gauge field configurations generated by the ETMC with N f = 2+1+1 dynamical quark flavours. In the gauge sector the Iwasaki action [13] is used. The up and down quarks reside in a mass-degenerate flavour doublet χ l = u d . They are described by the standard twisted mass lattice action [14] S is the Wilson Dirac operator with covariant forward and backward derivatives ∇ µ and ∇ * µ , respectively. The value of m 0 has been tuned to its critical value in order to ensure automatic O(a) improvement of all physical observables. The bare twisted quark mass is denoted by µ and τ 3 is the third Pauli matrix acting in flavour space. For the heavy quarks χ h = c s the twisted mass action for a non-degenerate fermion doublet [15] has been implemented with the same value for m 0 . τ 1 is the first Pauli matrix, µ δ denotes the mass splitting between the charm and the strange quarks, and µ σ the bare twisted mass with a twist in flavour space orthogonal to the one in the light sector. Details of the analysed ensembles are given in Table 1. They involve different volumes, three different lattice spacings, and up to five pion masses per lattice spacing such that finite-size effects, lattice artefacts, and the chiral extrapolation can be studied. Parameters of the N f = 2 + 1 + 1 flavour gauge field configurations that have been analysed in this work. β denotes the gauge coupling, a the lattice spacing, L a 3 × T a the spacetime volume, and m P S is the, unphysical, value of the pion mass. The values for m P S have been determined in [8]. L is the spatial extent of the lattices. The approximate lattice spacings given here are taken from a first analysis of the used gauge field configurations [16].
As mentioned before, the calculation of Π(Q 2 ) requires the evaluation of the vectorvector correlation function in Eq. (2.2). On the lattice, one can use the local vector current J L µ = χ(x)Q el γ µ χ(x) in complete analogy to the continuum with the charge matrix Q el = diag(2/3, −1/3), where χ stands either for the light or the heavy doublet. However, due to the lattice regularisation this current has to be renormalised multiplicatively and does not satisfy a Ward-Takahashi identity. To avoid these issues we instead use the conserved Noether current also known as point-split vector current at sink and source This current is defined for each flavour separately using the decomposition of the electromagnetic vector current of Eq. (2.1) into the quark currents J u , J d , J c , and J s Since the Noether current is only conserved for actions diagonal in flavour space, we employ an action for the heavy valence quarks different from the sea quark action given in Eq. (3.3), namely the so-called Osterwalder-Seiler (OS) action [17,18] How to achieve O(a) improvement in this situation is discussed in Ref. [18]. The bare twisted mass parameters for the valence strange and the charm quarks, µ s and µ c , are tuned in such a way that the physical values for 2m 2 K − m 2 PS and the D-meson mass, respectively, are reproduced. Here, m K denotes the kaon mass. This leads to the following values .
With our definition of the vacuum polarisation tensor involving the contact term we can now rely on the Ward identity being valid at non-zero lattice spacings a. Hence, the vacuum polarisation tensor on the lattice satisfies in momentum space The reason for the O(a) lattice artefacts is the lack of full Euclidean symmetry due to the lattice discretisation. Note that nevertheless, as discussed below, the vacuum polarisation function can be made O(a)-improved.

Vacuum polarisation
Plugging Eq. (2.1) into Eq. (2.2) and noting that there are no mixed-flavour contributions for the quark-connected Feynman diagrams, we can obtain the total vacuum polarisation tensor by adding the vacuum polarisation tensors of the single-flavour contributions including the appropriate charge factors Notice that we have defined the single-flavour contributions to Π tot µν (Q) excluding the factors of the squared charges. In the following we will generically denote them by Π µν (Q) and the respective vacuum polarisation function by Π(Q 2 ). The factor of Π ud µν (Q) is the sum of the squared charges of up and down quarks. Similarly, the total a hvp µ is procured as a sum of the single-flavour contributions including the charge factors.
For constructing the vacuum polarisation function Π(Q 2 ) from the vacuum polarisation tensor Π µν (Q) in Eq. (3.9), we define for the single-flavour contributions as the real part of the projection onto the tensor is defined by all momenta related by a Euclidean space-time symmetry transformation, although these are not all exact symmetries at non-zero lattice spacing and finite volume. It has been checked that this does not lead to observable effects in Π(Q 2 ) in the low-momentum region, Q 2 ≤ 2 GeV 2 , which is by far the most important for the determination of a hvp µ . In particular, the definition in Eq. (3.11) involves summing over Q and −Q which is required for achieving O(a) improvement.
Following a similar line of arguments as given in [19,20], it can be demonstrated that also for the case of the vacuum polarisation function defined as in Eq. (3.11) short-distance singularities do not spoil automatic O(a) improvement at maximal twist [21]. Therefore, continuum extrapolations are performed with an a 2 -term in Sec. 4.
Due to the renormalisation of the vacuum polarisation function given in Eq. (2.4), the calculation of a hvp µ for each flavour involves both an interpolation of Π(Q 2 ) in between discrete lattice momenta as well as an extrapolation to zero momentum. Hence, the choice of fit function is important, especially in the low-momentum region where the weight function in Eq. (2.6) is peaked, and its influence on the calculation of a hvp µ has to be investigated.
Our standard fit of the low-momentum dependence of the vacuum polarisation function involves two different terms. The first is inspired by vector meson dominance with M vector meson mass poles whereas the second term parametrises remaining deviations in the low-momentum region which extends up to a matching momentum Q 2 match . Here, m i denotes the mass of the vector meson states and f i their decay constants. They are determined before fitting the vacuum polarisation from the same vector correlation functions partially Fourier transformed in the spatial directions on the same bootstrap samples as will be described in the next section. Thus, the a j are the only parameters fitted here. For the matching of low and high momentum fit functions we have used Q 2 match = 2 GeV 2 . The high-momentum fraction of Π(Q 2 ) for Q 2 > Q 2 match is interpolated, inspired by perturbation theory, using a polynomial in Q 2 and a polynomial multiplied by a logarithmic term Here, the fit parameters are denoted b k and c l . The total vacuum polarisation function for each flavour is then obtained from with Θ(Q 2 ) denoting the Heaviside step function.  Table 1 for details).
For the light and the strange quark contributions to Π tot µν (Q) our standard fit is characterised by M = 1, N = 2, B = 4, and C = 1. An example fit for the vacuum polarisation for one of our lightest pion masses (B25.32t in Table 1) is given in Fig. 1. Since the curvature of the vacuum polarisation function for the charm quark contribution is much smaller, more parameters are needed in the high-momentum domain in order to ensure smooth contact between the fit functions in the low and high momentum regions. Thus we employ for the charm quark contribution M = 1, N = 2, B = 3, and C = 5 in our standard fit. The effect of varying the number of terms has been checked and will be reported when discussing the systematic uncertainties of our calculation.
To minimise finite-size effects we have excluded the points with the lowest lattice momentum in the vacuum polarisation fits which due to the big uncertainty of these points has only a very mild influence on the fits. This has been suggested by a treelevel study showing that except for the lowest momentum point all other data points for different volumes fall on top of each other.

Vector meson mass fits
Since the ρ-meson mass will be needed in the chiral extrapolation as will be described below and since Eq. (3.12) involves the masses and decay constants of the vector meson states, we determine these basic quantities from the same vector-vector correlators, i.e. using the vector currents as interpolating operators for the vector mesons. Employing the same bootstrap samples as in the vacuum polarisation fits, the uncertainties of the determination of the vector meson properties can be correctly propagated to a hvp µ . In fact, our way of chirally extrapolating the data to the physical point benefits from a cancellation of the error of the ρ-meson mass achieved in this way.
In principle, the vector mesons should be treated as resonances. However, due to angular momentum conservation the decay products, two pions in the case of the light quarks, can only be produced if the kinematical condition m V ≥ 2 m 2 PS + p 2 is satisfied with p = 0. Since on the lattice with finite spatial extent L, the momenta are quantised, the above condition becomes This condition is not fulfilled for all but one ensemble (D15.48 of Table 1) where the vector meson mass m V and the energy of the 2-pion state with non-zero momentum become consistent within errors. We nevertheless treat the lightest vector meson as a stable asymptotic state for all ensembles 1 and obtain the spectral information from the large-time behaviour of the correlator projected to zero momentum in order to utilise it in the vacuum polarisation fits, including also the error propagation to a hvp µ . Following the analysis described in [23] we adopt for our correlated fits to extract m V and f V from it. The factor three arises from the polarisation sum of the vector meson. Likewise fits including M − 1 excited state contributions are performed with in an appropriate fit range. The statistical uncertainties of the fit parameters are estimated using the bootstrap method. For single-state fits the initial timeslice for the fit should be large enough such that the first excited state is sufficiently suppressed. The final fitting timeslice should be small enough to avoid the noisiest part of the correlator which is obtained from simple point sources in our calculation. Taking those restrictions into account we have selected fixed time ranges in physical units for the single-state fits by requiring the mean χ 2 /dof for all the ensembles to be close to 1. To fit the ρ-meson properties our standard fit range is 0.7 fm < t < 1.2 fm. For the ss-state we have chosen 0.9 fm < t < 1.4 fm and for the J/Ψ fits 1.2 fm < t < 1.7 fm. In the appendix in Table 4 and Table 5 we list the results for m V and f V obtained with single-state fits in our standard fit ranges. For one value of the lattice spacing the results for the ρ-meson masses are also depicted in Fig. 2. For fits with M = 2 done to check systematic effects of choosing a MNBC fit function, we have used 0.3 fm < t < 1.2 fm in the light sector, 0.35 fm < t < 1.4 fm in the strange sector, and 0.4 fm < t < 1.7 fm in the charm sector. We include the systematic effect of choosing different ranges for fitting the vector meson properties in our total error budget which will be discussed below.

Chiral extrapolation of a µ
As in [4,5] we use the modified definition to determine a hvp µ with the same motivation as outlined there. H stands for some hadronic scale that can be determined at unphysical values of the pion mass m PS . It is required to have a well-defined limit at the physical pion mass, m π , denoted H phys which has to be known from experimental measurements for example. With such a choice, by definition a hvp µ → a hvp µ when H → H phys .
Inspired by the observation that the ρ-meson gives the dominant contribution to a hvp µ , we have chosen H = m V in the following. Here, m V denotes the ρ-meson mass for unphysical values of the light quark masses. As noted in Ref. [4] the ρ-meson mass determined from the time-dependent correlator attains unphysically high values at unphysical values of the pion mass which cannot reliably be extrapolated to the physical point. This can also be seen for our N f = 2 + 1 + 1 ensembles for one value of the lattice spacing in Fig. 2.
Also in the four-flavour case the condition that the ρ-meson mass attains its physical value when the pion mass does, seems to hold. Thus, as outlined in [4],  Table 1).
amounts to an absorption of a large part of the pion mass dependence of a hvp µ in the light sector and to an effective redefinition of the muon mass on the lattice given by Since we want to use a consistent definition of the so modified muon mass for all singleflavour contributions to a hvp µ , we also use the ρ-meson mass m V to redefine a hvp µ for the strange and the charm quarks.
It is interesting to note that in this way we introduce a non-linear pion mass dependence for the heavy flavour contributions. This is in contrast to the standard definition in Eq. (2.5) where the pion mass dependence of a hvp µ,s and a hvp µ,c is rather flat. However, since the heavy flavours constitute only a small fraction of the total a hvp µ , the net effect of the modified definition of a hvp µ on the lattice is that the total a hvp µ still exhibits a flat pion mass dependence.

The light quark contribution, a hvp µ,ud
Considering only the light currents, for which the sea quark action is identical to the valence action, provides the contribution of the up and down quarks to the total a hvp µ . The pion mass dependence of this light quark portion is shown in Fig. 3. Here we utilise our lattice redefinition a hvp µ,ud , i.e. Eq.   The value at the physical point obtained by the linear fit can be compared to the value obtained with only two dynamical quark flavours from our earlier lattice QCD analysis [4] a hvp µ,ud = 5.67(11) · 10 −8 (N f = 2 + 1 + 1) a hvp µ,ud = 5.72(16) · 10 −8 (N f = 2) (4.1) yielding fully compatible results. The difference between the error of the two results is that the N f = 2 + 1 + 1 uncertainty given above is only of statistical nature whereas the N f = 2 value involves an estimate of systematic effects. The above result has been obtained by fitting all data from the ensembles listed in Table 1 simultaneously as the present quality of our data does not allow to discriminate lattice artefacts in the light sector. This is shown in Fig. 4. Here, we have first extrapolated a hvp µ,ud linearly to the physical point fixing the value of the lattice spacing. The figure shows that all chirally extrapolated values agree within the errorbars. We can therefore use a constant extrapolation to zero lattice spacing giving a hvp µ,ud = 5.72(13) · 10 −8 which is compatible with the result quoted in Eq. (4.1). We also performed a combined fit in m 2 PS and a 2 to all the data in Fig. 3 yielding a coefficient of the a 2 -term compatible with zero. Hence, for the present level of precision of our data, a hvp µ,ud does not show any significant lattice spacing artefacts.

The three-flavour contribution, a hvp µ,uds
For the three-flavour contribution we use again H = m V as the hadronic scale in order to have a consistent redefinition of the muon mass on the lattice. It turns out that this leads to larger statistical uncertainties for the strange quark contribution than employing the standard definition of a hvp µ,s . In addition, the dependence of a hvp µ,s on the squared pion mass appears to be non-linear. However, since the light quark contribution constitutes by far the largest part of a hvp µ , we still obtain a mild pion mass dependence for a hvp µ,uds as can be seen when looking at the twisted mass points (upper set of data points with filled symbols) in  In this figure we also include data obtained with different fermion actions naturally possessing differing cut-off effects from the literature. The orange downward triangles are from [24] using clover-improved Wilson fermions with two dynamical light and a quenched strange quark. The green upward triangles and diamonds have been computed with N f = 2 + 1 dynamical domain wall fermions [25]. These results from the other groups employ the standard definition for a hvp µ given in Eq. (2.5). We therefore add also twisted mass points obtained with the standard definition (lower set of twisted mass points with open symbols). We find an overall agreement between the different lattice determinations for the raw data of a hvp µ,uds when the standard definition is used. Besides the aforementioned cut-off effects, slight differences can also originate from utilising different conditions to determine the strange quark mass used for computing the strange quark contribution to a hvp µ,uds . It is clear from Fig. 5 that the improved definition of a hvp µ leads to a smooth and linear extrapolation to the physical point. In contrast, the standard definition of a hvp µ needs a more complicated extrapolation resulting in a larger uncertainty. Since the strange quark is heavier than the light quarks, the size of lattice artefacts can be significantly enhanced for the strange quark contribution. In fact, if we take lattice artefacts into account performing a combined fit in m 2 PS and a 2 , we find a non-zero value for the coefficient of the a 2 -term. The presence of lattice artefacts can also be seen when looking at the lattice spacing dependence of a hvp µ,s at a fixed pion mass of about 320 MeV as shown in Fig. 6. Here, we clearly see that the limit a → 0 can no longer be obtained by a constant extrapolation.
As a result of this observation it is clear that we have to take an a 2 -term into account when we discuss the continuum and chiral extrapolation of a hvp µ,uds . We therefore use the following fit function to obtain a hvp µ,uds at the physical point In order to compare our three-flavour value to a result from a dispersive analysis, we need to disentangle the quark flavours. To this end, we have to reweight the total a hvp µ . There are different possibilities to carry out such a reweighting. We have decided to reweight the values given in [26] by the charges of the active flavours. This approach is based on the assumption of quark-hadron duality. Given the ambiguity of such an approach we indicate by the abbreviation "pheno" that a certain phenomenological analysis has been employed. Comparing our lattice result with this phenomenological extraction method leads to a hvp µ,uds = 6.55(21) · 10 −8 (N f = 2 + 1 + 1) a hvp µ,uds = 6.79(05) · 10 −8 (pheno) (4.3) where we find, at least within the errors, an agreement. Given the fact that our phenomenological value is certainly afflicted by an unkown systematic error, we consider it reassuring that our lattice QCD analysis can reproduce the phenomenological value at this level of accuracy. As mentioned in the introduction, this ambiguity in the comparison of lattice results and those utilising the dispersion relation can be removed by the inclusion of the charm quark in the calculation which we want to report on next.

The four-flavour contribution, a hvp
µ Adding the charm quark contribution according to Eq. (3.18) again using H = m V we are able to directly compare to experimental values and those from different dispersive analyses. Since the charm quark is even heavier than the strange quark, we again use a combined fit involving an m 2 PS -and an a 2 -term of the form stated in Eq. (4.2). In this way, we arrive at the picture shown in Fig. 8. Here, our result obtained in the continuum limit and at a physical value of the pion mass, represented by the red triangle, can now be unambiguously confronted with the corresponding one from a dispersive analysis [26]:  Comparing the value of the total a hvp µ now a convincing agreement between the two ways of determining this important quantity is found. However, it needs to be noted that at this point our result from twisted mass lattice QCD has a significantly larger error than the one from the dispersive analysis.
Our result also agrees with the value a hvp µ = 6.76 · 10 −8 obtained for five flavours with the help of Dyson-Schwinger equations in [27], where the systematic uncertainty of this number has been estimated to be about 10%.

Systematic effects
Systematic effects play a very important role in any lattice calculation and need to be controlled. We therefore provide in this section a comprehensive discussion of the various systematic uncertainties appearing in our calculation.
• Finite-size effects The systematic uncertainty of finite-size effects appears to be small in our com- Furthermore, there exist ensembles to explicitly check the volume dependence of a hvp µ for a pion mass of about 320 MeV. For the B35 ensembles (see Table 1) we have two different volumes (32 3 × 64 and 48 3 × 96) at our disposal. A comparison is given in Table 2 Table 2.
Comparison of light-quark contribution to a hvp µ and total a hvp µ from ensembles of different volumes. See Table 1 for a description of the ensembles used here.
We conclude that finite size effects are negligible compared to our statistical error and we therefore do not take them as a systematic error into account.

• Chiral extrapolation
Also the systematic uncertainty of the chiral extrapolation is small. The ensembles utilised to arrive at the result given in Eq. With this procedure, we do not find any significant differences in the values for a hvp µ,s and a hvp µ,c . For the light quark contribution we find a systematic shift due to excited state contaminations when including time slices corresponding to a shift of 0.1 fm to the left of our standard fit range. This can be seen in Fig. 9.  Additionally, we have checked that varying the matching momentum between 1 GeV 2 and 3 GeV 2 gives compatible results for a hvp µ as long as the transition between the fit functions in the low-and high-momentum regions is smooth. Another criterion has been that the coefficients of the fit polynomial in the high-momentum region, where more data is available, do not influence the coefficients of the fit polynomial in the low-momentum region. Applying these criteria all choices of different functions to combine the two momentum regions have not resulted in significant differences in the final values for a hvp µ compared to the Heaviside step function that we have used for our quoted result.
In [10] it was suggested to utilise so-called correlated [1,1] Padé approximants to fit the vacuum polarisation function in the low momentum region. We tried such fits using an upper integration scale of 1.5 GeV 2 . Compared to our standard M1N2 fit, with the same upper integration limit, we obtain compatible results for the light a hvp µ,ud employing the standard definition of Eq. (2.5). Performing Padé fits for the redefinition in Eq. (3.18) using H = m V determined from the vector meson twopoint function we, however, observe a larger error than from our M1N2 fits. This is caused by the fact that for the Padé fits the pole is fitted simultaneously and not taken as an external input parameter and thus the favourable error cancellation between the pole and the vector meson mass which holds for the M1N2 fits for the redefinition of Eq. (3.18) with H = m V does not occur. In fact, we find that fitting the vacuum polarisation function by the Padé ansatz the values of a hvp µ,ud obtained with the redefinition show up to three times larger uncertainties.

• OS matching uncertainties
Since we use the OS action in the strange and charm quark sector, different values for the corresponding quark masses could be used which, however, have to lead to the physical values of the Kaon and D-meson masses in the continuum and chiral limit. Varying the strange and charm quark masses within the uncertainties given for aµ s and aµ c in Sec. 3.1 has been found to be negligible. Likewise changing µ s to the value obtained from directly matching with the physical kaon mass gives a compatible result. The same is true when using the µ s and µ c values procured when allowing for a 2 -effects in the fit function employed in the matching.
• Different strange and charm sea quark masses Additionally to the choice of valence quark masses, our result might be influenced by sea quark masses which for some of the ensembles have not been tuned to their correct physical values. For details see [8]. By changing the mass splitting parameter µ δ of the twisted mass action for a non-degenerate fermion doublet Eq.   Table 3.
Comparison of single-flavour contributions and total a hvp µ from ensembles having different strange and charm sea quark masses.
their physical values. The new ensemble is called A100.24s whereas the old one is A100.24, sharing apart from µ δ the same parameters. They have both been tuned to maximal twist and possess a pion mass of about 500 MeV and a space-time volume in units of the lattice spacing of 24 3 × 48. Due to the large pion mass they are not included in the rest of our analysis. Using the same matching condition for the OS valence quarks, we arrive at consistent values for the single-flavour contributions to a hvp µ as can be seen in Table 3. Hence, we conclude that the impact of different sea quark masses in the heavy sector on a hvp µ is negligible.

• Disconnected contributions
This is a systematic effect we can currently not adequately quantify. There are, however, several reasons for assuming that the disconnected contributions are small. First of all, a dedicated study in the 2-flavour case has revealed them to be compatible with zero [4]. Note, however, that in Refs. [28,29] the impact of disconnected contributions has been estimated to be −10%, at least in the energy range 2m π < q < 3m π . Secondly, in the SU (3) flavour limit they are identically zero due to charge cancellation. Thirdly, the disconnected contribution arising from the charm quark has been computed in perturbation theory and shown to be suppressed by a factor

Conclusions
In this paper, we have performed a first calculation of the four-flavour leading-order QCD contribution to the muon anomalous magnetic moment. Such a four-flavour computation has the invaluable advantage that a direct comparison to a phenomenological extraction of this quantity can be achieved without any ambiguity when discriminating the different flavour contributions which appears unavoidably in a two-or three-flavour comparison.
In our work we found that also for the four-flavour situation considered here, the improved method of Ref. [4] leads to a smooth and well-controlled chiral extrapolation of a hvp µ . In addition to the chiral extrapolation, we performed a comprehensive analysis of systematic effects, such as finite lattice spacing and finite volume artefacts, excited state contaminations in the vector meson mass determination, different choices of the vacuum polarisation fitting function, and different choices of valence and sea strange and charm quark masses. From this set of systematic effects only the different choices of the fitting functions and excited state contaminations of the vector states lead to significant systematic uncertainties which we include in our final error estimate.
As our main result we provide a comparison to a dispersive analysis [26]: In Fig. 11 we also compare the outcome of our first-principle computation with a summary of other results obtained utilising the dispersion relation. Although our lattice QCD determination of a hvp µ shows an overall agreement with phenomenology, the lattice QCD result has clearly a significantly larger error being, however, already at the same order of magnitude.

HLS estimate
This work a hvp µ 7.4e-08 7.2e-08 7e-08 6.8e-08 6.6e-08 Figure 11. Comparison of our first four-flavour lattice result of a hvp µ with different results based on dispersion relations: Davier et al. [31], Jegerlehner and Szafron [32], Hagiwara et al. [33], and HLS [34] A substantial step forward to improve the lattice determination will be computations directly at the physical value of the pion mass. Such simulations, in combination with a significantly increased statistics including also isospin breaking and electromagnetic effects have the potential to reach or even beat the error from a dispersive analysis. In addition, results from lattice QCD will only rely on QCD alone and hence can provide a stringent test of the standard model. This opens the exciting possibility to eventually clarify whether the present discrepancy for the muon anomalous magnetic moment is indeed a sign of new physics.