Lattice calculation of the leading strange quark-connected contribution to the muon g − 2

We present results for the leading hadronic contribution to the muon anomalous magnetic moment due to strange quark-connected vacuum polarisation effects. Simulations were performed using RBC-UKQCD’s Nf = 2 + 1 domain wall fermion ensembles with physical light sea quark masses at two lattice spacings. We consider a large number of analysis scenarios in order to obtain solid estimates for residual systematic effects. Our final result in the continuum limit is aμ(2)had,s = 53.1(9)(− 3+ 1) × 10− 10.


Introduction
With an accuracy of the order of 1ppm [1], the anomalous magnetic moment of the muon, a µ , is one of the most precisely determined quantities in experimental as well as theoretical particle physics. Since the value of a µ is sensitive to potential new physics contributions (see e.g. [2]), the persistent 3 − 4 σ tension between experiment and theory is generating much interest within the particle physics community. The new g−2 experiment at Fermilab g − 2 (E989 collaboration) is expected to reduce the uncertainty in the experimental value by a factor of around four, down to 140ppb [3]. This puts much pressure on the theory community to match this precision. Table 1 lists our current knowledge of the Standard Model (SM) contributions to a µ [1]. The dominant source of the uncertainty comes from the leading order (LO) hadronic contribution, a (2)had µ , which we concentrate on here. This is followed closely by the lightby-light (LbL) contribution, a (3)had µ , for which we note the recent efforts in the computation of this quantity on the lattice [4][5][6].
The SM prediction for the LO hadronic contribution as stated in by relating the photon hadronic vacuum polarisation (HVP) to the cross section data for e + e − decays into hadrons using a dispersive integral over this data [7,8]. Despite providing an accurate determination of a (2)had µ , there are underlying difficulties in interpreting the cross section data and in combining individual data sets to yield the final result. Given the importance of a µ in building models of new physics above the electroweak scale, an entirely independent computation of a (2)had µ from first principles is highly desirable. As will be explained in detail in the next section, the basic building block of the lattice computations of a (2)had µ [9][10][11][12][13][14][15][16][17][18][19] is the 2-point correlation function of two electromagnetic currents. This splits into connected and disconnected Wick contractions which, as was argued in [20,21], all have their individual infinite volume and continuum limits. It has recently become increasingly apparent that tailor-made techniques in lattice QCD have to be devised for individual Wick contractions in order to achieve the required level of precision. Recently, for example, the full set of quark-disconnected contributions, so far believed to be the main obstacle in obtaining a lattice determination of a (2)had µ at the percent level, has been computed with much improved precision in lattice QCD with physical light quark masses [17]. The crucial step in this computation was identifying the dynamics mainly responsible for the disconnected contribution and tailoring corresponding lattice techniques [17,22].
In this spirit we here present results for a second building block toward the full lattice computation of a (2)had µ , namely the computation of the quark-connected strange contribution, a (2)had,s µ , in the continuum limit of Möbius domain wall fermion (MDWF) [23][24][25][26][27] lattice QCD with N f = 2 + 1. We use a variety of analysis techniques in order to test both the techniques themselves and their effect on the final value of a (2)had,s µ . In section 2 we discuss the lattice strategy for computing the HVP form factor, which we motivate as a crucial ingredient in the computation of a (2)had µ . In section 3 we present details of the data analysis techniques we have used and our final results for a

JHEP04(2016)063
2 Lattice computation of the HVP form factor Before describing our computation of the HVP form factor, it is worth motivating this computation. The contribution a (2)had µ can be related to the Euclidean space-time HVP in the following way [9]: where α is the QED coupling, Q is the Euclidean four-momentum of the intermediate photon,Π(Q 2 ) = Π(Q 2 ) − Π(0) is the renormalised HVP form factor, and f is the following integration kernel: and m µ is the mass of the muon. The HVP form factor is related to the electromagnetic current 2-point function in momentum space through the usual form factor decomposition where δ µν is the Euclidean metric. The HVP is therefore crucial in the computation of a (2)had µ .

