A novel computational paradigm for a precise determination of the hadronic contribution to (g_μ-2) from lattice QCD
Leonardo Giusti, Mattia Dalla Brida, Tim Harris, Michele Pepe
AA novel computational paradigm for a precisedetermination of the hadronic contribution to ( 𝒈 𝝁 − ) from lattice QCD Leonardo Giusti, 𝑎,𝑏, ∗ Mattia Dalla Brida, 𝑎,𝑏
Tim Harris 𝑎,𝑏 and Michele Pepe 𝑏 𝑎 Dipartimento di Fisica, Università di Milano-BicoccaPiazza della Scienza 3, I-20126 Milano, Italy 𝑏 INFN, Sezione di Milano-BicoccaPiazza della Scienza 3, I-20126 Milano, Italy
E-mail: [email protected], [email protected],[email protected], [email protected]
The hadronic contribution to the muon anomalous magnetic moment 𝑎 𝜇 = ( 𝑔 𝜇 − )/ 𝑎 𝜇 . ∗ Speaker © Copyright owned by the author(s) under the terms of the Creative CommonsAttribution-NonCommercial-NoDerivatives 4.0 International License (CC BY-NC-ND 4.0). https://pos.sissa.it/ a r X i v : . [ h e p - l a t ] J a n L-HVP
Leonardo Giusti
1. Introduction
The measurement of the muon anomalous magnetic moment 𝑎 𝜇 = . ( . ) × − by the E821 experiment has the remarkable precision of 0 .
54 parts per million (ppm) [1], and theon-going E989 experiment at FNAL is expected to reach the astonishing precision of 0 .
14 ppm bythe end of its operation [2]. The Standard Model (SM) prediction includes contributions from five-loop Quantum Electrodynamics, two-loop Weak interactions, the Hadronic leading-order VacuumPolarization (HVP) and Hadronic Light-by-Light scattering (HLbL) [3]. The overall theoreticaluncertainty is dominated by the hadronic part. So far, lacking precise purely theoretical computa-tions, the hadronic contributions have been extracted (by assuming the SM) from experimental datavia dispersive integrals (HVP & HLbL) and low-energy effective models supplemented with theoperator product expansion (HLbL). This leads to 𝑎 𝜇 = . ( . ) × − (0.37 ppm) [3],which deviates by 3 − .
6% to roughly 2%, seeRef. [3] and references therein, corresponding to an overall error on 𝑎 𝜇 which is still 5-15 timeslarger than the anticipated uncertainty from E989. The main bottleneck [3] for matching thatprecision is the large statistical error in the Monte Carlo evaluation of the required correlationfunctions. We have recently proposed [4] to solve this problem by a novel computational paradigmbased on multi-level Monte Carlo integration in the presence of fermions [5, 6]. With respect tothe standard approach, this strategy reduces the variance exponentially with the temporal distanceof the fields. In this first feasibility study we focused on the HVP, but the strategy is general and itcan be applied to the HLbL, the isospin-breaking and electromagnetic contributions as well.
2. The signal-to-noise problem
The HVP can be written as 𝑎 HVP 𝜇 = (cid:16) 𝛼𝜋 (cid:17) ∫ ∞ 𝑑𝑥 𝐾 ( 𝑥 , 𝑚 𝜇 ) 𝐺 ( 𝑥 ) , (1)where 𝛼 is the electromagnetic coupling constant, 𝐾 ( 𝑥 , 𝑚 𝜇 ) is a known function increasingquadratically at large 𝑥 , 𝑚 𝜇 is the muon mass, and 𝐺 ( 𝑥 ) is the zero-momentum correlationfunction 𝐺 ( 𝑥 ) = ∫ 𝑑 x (cid:104) 𝐽 𝑒𝑚𝑘 ( 𝑥 ) 𝐽 𝑒𝑚𝑘 ( )(cid:105) (2)of two electromagnetic currents 𝐽 𝑒𝑚𝑘 = 𝑖 (cid:205) 𝑁 𝑓 𝑖 = 𝑞 𝑖 ¯ 𝜓 𝑖 𝛾 𝑘 𝜓 𝑖 , for unexplained notation see Ref. [4]. Herewe consider 𝑁 𝑓 =
3, the 3 lighter quarks of QCD with the first 2 degenerate in mass, so that 𝐺 ( 𝑥 ) = 𝐺 conn 𝑢,𝑑 ( 𝑥 ) + 𝐺 conn 𝑠 ( 𝑥 ) + 𝐺 disc 𝑢,𝑑,𝑠 ( 𝑥 ) . (3)The light-connected Wick contraction 𝐺 conn 𝑢,𝑑 ( 𝑥 ) and the disconnected one 𝐺 disc 𝑢,𝑑,𝑠 ( 𝑥 ) are the mostproblematic contributions with regard to the statistical error. In standard Monte Carlo computations,the relative error of the former at large time distances | 𝑥 | goes as 𝜎 𝐺 connu , d ( 𝑥 )[ 𝐺 connu , d ( 𝑥 )] ∝ 𝑛 𝑒 ( 𝑀 𝜌 − 𝑀 𝜋 ) | 𝑥 | , (4)2 L-HVP
Leonardo Giusti Λ Λ Λ Λ J emk J emk Ω Ω time s p a c e Figure 1:
Domain decomposition of the lattice adopted here. Periodic (anti-periodic) boundary conditionsin the time direction are enforced for gluons (fermions). where 𝑀 𝜌 is the lightest asymptotic state in the iso-triplet vector channel, and 𝑛 is the number ofindependent field configurations. Therefore the computational effort, proportional to 𝑛 , of reachinga given relative statistical error increases exponentially with the distance | 𝑥 | . For the disconnectedcontribution 𝐺 disc 𝑢,𝑑,𝑠 ( 𝑥 ) , the situation is worse since the variance is constant in time and thereforethe coefficient multiplying | 𝑥 | is larger. At present this exponential increase of the relative error isthe barrier which prevents lattice theorists to reach a per-mille statistical precision for the HVP.
3. Multi-level Monte Carlo
Thanks to the conceptual, algorithmic and technical progress over the last few years, it is nowpossible to carry out multi-level Monte Carlo simulations in the presence of fermions [5, 6]. Thefirst step in this approach is the decomposition of the lattice in two overlapping domains Ω and Ω , see e.g. Fig. 1, which share a common region Λ . The latter is chosen so that the minimumdistance between the points belonging to the inner domains Λ and Λ remains finite and positive inthe continuum limit. The next step consists in rewriting the determinant of the Hermitean massiveWilson-Dirac operator 𝑄 = 𝛾 𝐷 asdet 𝑄 = det ( − 𝑤 ) det 𝑄 Λ det 𝑄 − Ω det 𝑄 − Ω , (5)where 𝑄 Λ , 𝑄 Ω , and 𝑄 Ω indicate the very same operator restricted to the domains specified by thesubscript. They are obtained from 𝑄 by imposing Dirichlet boundary conditions on the externalboundaries of each domain. The matrix 𝑤 is built out of 𝑄 Ω , 𝑄 Ω and the hopping terms of theoperator 𝑄 across the boundaries in between the inner domains Λ and Λ and the common region Λ [6]. The denominator in Eq. (5) has already a factorized dependence on the gauge field sincedet 𝑄 Λ , det 𝑄 − Ω and det 𝑄 − Ω depend only on the gauge field in Λ , Ω and Ω respectively. In thelast step, the numerator in Eq. (5) is rewritten asdet ( − 𝑤 ) = det [ − 𝑅 𝑁 + ( − 𝑤 )] 𝐶 (cid:206) 𝑁 / 𝑘 = det (cid:2) ( 𝑢 𝑘 − 𝑤 ) † ( 𝑢 𝑘 − 𝑤 ) (cid:3) , (6)where 𝑢 𝑘 and 𝑢 ∗ 𝑘 are the 𝑁 roots of a polynomial approximant for ( − 𝑤 ) − , the numerator is theremainder, and 𝐶 is an irrelevant constant. The denominator in Eq. (6) can be represented by anintegral over a set of 𝑁 / Λ and Λ inherited from 𝑤 . When the polynomial approximation is properlychosen, the remainder in the numerator of Eq. (6) has mild fluctuations in the gauge field, and it isincluded in the observable in the form of a reweighting factor.3 L-HVP
Leonardo GiustiA simple implementation of these ideas is to divide the lattice as shown in Fig. 1, where Λ and Λ have the shape of thick time-slices while Λ includes the remaining parts of the lattice. Theshort-distance suppression of the quark propagator implies that a thickness of 0 . Λ is good enough, and is not expected to vary significantly with the quarkmass. This is the domain decomposition that we use for the numerical computations presented here.The Monte Carlo simulation is then performed using a two-level scheme. We first generate 𝑛 level-0 gauge field configurations by updating the field over the entire lattice; then, starting from eachlevel-0 configuration, we keep fixed the gauge field in the overlapping region Λ , and generate 𝑛 level-1 configurations by updating the field in Λ and in Λ independently thanks to the factorizationof the action. The resulting gauge fields are then combined to obtain effectively 𝑛 · 𝑛 configurationsat the cost of generating 𝑛 · 𝑛 gauge fields over the entire lattice. Previous experience on two-level integration suggests that, with two independently updated regions, the variance decreasesproportionally to 1 / 𝑛 until the standard deviation of the estimator is comparable with the signal,i.e. until the level-1 integration has solved the signal-to-noise problem. From Eq. (4) we thus inferthat the variance reduction due to level-1 integration is expected to grow exponentially with thetime-distance of the currents in Eq. (2).
4. Lattice computation
In order to assess the efficiency of two-level Monte Carlo integration, we simulated QCD withtwo dynamical flavours supplemented by a valence strange quark on a lattice of size 96 × witha spacing of 𝑎 = .
065 fm, and with a pion mass of 270 MeV. The domains Λ and Λ are theunion of 40 consecutive time-slices, while each thick time-slice forming the overlapping region Λ is made of 8 time-slices. The determinants in the denominator of Eq. (5) are taken into account bystandard pseudofermion representations, while the number of multi-bosons is fixed to 𝑁 =
12. Thevery same action and set of auxiliary fields are used either at level-0 or at level-1. The reweightingfactor is estimated stochastically with 2 random sources, which are enough for its contribution tothe statistical error to be negligible. We generate 𝑛 =
25 level-0 configurations, and for each ofthem, we generate 𝑛 =
10 configurations in Λ and Λ . Further details on the algorithm andits implementation can be found in Refs. [4–6]. To single out the reduction of the variance dueonly to two-level averaging, we carry out a dedicated calculation of correlation functions. Wecompute the light-connected contraction by averaging over 216 local sources put on the time-slice 𝑦 / 𝑎 =
32 of Λ which is at a distance of 8 lattice spacings from its right boundary and, as usual,by summing over the sink space-position. We determine the disconnected contraction by averagingeach single-propagator trace over a large number of Gaussian random sources, namely 768, so tohave a negligible random-noise contribution to the variance [4, 7].The variance of the light-connected contribution as a function of the distance from the source isshown on the left plot of Fig. 2. For better readability only the time-slices belonging to Ω are shown,i.e. those relevant for studying the effect of two-level integration given the source position. Data arenormalized to the variance obtained with the same number of sources on CLS configurations whichwere generated with a conventional one-level HMC. The exponential reduction of the variance withthe distance from the source is manifest in the data, with the maximum gain reached from 2 . 𝑛 =
10. The loss of about a factor between 2 and 3 with respect to the best possiblescaling, namely 𝑛 , either for 𝑛 = L-HVP
Leonardo Giusti . . . . . . { σ / σ H M C } ( x , y ) x − y (fm) G conn u,d n = 1 n = 3 n = 3 ∗ n = 10 10 . . . a σ ( x , y ) x − y (fm) (cid:0) απ (cid:1) { KG conn u,d } × n = 1 n = 3 n = 3 ∗ n = 10 Figure 2:
Left: variance of the light-connected contraction as a function of the difference between the time-coordinates of the currents for 𝑛 = , ,
10. Data are normalized to the analogous ones computed on CLSconfigurations generated by one-level HMC. Dashed lines represent the maximum reduction which can beobtained by two-level integration, namely 1 / 𝑛 , in the absence of correlations between level-1 configurations.Grey bands indicate the thick time-slices where the gauge field is kept fixed during level-1 updates. Right:variance of the light-connected contribution to the integrand in Eq. (1). residual correlation among level-1 configurations. The power of the two-level integration can bebetter appreciated from the right plot of Fig. 2, where we show the variance of the light-connectedcontribution to the integrand in Eq. (1) as a function of the time-distance of the currents. The sharprising of the variance computed by one-level Monte Carlo ( 𝑛 =
1, red squares) is automaticallyflattened out by the two-level multi-boson domain-decomposed HMC ( 𝑛 =
10, blue triangles)without the need for modeling the long-distance behaviour of 𝐺 conn 𝑢,𝑑 ( 𝑥 ) .To further appreciate the effect of the two-level integration, we compute the integral in Eq. (1)as a function of the upper extrema of integration 𝑥 max0 which we allow to vary. For 𝑛 =
1, theintegral reads 446 ( ) and 424 ( ) for 𝑥 max0 = . . 𝑛 =
10 theanalogous values are 467 . ( . ) and 473 . ( . ) . While with the one-level integration the errors onthe contributions to the integral from 0 to 2 . . .
5. Results and discussion
Our best result for the light-connected contribution to the integrand in Eq. (1) is shown onthe left plot of Fig. 3 (red squares). It is obtained by a weighted average of the above discussedcorrelation function computed on 32 point sources per time-slice on 7 time-slices at 𝑦 / 𝑎 = { , , , , , , } and on 216 sources at 𝑦 / 𝑎 =
32. We obtain a good statistical signalup to the maximum distance of 3 fm or so. The strange-connected contraction 𝐺 conn 𝑠 ( 𝑥 ) is muchless noisy, and it is determined by averaging on 16 point sources at 𝑦 / 𝑎 =
32. Its value, shownon the left plot of Fig. 3 (blue circles), is at most one order of magnitude smaller than the lightcontribution, and it has a negligible statistical error with respect to the light one. The best resultfor the disconnected contribution has been computed as discussed in the previous section, and it isshown in the left plot of Fig. 3 as well (green triangles). It reaches a negative peak at about 1 . . L-HVP
Leonardo Giusti . . . x (fm) (cid:0) απ (cid:1) { KG } ( x ) × (fm − ) G conn u,d G conn s G disc u,d,s × − . . . x max0 (fm) a HVP µ ( x max0 ) × G conn u,d G conn s G disc u,d,s × Figure 3:
Left: best results for the contribution to the integrand in Eq. (1) from the light-connected (redsquares), strange-connected (blue circles) and disconnected (green triangles) contractions as a function ofthe time coordinate. Right: best results for the contributions to 𝑎 HVP 𝜇 from light-connected (red squares),strange-connected (blue circles), and disconnected (green triangles) contractions as a function of 𝑥 max0 . strange-connected (blue circles), and disconnected (green triangles) contributions to 𝑎 HVP 𝜇 · as a function of the upper extrema of integration 𝑥 max0 in Eq. (1). The light-connected part startsto flatten out at 𝑥 max0 ∼ . 𝑥 max0 = . . ( . ) . The value of the strange-connected contribution is 52 . ( ) at 𝑥 max0 = . 𝑥 max0 ∼ . − . ( ) . For 𝑥 max0 = . . 𝑥 max0 = . 𝑥 max0 = . 𝑎 HVP 𝜇 = . ( . ) · − .In this proof of concept study we have achieved a 1% statistical precision with just 𝑛 · 𝑛 = 𝑎 HVP 𝜇 is reachable with multi-level integration by increasing 𝑛 and 𝑛 by a factor ofabout 4–6 and 2–4 respectively. When the up and the down quarks becomes lighter, the gain due tothe multi-level integration is expected to increase exponentially in the quark mass, hence improvingeven more dramatically the scaling of the simulation cost with respect to a standard one-level MonteCarlo. The change of computational paradigm presented here thus removes the main barrier formaking affordable, on computers available today, the goal of a per-mille precision on 𝑎 HVP 𝜇 . References [1] G. W. Bennett, et al., Phys. Rev. D73 (2006) 072003, arXiv:hep-ex/0602035 .[2] J. Grange, et al., arXiv:1501.06858 .[3] T. Aoyama, et al., arXiv:2006.04822 .[4] M. Dalla Brida, L. Giusti, T. Harris and M. Pepe, arXiv:2007.02973 .[5] M. Cè, L. Giusti, S. Schaefer, Phys. Rev. D93 (9) (2016) 094507, arXiv:1601.04587 .[6] M. Cè, L. Giusti, S. Schaefer, Phys. Rev. D95 (3) (2017) 034503, arXiv:1609.02419 .[7] L. Giusti, et al., Eur. Phys. J. C79 (7) (2019) 586, arXiv:1903.10447arXiv:1903.10447