General lattice methodology
To compute the quark-connected HVP form factor on the lattice, we choose the following discrete version of the electromagnetic current 2-point function: where V f µ is the conserved vector current for some choice of lattice action, V f ν =q f γ ν q f is the non-conserved local vector current, the subscript f indexes quark flavours, Q f is the electric charge of flavour f in units of the positron charge, and Z V is the renormalisation constant for V f ν . For this specific choice of currents, one obtains a Ward identity similar to the continuum one where ∂ * µ is the backward finite difference operator. This identity guarantees that the shortdistance divergences in C µν (x) for x → 0 are at most logarithmic and can be regulated by using the usual subtracted form factorΠ(Q 2 ) = Π(Q 2 ) − Π(0).
In practice, we evaluate the current 2-point function using either a point source (i.e. as in eq. (2.5)) or a complex-valued Z 2 wall source [28][29][30], which performs a stochastic average JHEP04(2016)063 on the local current spatial position. Noise sources have been known to improve the signalto-noise ratio over point sources at the same computational cost [30], and we will provide such a comparison in the next section.
We also employ the following modification of the definition of the HVP tensor in eq. (2.3), where a is the lattice spacing. The additional term on the right hand side corresponds to a zero-mode subtraction (ZMS) [31]. In the infinite volume theory one can show that this zero-mode vanishes by Lorentz symmetry. However, in finite volume, where momentum is discretised, the volume sum of C µν does not have to vanish. 1 As will be discussed later, we find that this procedure greatly improves the signal-to-noise ratio for Π µν (Q), in particular at low-Q 2 . On a finite lattice Lorentz symmetry is broken into a finite symmetry group. As a result, the tensor decomposition in eq. (2.4) receives additional contributions whereQ = 2 a sin(aQ/2) is the usual lattice momentum. The ellipsis denotes a series of terms individually proportional to a product of Q n µ Q m ν for some odd integers n and m and µ Q n µ with n an even integer. These contributions are hyper-cubic covariant expressions that are not Lorentz covariant: they vanish in the simultaneous continuum and infinite volume limit where Lorentz symmetry is restored. Contributions containing Q n µ Q m ν are sensitive to the anisotropy of the momentum Q and can be removed exactly by only considering momenta where Q µ = 0 or Q ν = 0 [13]. In all the following, we will only consider momenta with a vanishing spatial part, and we define our lattice HVP form factor function as follows: where the index j runs over spatial directions only.

Ensemble properties
We present results on two dynamical ensembles with near-physical quark masses and 2+1 dynamical flavours of domain wall fermions (DWF) [23,24]. Our formulation of DWF uses a Möbius action with an H T kernel to improve the sign function approximation as described in [25][26][27], and we hence refer to this formulation as MDWF. The nice property of this choice of discretisation is a continuum-like chiral symmetry, which produces automatic O(a)-improvement. The explicit form of V f µ for this action can be found in [34]. The ensembles, which are described in detail in [34], have been generated with the Iwasaki gauge action, and their basic properties are listed in   parameters t 0 and w 0 , the inverse lattice spacing was computed for these ensembles using as hadronic input the masses of the pion, kaon and omega baryon [34]. As indicated by the kaon masses in table 2, which deviate from the value of 495.7 MeV taken as the target physical value in [34], each of the ensembles we have used has slight mistunings in the masses of the strange quarks in the sea. To account for this we performed two sets of strange measurements on each ensemble: one unitary and one partially quenched. A summary of our measurements can be found in table 3.

Comparative study of point and stochastic sources
For our valence measurements we again used MDWF. We initially performed inversions on both stochastic Z 2 wall and point sources. We accelerated our inversions using the HDCG algorithm [35]. For our unitary measurements on the 48I ensemble we performed inversions using Z 2 wall sources on every other timeslice, making 48 measurements per configuration, whilst for point sources we performed inversions on every eighth timeslice, making 12 measurements per configuration. In the point source case we located the source at the spatial origin of each timeslice. A similar set of measurements for the 64I ensemble can also be found in  Table 4. Values of point (Q 2 min )/ Z2 (Q 2 min ) under various analysis conditions as computed on the 48I ensemble, where is defined in equation (2.10). Here Q 2 min refers to the lowest non-zero value of Q 2 . Note that the Z 2 wall sources only provide an improvement over point sources for the same computational cost (i.e. the first row of this table) when the ZMS procedure is applied. as defined by where ∆Π(Q 2 ) denotes the statistical error in Π(Q 2 ). We compared this quantity for the two different source types at the lowest non-zero Q 2 , which we denote Q 2 min . As will become clear later, this is the region that contributes predominantly to a (2)had,s µ due to the diverging nature of f (Q 2 ) as Q 2 → 0. We also compared the effect of the ZMS technique on the error at the smallest non-zero Q 2 . Table 4 shows the factors of improvement of the Z 2 wall source data over the point source data, as well as the effects of ZMS and the number of timeslices used. ZMS allows Z 2 wall sources to out-perform point sources in the low-Q 2 region, reducing Z 2 (Q 2 min ) by a factor of about 87 in the equal cost case on the 48I ensemble. For this reason the remainder of this paper will use results exclusively from our measurements on Z 2 wall sources.
3 Computation of a (2)had,s µ In this section we describe how we compute a (2)had,s µ from the HVP form factor discussed in the previous section. We begin by describing two strategies for performing the integral in equation (2.1), namely the hybrid method and sine cardinal interpolation (SCI). This is followed by a description of our continuum and quark mass extrapolations. We conclude by summarising our systematic error estimation and presenting our final result.

Hybrid method
The I = 1 contribution to equation (2.1) is highly peaked near Q 2 ∼ m 2 µ /4, and a (2)had,s µ is expected to be similarly dominated by contributions from the low-Q 2 region. This presents a challenge for any lattice computation of this quantity, since lattice momenta are generally quantised due to the imposition of a finite volume with periodic boundary conditions. In this particular case we are restricted to Q 0 = 2πn 0 T , where n 0 is an integer and −T /2 ≤ n 0 < T /2. This means the lowest non-zero Q 2 we can achieve with the two ensembles available to us is approximately 0.013 GeV 2 , or approximately 1.2m 2 µ . We must hence employ some parametrisation or model to approximate Π(Q 2 ) at small Q 2 .
To this end, we use the hybrid method as described in [36]. This method consists of partitioning the integrand in equation (2.1) into three non-overlapping adjacent regions using cuts at low-and high-Q 2 (Q 2 low and Q 2 high ) (see figure 1). The integrand is then computed JHEP04(2016)063

Numerical Integration
Model Perturbation Theory for the three regions in different ways. The low-Q 2 region is integrated by constraining some parametrisation of Π(Q 2 ) using the data computed on the lattice. This parametrisation is then used to compute Π(0) and thenceΠ(Q 2 ). This result is then combined with the kernel f (Q 2 ) to produce the integrand of interest, which is then integrated numerically. The mid-Q 2 region is integrated directly by multiplying the lattice data by f (Q 2 ) before using some numerical integration method such as the trapezium method. Finally, the integral over the high-Q 2 region is computed using perturbation theory. This last step is performed by using the 3-loop expression [37] for the HVP form factor combined with our previous result [34] for the strange quark mass in the MS scheme at 3 GeV. The three-loop expansion of the perturbative expression is more than adequate for the purposes of the present calculation, since the higher order corrections are negligible in this context, and the perturbative contribution typically accounts for 0.1% of the value of a (2)had,s µ . When performing the integration of the mid-Q 2 region, if either of the specified values of Q 2 low and Q 2 high is not aligned with any values computed using the lattice, then a simple linear interpolation is performed to compute a value ofΠ(Q 2 low ) orΠ(Q 2 high ). Using the hybrid method, we can minimise systematic errors arising from the use of a parametrisation of Π(Q 2 ) when extrapolating to Q 2 = 0. The magnitudes of the curvatures in Π(Q 2 ) decrease monotonically with increasing Q 2 , whilst the statistical errors in Π(Q 2 ) decrease. There is hence an incentive to reduce Q 2 low in order to minimise the systematic error arising from the use of a parametrisation for the HVP. However, there is also an incentive to increase Q 2 low to increase the amount of data available to the parametrisation, improving the statistical error on the low-Q 2 integral. The dispersive model study of [36] provides some useful guidance on the selection of Q 2 low . In response we have performed our analysis with various Q 2 low in an attempt to ascertain the effect of varying this parameter on the final result. Based on dispersive model studies of the type recommended in [38],

JHEP04(2016)063
all of the fits entering our final assessment, with the exception of the R 1,1 Padé fits where Q 2 low = 0.7 and 0.9 GeV 2 , would be acceptable for use in the isovector channel. With the strange HVP form factor exhibiting significantly less curvature than the light quark HVP form factor, larger Q 2 low will be usable for a given parametrisation in the strange case. As we will show, the excellent agreement of the results obtained from the R 1,1 Padé fits where Q 2 low = 0.7 and 0.9 GeV 2 with those of the other fits confirms this expectation.
We use a variety of parametrisations to integrate the low-Q 2 region in the hope of determining the systematic uncertainty arising from this method. In addition, we have used two methods to constrain these parametrisations. We discuss these aspects of our analysis in the following subsections.

Low-Q 2 parametrisations
We use two classes of parametrisations for the low-Q 2 region when performing the integral in equation (2.1): Padé approximants and conformal polynomials. It has been shown that both of these representations of the HVP converge to the HVP as successive terms are added to them [36,39]. In this sense they are independent of any phenomenological model.
The Padé approximants are motivated by the once-subtracted dispersion representation of the HVP [40], i.e.
where ρ(t) is the vector spectral density. Using a Stieltjes transformation it can be shown that Φ can be expressed as a continued fraction of Stieltjes functions. This function can in turn be approximated by Padé approximants that converge to Φ(Q 2 ) as more values of Q 2 and Φ(Q 2 ) are used to constrain the Padés. The Padés have poles on the negative real axis, and so we choose to write them as follows [40]: where a i , b i , Π 0 and c are parameters to be determined. The dispersive model study of the I = 1 contributions in [36] suggests that, for the Q 2 low we intend to work with, the R 1,1 and R 1,2 forms will provide an accuracy below ∼ 1%.
The conformal polynomials are motivated by a desire to improve the convergence properties of the Taylor series of Π, which is only convergent for Q 2 less than the square of the two-particle mass threshold, E min . We employ the standard conformal transformation approach to map the Q 2 -plane onto the unit disc, i.e. we introduce where E is some energy parameter with the requirement E < E min . This results in the Q 2 -plane, excluding the real interval (−∞, −E 2 ), being mapped onto the interior of the unit disc, with the interval (−∞, −E 2 ) being mapped onto its boundary. Provided that

JHEP04(2016)063
E remains below the two-particle mass threshold, a Taylor expansion of Π in w will be convergent for Q 2 ≥ 0. Our truncated conformal polynomial ansätze of degree n are hence described by where p k and Π 0 are parameters to be determined. Drawing on [36], we expect third-and fourth-order polynomials to be adequate in describing the lattice data at low-Q 2 .

Matching at low-Q 2
We use two techniques to constrain the low-Q 2 parametrisations: a fit using χ 2 minimisation and discrete time moments. The latter is a discrete version of the moments method described in [18]. The χ 2 minimisation involves a fit where the covariance matrix is approximated by its diagonal, i.e. the fit is uncorrelated. In principle different values of the HVP form factor are strongly correlated because they originate from the same data. In practice we found the correlated fit impossible to perform, the covariance matrix being singular at the present level of precision. Further to this, we also found that the eigenvalue spectrum of the covariance matrix did not allow for the elimination of singular values from the matrix whilst preserving the essential information contained within it. We hence found that using the pseudo-inverse of the covariance matrix was not advantageous when compared to replacing the covariance matrix with its diagonal. The χ 2 minimisation lends weight to points in the computed HVP with a smaller statistical error at larger values of Q 2 .
The moments method exploits the relationship between the HVP form factor and the diagonal components of the lattice space-averaged current-current correlator Taking the 2n-th discrete central derivative in direction 0 at Q 0 = 0 allows us to writē where∂ Q 0 is a general central discrete derivative operator. In this particular analysis we use a central discrete derivative improved to O(a 2 ). We then insert one of the above analytical ansätze for the HVP form factor, setting up a system of non-linear equations that we solve numerically to determine the ansatz parameters. When performing the moments method we use a representation of the HVP that is a function ofQ 2 . However, within the moments method, derivatives are taken with respect to the Fourier momentum Q 0 and notQ 0 . We observed a marked reduction in the cut-off dependence of a (2)had,s µ in response to this change in momentum definitions. Within the determination of the ansatz parameters, the low-Q 2 cut is not used as an input for this technique, so the resulting parameters do not depend on the low cut used in the hybrid method [18]. Resulting parametrisation after matching parametrisation R 1,1 using a χ 2 fit. This curve is typical of the parametrisations generated using the various analytical expressions and matching methods described in this paper. We find that these results typically pass within negligible distance of the lattice data point central values. Figure 2 shows a typical parametrisation resulting from the techniques and parametrisations described above. The HVP data in these plots is computed on the 48I ensemble using the unitary strange quark masses. We find that both matching techniques produce parametrisations that differ negligibly from the lattice Π(Q 2 ) data for Q 2 ≤ Q 2 low .

Integrating the low-and mid-Q 2 regions
The numerical evaluation of (2.1) is problematic, as the integrand is highly peaked near Q 2 = 0. To overcome this difficulty we perform a change of variables which allows us to rewrite the low-and mid-Q 2 portions of the integral as An example of the resulting integrand is given in figure 3. In this case an R 11 parametrisation was used and the matching was performed using discrete moments with a low-Q 2 cut of 0.7 GeV 2 . This figure highlights the peak in the low-Q 2 region, which can significantly affect the final value of a (2)had,s µ if it is poorly constrained.

Sine cardinal interpolation
One alternative to the hybrid method is computing the HVP directly at an arbitrary momentum by performing the Fourier transform in equation (2.7) at said momentum [41]. Whereas before we used Q 0 = 2π T n 0 with n 0 ∈ Z, −T /2 ≤ n 0 < T /2, we now let . Compared to this equation the integrand in the plot has been multiplied by the factor of 4α 2 for consistency with eq. (2.1). The parametrisation is achieved using discrete moments to constrain R 1,1 . The red lattice data points are computed using unitary strange data on the 48I ensemble. Note that, despite the legend, the blue curve has not actually been fitted directly to the lattice data points. Rather, the HVP parametrisation has been constrained before multiplying it with the integration kernel in eq. (2.2).
anywhere on the real half-closed interval [−T /2, T /2). This allows for the computation of a (2)had,s µ without using a parametrisation of the HVP. Because of its connection to sampling theory [33], we call this technique sine cardinal interpolation (SCI). This interpolation of the discrete value of the HVP tensor in the calculation of a (2)had,s µ is a source of finite-time effects, which can be shown to decay exponentially with m π T [33].
Using this technique, we compute the HVP at arbitrary momenta up to Q 2 high , after which the perturbative result is used. To compute a (2)had,s µ from (2.1), the integration up to Q 2 high is performed in a similar way to what is described in section 3.1.3.

Physical mass and continuum extrapolations
We extrapolate to both the continuum limit and the physical strange quark mass using the values of a where am res is the residual mass arising from residual chiral symmetry breaking in MDWF, and am phys s is the lattice strange quark mass required to give the target kaon mass for the ensemble in question, as specified in [34] and table 2. Because we are using the MDWF action, which is O(a) improved, we can neglect cut-off effects of this order when extrapolating to the continuum limit. To account for errors in the physical value of the strange quark mass, we use a Gaussian distribution to sample this value for each ensemble using  the error specified in [34] and table 2. We perform a correlated fit using the four values of a (2)had,s µ computed from our two ensembles in the unitary and partially quenched theories. We also attempted a physical point extrapolation where we forced the value of α in eq. (3.9) to equal zero, meaning we performed a constant fit in a 2 . We found that it was not possible to exclude this ansatz on the basis of the resulting χ 2 or p-value. However, there is no theoretical justification for the absence of an a 2 dependence within a (2)had µ for MDWF. On this basis, and since the constant fit with α = 0 could artificially decrease the error in the extrapolated value of a (2)had,s µ , it is necessary to include the a 2 term in our fit ansatz. Figure 4 illustrates examples of our continuum and strange quark mass extrapolations. In the left-hand plot the lattice data has been projected into the physical strange quark mass limit, meaning we have subtracted variations arising from the strange quark mass. In the right-hand plot, we have projected the lattice data into the continuum limit in a similar manner. To produce these particular plots we used the P 0.6GeV 3 parametrisation, which was constrained using discrete moments. The low cut in this case was 0.7 GeV 2 . We found a strong dependence of a (2)had,s µ on the strange quark mass, to the extent that the sign on α changed in response to the inclusion of the partially quenched data points (see figure 4). This had the effect of shifting the final value of a (2)had,s µ from approximately 50 × 10 −10 to 53 × 10 −10 .

Statistical error propagation
This analysis relies on various measurements computed as part of global chiral fits to results on a number of different DWF ensembles [34]. Of particular note is the lattice JHEP04(2016)063 spacing, which is required to reconcile the dimensionful muon mass with the dimensionless lattice momenta used in the integration kernel f . In order to account for potential non-Gaussianity, this was sampled from the global fits jackknife samples used in [34]. We found that the inclusion of the lattice spacing error increased the error in the final value of a (2)had,s µ significantly, since the peak in the integrand (see figure 3 for example) depends strongly on the muon mass.
In addition, for Z V we drew random samples from a Gaussian distribution for each bootstrap sample. Since the statistical error on Z V is small (0.04% for the 48I ensemble and 0.02% on the 64I ensemble), we assume the original data set follows a Gaussian distribution.

Systematic error estimation
We use a variety of analysis techniques in order to determine the systematic error in the value of a (2)had,s µ arising from the choice of a particular technique. Although different in some aspects, this method is motivated by the frequentist approach developed in [42].
We picked energy thresholds of 0.5 and 0.6 GeV for the chosen conformal polynomials as we believed these to be below the two particle energy threshold, and we wished to study the effect of the variation of this quantity on the final value of a (2)had,s µ . The Padé approximants and the conformal polynomials have been shown to converge to the HVP in the limit of infinitely many parameters [36,39]. We observed that the result for a (2)had,s µ underwent a saturation as more terms were added to these parametrisations, although only at the level of the statistical error. We took this as a possible manifestation of the aforementioned behaviour. As a result, we chose to rely only on the two higher order parametrisations to approximate the low-Q 2 region. These are expected to be closest to the physical value and agree well with the recommendations of [36].
We used three different low cuts: 0.5 GeV 2 , 0.7 GeV 2 and 0.9 GeV 2 . These were selected such that we had sufficient degrees of freedom to perform a χ 2 fit for all the parametrisations described above. We initially experimented with three high cuts: 4.5 GeV 2 , 5.0 GeV 2 and 5.5 GeV 2 . We selected these at 0.5 GeV 2 spacings to allow sufficient variation in the cut so that the perturbative contribution could vary. However, it became apparent that the high cut made negligible difference to the final value of a (2)had,s µ (less than 0.1% of the final value), and so ultimately we chose a single high cut at 5.0 GeV 2 .
We also varied the numerical technique used to integrate the mid-Q 2 region when implementing the hybrid method. We studied the effect of using the trapezium rule and Simpson's rule.
Finally, we used both discrete time moments and a χ 2 minimisation to determine the extent to which the low-Q 2 matching technique affected a

JHEP04(2016)063
In the case of sine cardinal interpolation we used a step in n 0 of 0.005, with the same high cut as used in the hybrid method (5.0 GeV 2 ). We found this step size was sufficient to produce a value of a (2)had,s µ with an integration error that was negligible compared to our statistical error.
In total we used 73 different methods to determine a (2)had,s µ . We display stacked histograms of these values in figure 5, colour coded according to which aspect of the analysis is being varied. These various values enable us to gauge the systematic error arising from our choice of analysis technique. We compute the overall central value by taking the median of the central values from each of the 73 analyses and take the statistical error as being the bootstrap error for the analysis corresponding to this value. The systematic error is then computed by taking the difference between this central value and the smallest and largest of the 73 analysis central values. This gives us an asymmetric determination of the systematic uncertainty in the final value. From panel (a) in figure 5, it is apparent that much of this asymmetry comes from the P 0.5GeV 3 parametrisation constrained using a χ 2 fit. One feature immediately apparent in figure 5 is the apparent lack of sensitivity of the final value of a (2)had,s µ to a particular analysis technique, especially when compared to the overall statistical error. Indeed, the only set of analyses that could be considered outliers is approximately 0.25σ from the band of central values around 53 × 10 −10 .
We find that the values of a (2)had,s µ computed using the discrete moments matching method are more consistent with one another than those computed using a χ 2 fit. There are two reasons for this. First, the discrete moments method does not depend on the value of the low-Q 2 cut, meaning that the parameters for a particular parametrisation will be the same as the low cut is varied. Second, the moments method relies on expanding the HVP parametrisation as a Taylor series around Q 2 = 0. As a result, the parameters are more sensitive to variations in the HVP at low-Q 2 . This can be contrasted with the χ 2 fit strategy, where the points at larger Q 2 have a smaller statistical error and so contribute more to the χ 2 , playing a larger role in constraining the low-Q 2 parametrisation than those at small Q 2 . This is not to say that the moments will always produce an excellent parametrisation of the HVP at all Q 2 , but given that the integrand in eq. (2.1) is highly peaked at low Q 2 , any deviation from the true HVP at large Q 2 by one of these parametrisations will be suppressed by the integration kernel f .
Finally, the central value of a (2)had,s µ computed using SCI shows good agreement with those computed using the other analysis methods.
We expect finite-volume (FV) effects to be very small for the strange HVP, for the following reason. Although NLO chiral perturbation theory (ChPT) does not provide a good low-Q 2 representation of the fully subtracted HVP form factor, its two-pion loop contribution has recently been shown to reproduce observed FV effects rather well [32]. This observation is not totally unexpected, since contributions from the lowest-lying states (in this case, two pions states) are expected to dominate FV effects, and such contributions are present already at NLO. This is in contrast to the resonance contributions which, (2)had,s µ computed using the various analysis techniques, colour coded by the low-Q 2 parametrisation (panel (a)), the method used to match these parametrisations (panel (b)), the low-Q 2 cut (panel (c)) and the numerical method used to integrate the mid-Q 2 region (panel (d)). The large grey band illustrates the final statistical error in our result. current and its subsidiary strange component from coupling to two pions in the isospin limit. As a result, two-pion-induced FV effects are absent, for example, from the full connected plus disconnected strange current contribution to the LO HVP. This is not to imply that two-pion FV effects are generally negligible; though absent from the isoscalar contribution to the LO HVP, they are certainly present in its isovector contribution, and have to be dealt with there. We therefore expect the leading finite volume effect in the strange case to be negligible as a result of an exponential suppression like e −m K L , where m K L ≈ 13.8 for our ensembles. This situation is entirely different from the case of the light contribution, where pion-induced FV effects, which are expected to have an e −mπL rather than e −m K L suppression, will be significantly larger, and appear non-negligible for the volumes currently available to us.

Final results
Following the procedure outlined in section 3.4.2, gives us our final result: a (2)had,s µ = 53.1(9)( +1 −3 ) × 10 −10 , (3.10) where the first error is that arising from our statistical uncertainty and the second is that arising from our systematic error determination. The central value and statistical error here correspond to the analysis using a P 0.5 GeV 4 low-Q 2 parametrisation constrained using the discrete moments method with a low cut of 0.5 GeV 2 and a high cut of 5.0 GeV 2 . In this case Simpson's rule was used for the mid-Q 2 region. If we were to omit the P 0.5GeV n parametrisations from the group of analyses used to determine the systematic error, we would expect it to be much more symmetric: a (2)had,s µ = 53.1(9)(1) × 10 −10 . (3.11) In both cases we find that our overall error (approximately 2%) is dominated by our statistics, which illustrates the robustness of the various analysis techniques. This uncertainty is small enough to allow for a future evaluation of the total a (2)had µ with sub-percent precision. In addition, our final value is in good agreement with HPQCD, who quote 53.4(6) × 10 −10 as their final value 2 [18].

Conclusion
We have computed the quark-connected strange contribution to the anomalous magnetic moment of the muon using scaled Shamir domain wall fermions with physical quark masses. We have used a variety of analysis techniques, most notably the hybrid method, and a variety of low-Q 2 parametrisations in order to gauge the systematic uncertainty arising from the selection of any particular analysis technique. Our use of the hybrid method allows us to overcome the systematic effects associated with using a low-Q 2 parametrisation of the HVP at values of Q 2 large enough to use the perturbative expression. We have focused on using Padé approximants and conformal polynomials for our low-Q 2 parametrisations, since these are well-motivated and model independent. These and other variations in our analysis allow us demonstrate the insensitivity in the final value of a (2)had,s µ to these variations. Our final result, as stated in equation (3.10), is in good agreement with the value quoted by HPQCD [18]. Furthermore, the final error in our result is well within the limits required to produce a value of a (2)had µ to rival that produced by current phenomenological methods. Our research into the computation of the light contribution to a (2)had µ is ongoing. Once again we expect that a computational strategy will need to be tailored to reduce the statistical error in this result and put us in a position to compete with the phenomenological 2 HPQCD use a slightly different definition of the isospin symmetric kaon mass, which differs from ours due to electromagnetic effects. We also performed our extrapolation using the same convention as HPQCD. We observed a relative deviation of a fraction of a percent, which is compatible with the O(α) effects expected from this change in convention and represents an insignificant per mille correction to the total value of a . Beyond this, there will be the eventual need to include isospin breaking effects in our results. Finally, our results for the connected HLbL contribution computed at physical pion mass are encouraging [43], and studies of disconnected contributions and finite volume and non-zero lattice spacing effects are underway.