The charm-quark contribution to light-by-light scattering in the muon $(g-2)$ from lattice QCD

We compute the hadronic light-by-light scattering contribution to the muon $g-2$ from the charm quark using lattice QCD. The calculation is performed on ensembles generated with dynamical $(u,d,s)$ quarks at the SU(3)$_{\rm f}$ symmetric point with degenerate pion and kaon masses of around 415 MeV. It includes the connected charm contribution, as well as the leading disconnected Wick contraction, involving the correlation between a charm and a light-quark loop. Cutoff effects turn out to be sizeable, which leads us to use lighter-than-physical charm masses, to employ a broad range of lattice spacings reaching down to 0.039 fm and to perform a combined charm-mass and continuum extrapolation. We use the $\eta_c$ meson to define the physical charm-mass point and obtain a final value of $a_\mu^{\rm HLbL,c} = (2.8\pm 0.5) \times 10^{-11}$, whose uncertainty is dominated by the systematics of the extrapolation. Our result is consistent with the estimate based on a simple charm-quark loop, whilst being free of any perturbative scheme dependence on the charm mass. The mixed charm-light disconnected contraction contributes a small negative amount to the final value.


I. INTRODUCTION
The anomalous magnetic moment of the muon, a µ ≡ (g − 2) µ /2, is one of the most precisely measured quantities in fundamental physics. Currently, the experimental world average [1,2] and the theoretical evaluation of the 2020 White Paper (WP) [3] based on the Standard Model (SM) of particle physics are in tension at the 4.2σ level. The theory uncertainties are entirely dominated by the hadronic contributions. Surprisingly, a lattice-QCD based calculation [4] of the leading hadronic contribution finds a larger value than the dispersion-theory based estimate of the WP, which would bring the overall theory prediction into far better agreement with the experimental value of a µ . Thus it will be vital to resolve the tension between the different determinations of the leading hadronic contribution in order to strengthen the unique test of the SM offered by the anomalous magnetic moment of the muon.
A subleading hadronic contribution to a µ , the hadronic light-by-light (HLbL) contribution, also contributes sizeably to the error budget of the SM prediction. The HLbL contribution is significantly more complex to evaluate than the leading hadronic contribution; however, because it is suppressed by an additional power of the fine-structure constant α, it only needs to be determined at the ten percent level. The HLbL contribution, too, has been evaluated using either dispersive methods [3] or lattice QCD [5,6]. In this case, good agreement is found among the three evaluations within the quoted uncertainties.
One missing ingredient in the otherwise complete HLbL calculation of [6] is the contribution of the charm quark. The present paper addresses this missing contribution. Since the charm quark is much heavier than the muon, on general grounds [7,8] one expects this contribution to be in a regime where it is roughly proportional to m 2 µ /m 2 c . In phenomenological estimates, it has been evaluated using the prediction based on Quantum Electrodynamics (QED), amended for the appropriate charge and colour factors. We quote the value and uncertainty from the 2020 White Paper [3], a HLbL,c µ (WP) = (3 ± 1) × 10 −11 . (1) The main goal of this paper is thus to test the prediction (1) using lattice QCD, in case a qualitative effect might have been missed. Certainly, this contribution is small compared to the overall uncertainty 43×10 −11 of the WP prediction for a µ , however the other uncertainties are also expected to shrink, especially if the issues in the leading hadronic contribution can be resolved. Our second motivation for addressing the charm HLbL contribution from first principles is to answer the qualitative question whether approximating this contribution via a simple quark loop is adequate. In lattice QCD, the calculation involves computing charm propagators on an ensemble of non-perturbative background SU(3) gauge fields. If the simple quark-loop picture is approximately correct, the details of this gauge field should not matter much, and the charm propagators can be replaced by free Dirac propagators. In this case, the sensitivity to the sea quarks enters (at the earliest) at quadratic order in α s (m c ), the strong coupling constant at the scale of the charm mass. It is largely for this reason that we will focus on the SU(3) f -symmetric mass point with m π = m K 415 MeV, enabling us to reach sufficiently fine lattices at a moderate computational cost.
A further aspect of the quark-loop picture is that the various disconnected diagrams entering the HLbL amplitude are expected to be small. In contrast, if the η c pole exchange or D meson loops played a sizeable role in the charm-quark contribution, the leading disconnected charm contribution, consisting of a charm loop and a light-quark loop, each attached to two electromagnetic currents, would be sizeable (in analogy to the analyses in [9] and appendix A of Ref. [6] for the three-flavour case). We recall that for the light quarks, individual mesons, especially the pseudoscalars π 0 , η, η , contribute substantially to a HLbL µ , even at the aforementioned SU(3) f -symmetric point [10]. In lattice QCD, we can quantitatively test the relevance of the disconnected contributions. This paper is organized as follows. We describe our lattice setup, the tuning of the charm quark mass and our specific representation of a HLbL µ in section II. Section III provides some basic theory expectations concerning the connected and leading disconnected contributions involving a charm quark. Section IV presents our lattice results on the connected contribution for a sequence of increasing charm-quark masses, and section V contains our results at the target charm mass for the leading topology of disconnected diagrams. We provide our final result and conclude in section VI. Appendix A describes a test of our methods at a heavy quark mass in lattice QED, while appendix B contains tables of results for the connected charm contribution on individual ensembles and appendix C a representative set of fit results.  (3) f -symmetric ensembles used in this work. Each ensemble is parametrized by the gauge coupling parameter β ≡ 6/g 2 0 , the (u, d, s)-quark hopping parameter κ, the lattice size, and the temporal boundary condition. The lattice spacings a were determined in Ref. [12], apart from A653 and J500, where the lattice spacings were estimated from the ratio of the Wilson flow parameter t 0 ; the errors on the lattice spacing for these two ensembles (in bold) are simply estimated by scaling of the total error of the neighboring lattice spacings. Their pion masses (marked with asterisk) have been measured independently for this work.

II. LATTICE SETUP
We have performed lattice-QCD calculations on gauge ensembles provided by the Coordinated Lattice Simulations (CLS) initiative [11], which have been generated using three flavours of non-perturbatively O(a)-improved Wilson-clover fermions and with the tree-levelimproved Lüscher-Weisz gauge action. As in Ref. [10], where we computed the (u, d, s) quark contribution, we consider only ensembles realizing exact SU(3) f -symmetry. On these ensembles, the mass of the octet of light pseudoscalar mesons is approximately 415 MeV. The parameters of these ensembles, which correspond to six different values of the lattice spacing, are summarized in Table I. A. Calibrating the charm mass and current The connected contribution to a HLbL µ and the two-point correlation function ofcγ 5 c were computed on all ensembles of Table I for several (5 or 7) values of the charm-quark bare subtracted mass am c = (κ −1 c − κ −1 crit )/2, with values of κ c chosen to interpolate between the physical strange and charm hopping parameters. A determination of the latter is available from Ref. [13], obtained by tuning the D s meson mass to its physical value. For the (dominant) connected contribution however, we choose the physical charm-mass point as the one defined by the physical value of the η c meson mass. When we determine the η c mass, we do not include the disconnected diagram in the two-point function of the charm pseudoscalar density. This procedure corresponds to using the operatorc γ 5 c (am c = am c ), where the degenerate quark flavours c and c are both treated at the partially-quenched level. It should be noted that the tuning of Ref. [13] by the D s yields a heavier-than-physical η c meson at our SU(3) f -symmetric point. This comes from the quark masses at the latter point being lighter than the physical strange quark, and a D s -tuning de facto absorbs this effect into the charm-quark mass [14].
The reason for using lighter-than-physical charm quark masses is that we expect discretisation effects to become more and more significant when the charm mass increases. For a rough estimate of the typical size of discretisation effects, [14] found that the effective speed of light (as defined by the dispersion relation of a meson) for physical-mass charm quarks at worst deviates from unity by 20% in our setup.
The finite renormalisation factor Z c V (g 0 , am c ) for the local charm currentcγ µ c was determined by requiring the corresponding charge of the ground-state meson created when c γ 5 c acts on the vacuum to be unity. The meson correlators were computed using Z 2 × Z 2 stochastic wall sources [15,16]. The quark-mass dependence of Z c V (g 0 , am c ) is quite strong, especially at coarse lattice spacings. Since this factor enters to the third power into our final result, we determine it directly for every one of the bare quark-mass values. This is the same procedure that was implemented for the charm renormalisation in [13]. 1 B. Computing the charm contribution to a HLbL µ We apply the formalism described and used in [6,10] and therefore only recall the main aspects. The starting point of our calculation is the master formula 2 a HLbL Here e 2 /(4π) = α QED is the fine-structure constant and m µ the muon mass. The QED kernelL has been computed in the continuum [17] and represents the contributions of the photon and muon propagators and vertices in the diagrams of Fig. 1. There is a lot of freedom to alter the kernel without changing a HLbL µ in the continuum and in infinite volume. Specifically, we use the kernelL Λ defined in [10] with Λ = 0.40 throughout. The tensor i Π is a Euclidean hadronic four-point function with one of its vertices weighted linearly in one of its coordinates, ( The field j µ (x) appearing above is the hadronic component of the electromagnetic current, Here we focus on the contributions involving the charm current,cγ µ c. The QCD fourpoint function receives contributions from five classes of Wick contractions. First, we will focus on the fully-connected charm contribution, which involves four charm currents; for this contribution, we apply Eq. (7) of Ref. [6] with the flavour index set to charm, j := c. Second, we will consider the disconnected contributions involving two quark loops, each of which containing two vector vertices, with either one or both loops consisting of charm propagators. Here we apply Eq. (11) of Ref. [6] with the flavour indices i, j running over {u, d, s, c} under the constraint that at least one of them take the value c. The connected and (leading) disconnected contributions are illustrated in Fig. 1. 1 In that paper, the Z c V values quoted in its appendix A are erroneously described as stemming from the charm number of the D s meson. 2 See however the text below Eq. (4) for references to the precise formulae used in the present calculation.

III. THEORY EXPECTATIONS
The simplest prediction for the light-by-light contribution of a heavy 'charm' quark to 10 11 a µ relies on the analytic QED result originally applied to the τ lepton loop [18,19]. Taking into account the colour factor N c = 3 and the charge factor (2/3) 4 , it is given by the function withm Q is the heavy-quark mass in GeV. Already by m Q = 0.75 GeV, the O(m −4 Q ) terms only represent a reduction of the leading term by five percent. These terms certainly represent a small correction for m Q around the physical charm mass. Here we have dropped known higher-order terms in 1/m Q . We will take the function h(m Q ) as a baseline for comparison with our lattice results for the fully connected charm contribution.
For the (2+2) disconnected contribution involving one charm and one light-quark loop, it is less straightforward to make a 'baseline' prediction. The scalar-QED prediction for the contribution to a HLbL µ of the D ± meson loop is −0.33 × 10 −11 [19], to be roughly doubled in order to include the D s loop. Taking into account the charge factor of 2·3·Q 2 c (Q 2 u +Q 2 d +Q 2 s ) = 144 81 relevant for the charm-light (2+2) contribution (see [6], appendix A 3 ), one arrives at the prediction of a 2+2:lc µ = −0.58 × 10 −11 when treating the D + , D 0 , D s meson loops within scalar QED. 4 The absolute value of this prediction is surely an overestimate, given that electromagnetic form factors of the D mesons should suppress this prediction substantially: in the case of the pion loop, the suppression factor is almost three, and for the kaon almost ten [3,19]. All in all, these considerations finally lead us to expect an order of magnitude of (−0. In addition to the short-distance effect estimated above, the charm-light disconnected diagrams also involve a longer-distance contribution, whose size it is useful to estimate by theory arguments, given the difficulty of measuring the correlation function in the infrared. The intuitive idea is that the heavy-quark loop shrinks almost to a point in coordinatespace 5 , acting effectively like a local gauge-invariant gluonic operator from the point of view 3 The analysis of this reference can be applied here due to the u, d, s quarks being degenerate and the charm quark being quenched. 4 Note that the D 0 loop contributes to a 2+2:lc µ , even though at the scalar-QED level it does not contribute to a HLbL µ , due to it cancelling between different topologies. 5 This picture holds when the vertices of the light-quark loop are at a distance much greater than (2m c ) −1 from the charm loop. of the 'low-energy effective theory', which is QCD with (u, d, s) quarks. This picture can be formalized by writing an effective Lagrangian for the effective coupling induced between two photons and gluonic fields, much as in the classic work of Euler and Heisenberg [20]. This effective Lagrangian L (c) 2γ2g has been calculated long ago [21]; each term of the Lagrangian contains two photonic and two gluonic field strength tensors. From here, one infers the operator equation A µ being the photon field, which shows that the charm loop acts at low energies like a set of gluonic operators such as α s G a µν G a µν or α s G a µνG a µν . The main observation is that, on dimensional grounds, the effective Lagrangian is overall multiplied by a 1/m 4 c factor, indicating a strong suppression.
The argument above shows that a light flavour-singlet meson such as the scalar f 0 or the pseudoscalar η can propagate between the charm loop and the light-quark loop, albeit with a very suppressed coupling to the charm loop. To get an estimate of this contribution, which is long-range in comparison to the length-scale (2m c ) −1 , we use Eq. (6) to find out roughly how much the charm part of the electromagnetic current by itself contributes to the η transition form factor (TFF). Note that this contribution is independent of the photon virtualities, as long as these are small. Using the estimate 0|α s GG|η ≈ 0.5 GeV 3 based on Ref. [22], while the TFF normalisation amounts to |F η γγ (0, 0)| 0.34 GeV −1 (see for instance [23]), we obtain a contribution of about 8 × 10 −4 GeV −1 to F η γγ from the charm current. Since the η exchange contributes about 14.5 × 10 −11 to a HLbL µ [3], proportionally to its TFF at each end of η propagator, we arrive at the order-of-magnitude estimate of 0.01 × 10 −11 for the contribution to a HLbL µ of the η in the (2+2) charm-light diagrams. Even with a potential logarithmic enhancement [24], this is much smaller than our final uncertainty and cannot presently be resolved in our lattice calculations.
In addition to the Wick-contraction topologies considered above, the (3+1) topology with the single-current loop consisting of a charm propagator deserves some attention, since this contribution is neither SU(3) f nor 1/N c suppressed relative to the (2+2) topology, N c being the number of colours. In perturbation theory, the (3+1) contribution starts at O(α 3 s ) rather than at O(α 2 s ), while involving the same minimal number of charm propagators. Furthermore, the quark-charge and multiplicity factors numerically suppresses this contribution by a relative factor of three 6 since it is weighted by 4 · (Q 3 u + Q 3 d + Q 3 s )Q c = 48/81, while the charm-light (2+2) diagrams are weighted by 144/81, as noted above. A factor of three suppression relative to the (2+2) charm-light contribution is thus expected.

IV. LATTICE RESULTS FOR THE CONNECTED CONTRIBUTION
As a way of validating our computational methods, our analysis has been guided by a lepton-loop calculation, much like in Ref. [10]: in Appendix A we investigate the applicability of our QED-kernel implementation at particularly heavy scales by comparing the leptonloop contribution to a LbL µ to the known analytical expression [18]. While the agreement 6 Within the scalar QED framework, the two topologies by themselves contain equal and opposite contributions from the D meson loops, since they cancel in cγ µ cūγ ν uūγ ρ u f =u,d,s,cq f γ λ q f , given that the charge of D mesons is zero under the total quark number current f =u,d,s,cq f γ λ q f . is acceptable at fairly heavy lepton mass, the study suggests that cut-off effects will be significant and working at unphysically-light charm mass might allow for a better handle on these effects. The physical charm mass will therefore be approached via a simultaneous extrapolation in the quark mass and in the lattice spacing.

A. Results at individual quark masses
For the connected part of a HLbL,c µ , we have performed computations with the vector current connected to the external on-shell photon (the z-vertex in Eq. (3)) being either symmetrisedconserved 7 or local, while the rest of the currents are kept local. For each ensemble, we have tuned κ c to get five to seven different η c -masses, ranging from around 1.3 to 2.6 GeV. In order to better control rotational-symmetry breaking effects (and keep the higher-order lattice artifact coefficients the same) we will only use f (|y|) along the lattice direction (1,1,1,1) for all ensembles.   2 shows an example of our data for the A653 ensemble. The integrand is steeplypeaked at short distances and becomes more so at heavier quark masses (smaller κ c ). As can be seen from the partially-integrated results, even the lightest charm-quark-mass lattice data used here completely saturates the integral and therefore there is no need to perform any tail-extension procedure, and just the lattice (trapezoid-rule) integral suffices. We also note that the overall integrand and integral becomes substantially smaller as κ c decreases, representing the fact that this integral must vanish in the limit κ c → 0. There is a strong negative tail in the integrand causing a fairly significant cancellation for the overall integral, which becomes smaller as the charm-mass decreases. At very low κ c on coarse lattices it is unlikely that we will be able to properly resolve the peak of the integrand and end up with a lower estimate due to the negative tail cancelling against the peak contribution more than it should. As we move to finer lattices and the resolution at low |y| improves, we resolve the peak structure much better, as illustrated in Fig. 3. The results are given in Tabs. III -VIII of Appendix B and summarized in Fig. 4. Expectations are that a µ scales with m 2 µ /m 2 heavy [7,8], so it is instructive to focus on the dependence of a HLbL,c µ on 1/M 2 ηc . The data show a clear monotonic decrease as 1/M 2 ηc is decreased toward its physical value, starting (for the lightest charm quarks) at or above the WP prediction and ending (for the heaviest charm quarks) at or below the WP value. At similar η c mass, the data have a large spread between the coarsest and finest ensemble, indicating strong discretisation effects.
At this point, it is useful to compare the two choices of discretisations for the currents: the spread is larger in the local-local data than in the local-conserved data. Furthermore, the curvature in 1/M 2 ηc has a stronger dependence on the lattice spacing in the local-local data. In addition, the fact that the coarse local-local data at large M ηc become negative makes it more difficult to describe the data using a fit ansatz. For these reasons, we decide to base our determination of a HLbL,c µ solely on the analysis of the local-conserved data.

C. Extrapolation to the continuum and to the physical charm mass
Due to the heaviness of the valence charm quark, the intermediate states that could potentially contribute to the correlation function in question should be much suppressed at large distances; see the discussion in section III. Indeed, this can be seen by the saturation of the tail of the lattice integrand (Figs. 2 and 3). For this reason, in the approach to the physical point, we assume that the finite-size effects are minor and only extrapolate in the η c -meson mass and lattice spacing a. The statistical error on each individual data point is at the percent-level, which is comparable to the quoted error on the lattice spacings given in Tab. I; therefore, it is crucial to include the error on the lattice spacing while performing an extrapolation to the physical point.
To this end, a global fit is performed based on a Bayesian approach [26], where we promote each lattice spacing to a fit-parameter and associate to it a Gaussian prior with the central value and the width taken to be the quoted central value of the lattice spacingā and its error ∆a respectively. Although the parameter space is small, constructing a fit-ansatz with a χ 2 /dof on the order of unity is in fact not an easy task. After various attempts, we have identified two classes of ansätze which are able to describe the data with reasonably good χ 2 /dof. The most restrictive constraint that we deem important to fulfill is the m 2 µ /m 2 heavy scaling of a HLbL,c µ in the presence of a heavy scale [7,8]. It is natural to first consider the η c -meson mass for such a scale. A challenging part of the construction of fit-ansätze is to handle the apparent non-linear behavior in 1/M 2 ηc of the data (see Fig. 4), which gradually gets milder as we go down in the lattice spacing. This motivates our first class of ansätze, the P-class, which consist of linear combinations of a leading term in 1/M 2 ηc and terms in a n f (M ηc ) with n ∈ N * and f an elementary function, treating the non-linearity in 1/M 2 ηc as a lattice artifact.
Another way to account for the m 2 µ /m 2 heavy scaling is to use the charm quark mass as a heavy scale. A rough estimate in the non-relativistic limit is that the η c -meson mass should be equal to twice the charm quark mass, up to small relative corrections. Based on this observation, we define the R-class of fit ansätze, consisting of rational functions: where P and Q are polynomials in both a and M ηc and C is a constant. In principle, C can also have non-trivial dependence on a and on M ηc ; however, introducing additional parameters to describe this dependence turns out to be unnecessary, as the non-linearity of the data can already be well captured with the form in Eq. (7). With the aforementioned two fit-ansatz classes, it remains nevertheless difficult to get reasonable χ 2 /dof with the whole available dataset. In fact, this is not very surprising, as the resolution of the peak of the integrand becomes worse as κ c and a become small (see Fig. 2). Therefore, it is necessary to allow for various cuts to the data. At the same time, as we would like to reach as heavy as possible M ηc masses in order to have a better control over the extrapolation to its physical value, it is preferable to discard as few data points as possible. A lattice study in the pure QED case presented in App. A shows that our setup should be valid up to a charm-quark mass of at least 20/3 times that of the muon, with well-controlled cut-off effects. Based on the latter and with a simple linear relation between the physical η c mass and the MS-mass of the charm quark [27], we demand an admissible fit to be able to cover the data points in the range of 1/M 2 ηc > 0.20 GeV −2 . Our fitting strategy goes as follows: We build fit ansätze from either the P-or the R-class as explained earlier. To avoid overfitting, the number of fit parameters is limited to five. Apart from terms in a n M m , we have also tried logarithmic terms in a or M in order to allow for different types of curvature. The four datasets we consider are (D2, D3, and D4 are defined as D1 with extra omissions): Given that our final uncertainty estimate is dominated by systematics due to the choice of fit ansatz and that attempts at correlated fits yielded a poor fit quality, we choose to neglect correlations between different M ηc on the same ensemble. Although this harms the statistical interpretation of χ 2 and p-values computed in the standard way, we nevertheless use these to judge relative fit quality. Our criterion for an admissible fit is one with a p-value between 0.05 and 0.95, for which the extrapolated a µ and the p-value are stable under variation of the dataset choice. We have tested various fit ansätze from both the P-and R-classes and found that the following five-parameter fits are able to describe our data with the quality requirements A further important feature of these fits is that they qualitatively follow the trend of the data even in the region 1/M 2 ηc < 0.2 GeV −2 . As a general feature, the P-class ansätze tend to lead to larger results for a µ (0, M Phys ηc ) as compared to the R-class. As there is no exclusive theoretical argument for the finite-latticespacing behaviour of these functions, and our data seem not to be able to unambiguously rule out any of these classes, our decision is to include the fit results with good χ 2 /dof from both of them (see Tab. IX). More specifically, our final result is the average of our lowest (Fit 1, D3 : 2.64(4)) and our largest (Fit 5, D2 : 3.47(3)) values and we assign a generous systematic error estimate by quoting half the difference of the two, which brings us to our estimate for the connected contribution, a HLbL,c µ ,conn. = 3.1(4) × 10 −11 .
Our error on this quantity is entirely dominated by the systematic error from our modeling of its dependence on a and M ηc .

D. Comparison to the QED-based prediction
To close the study of the connected contribution, we compare our result for a HLbL,c µ ,conn.
to the charm-quark loop evaluated analytically within QED as given in Eq. (5). To make contact with that expression, we need to specify the relationship between the η c mass and the charm-quark mass. As explained while discussing the R-class of fit-ansätze, we assume the η c mass to be twice the charm quark mass plus an almost charm-quark-mass independent mass-shift within a given window of M ηc . We estimate the mass-shift using the MS charmquark mass and the physical M ηc and assign a five-percent uncertainty to this quantity. The prediction from this prescription is displayed in Fig. 6 together with our fit results. The difference between truncating at O(1/m 2 c ) and at O(1/m 4 c ) is tiny compared to the uncertainty that we assign to the mass-shift inferred from our prescription. It is worth noting that, even though the QED-based prediction gives a result that falls in the bulk of our estimate Eq. (9) at the physical charm mass, the milder curvatures in 1/M 2 ηc of the representative fit results suggest that non-perturbative effects are still significant at intermediate masses. : comparison between our continuum-extrapolated results using Fit 1 (lower, in magenta) and Fit 5 (higher, in red) based on dataset D1 and the QED-based prediction (in blue). The band on the latter indicates a ±5% change in the mass-shift relating the charm-quark mass to the η c mass (see text). The dotted vertical line indicates the upper bound for the η c mass included in the fits. Our estimate for the connected contribution, Eq. (9), is marked with 'Conn.' (a horizontal offset is applied for visibility).

V. THE DISCONNECTED CONTRIBUTION
The disconnected parts of the charm contribution are expected to be very small. From the outset, we neglect the (3+1),(2+1+1), and (1+1+1+1) Wick-contraction topologies, based partly on them being consistent with zero for the light quark contribution, as found in [6], and partly on the arguments laid out in section III. This leaves us with the (2+2) topology, which is a sizeable contribution in the light-quark a HLbL µ result. This contribution can be broken into the mixed charm-light and the charm-charm contributions, with the former (by analogy with the strange sector) expected to be the major contribution.
As the disconnected contribution is still an expensive calculation, we have limited ourselves to a single charm-quark mass determined by κ c from the D s -tuning of Ref. [13]. This tuning is suboptimal for our present purposes, an aspect we return to below. We will also use the Z c V values from [13], except for ensemble A653, where we computed the renormalisation factor ourselves. We employ exclusively local vector currents and restrict ourselves to the ensembles N300, N202, B450, and A653, reusing data for the light-quark loop from Ref. [10].
A plot of the partially-integrated charm-light and charm-charm disconnected contributions for ensemble N202 is shown in Fig. 7. It is clear that both of these contributions are noisy, small, negative, and very short-distance. Again this means we can use the lattice integral directly for our final result, and we use a simple constant fit to the partially-integrated result for our final determination. Based on the numerical evidence from Figs. 2-3 that the connected contribution becomes very short-ranged as the charm mass is increased towards its physical value, as well as the theoretical arguments of section III, we start this fit between 0.4 and 0.5 fm. Tab. II shows our results for this procedure and we see that A653 is an extreme outlier in the charm-light contribution. The other, finer, ensembles yield values much smaller and consistent with one another. We decide to omit this coarse ensemble entirely and fit the remaining charm-light data to a straight line in the variable a 2 . This leads to the result a HLbL,c,(2+2) µ = −0.28(21) × 10 −11 . We now come back to the issue of the tuning of the charm quark mass. The CLS ensembles we are using are designed to have the trace of the quark mass matrix equal to its physical value, to a rather good approximation [12]. We remind the reader that for the connected contribution, we chose the charm-quark mass such that the physical η c mass is reproduced. With this choice, the dependence of charm correlators on the SU(3) f breaking parameter [m s − (m u + m d )/2] is expected to be small, being a pure sea quark effect. As a consequence, the extrapolation to physical (u, d, s) quark masses is expected to be very mild. This is not the case if we tune the massM D of the triplet of D mesons at our SU (3)  GeV of the physical D meson masses, then we again avoid a valence-quark effect in the approach to physical quark masses. It is also interesting to ask, how different a tuning this represents as compared to the tuning via the η c mesons mass. We have found that the η c meson mass, extrapolated to the charm mass whereM D = (M D ) phys , amounts to 2.97(4) GeV, which is consistent with its physical value. This is an indication that sea quark effects are indeed small in the charm sector.
These observations lead us to apply a small correction to the charm-light disconnected contribution, to bring it to the point whereM We neglect the charm-charm contribution as its contribution is far smaller than our final  error for the charm-light.

VI. DISCUSSION OF THE RESULTS AND CONCLUSION
We have determined the charm-quark contribution to hadronic light-by-light scattering in the anomalous magnetic moment of the muon. We find that the lattice determination of this quantity is challenging, specifically in the modeling of the connected contribution's discretisation effects: the associated systematic error dominates our final error budget. As expected from the charm-loop picture, the connected contribution turns out to be the most significant overall. We find the charm-light disconnected contribution to be negative and much smaller in magnitude than the fully-connected contribution, amounting to a 10% correction with a large uncertainty. The charm-charm disconnected contribution is entirely negligible and we expect all higher-order contributions to be equally insignificant.
Before quoting our final result for the charm contribution to a HLbL µ , we address the question of its dependence on the (u, d, s) quark masses. The fact that several aspects of our lattice results can be understood via the the charm-quark loop picture is an indication that this dependence must be modest, and we may attempt to estimate its order of magnitude via the ambiguity induced by the choice of the charm-quark tuning condition away from the physical (u, d, s) quark-mass point. We saw in section V that tuning the average D + , D 0 and D s mass to its physical value was equivalent, within our uncertainties, to tuning the η c mass to its physical value. Still, we estimate that the connected contribution would potentially be modified by 2% had we chosen the alternative tuning. Another estimate can be based on the idea that the charm contribution is proportional to the sum of the inverses of the charged D-meson squared masses. That sum differs again by about two percent from the square inverse of the average D + , D 0 and D s mass. This argument suggests an absolute systematic error of 0.06 × 10 −11 , which we conservatively inflate to 0.12 × 10 −11 and add in quadrature to the other uncertainties below.
Our full result from adding Eqs. (9) and (10) This result is completely consistent with the 2020 White Paper estimate of (3±1)×10 −11 [3], and has half its uncertainty. Combining Eq. (11) with our previous result from the light and strange contributions of a HLbL,ls µ = (106.8 ± 15.9) × 10 −11 [6] obtained with dynamical (u, d, s) quarks yields a fully non-perturbative determination of a HLbL µ , including all relevant contributions. A last effect not yet accounted for is the charm sea-quark effect on the light-quark contributions, as for instance the D + meson loop can contribute to the connected four-point function of the down quark in a calculation with dynamical (u, d, s, c) quarks. Within a scalar-QED treatment of the D meson, we have however estimated this effect to be below 0.1 × 10 −11 . Therefore we neglect the charm sea quark effects and arrive at a HLbL µ = (109.6 ± 15.9) × 10 −11 .
This concludes our first-generation calculation of hadronic light-by-light scattering in the muon (g − 2).

ACKNOWLEDGMENTS
We thank Andreas Nyffeler for a fruitful collaboration on computing the QED kernel used here. This work is supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme through grant agreement 771971-SIMDAMA, and through the Cluster of Excellence Precision Physics, Fundamental Interactions, and Structure of Matter (PRISMA+ EXC 2118/1) within the German Excellence Strategy (Project ID 39083149). The project leading to this publication has also received funding from the Excellence Initiative of Aix-Marseille University -A*MIDEX, a French "Investissements d'Avenir" programme, AMX-18-ACE-005. JRG acknowledges support from the Simons Foundation through the Simons Bridge for Postdoctoral Fellowships scheme. Calculations for this project were performed on the HPC clusters "Clover" and "HIMster II" at the Helmholtz-Institut Mainz and "Mogon II" at JGU Mainz. Our programs use the deflated SAP+GCR solver from the openQCD package [28], as well as the QDP++ library [29]. We are grateful to our colleagues in the CLS initiative for sharing ensembles.
Appendix A: Methodology test for a heavy lepton With our implementation of the QED coordinate-space kernel, we have been able to reproduce various known light-by-light contributions in the continuum [10,30,31] at the one-percent level. The tests performed so far concern physics involving particles with masses on the same order as the muon mass. As our implementation of the QED-kernel relies on interpolating weight functions that are precomputed on a grid [17], it is important for the goal of this paper to test how robust this implementation is for computing contributions from more massive particles.
As an example of a calculation performed entirely in the continuum, we quote the result we obtain with the kernelL (Λ=0.40) and 'method 2' for the lepton-loop contribution with m /m µ = 4, namely a HLbL µ = (42.1 ± 0.5) × 10 −11 ; the exact result is (43.175...) × 10 −11 . While this precision is sufficient for our present purposes, it is clear that, using continuum propagators for the lepton loop, the quality and stability of the coordinate-space results degrade when the lepton mass increases. 8 In order to validate our computational setup, we turn to a test that is much closer to the procedure we used for the charm-quark contribution in lattice QCD. We have computed the lepton-loop contribution to a HLbL µ using lattice fermion propagators at a mass-scale 8 In our estimate, we regulate the numerics by setting the integrand to zero when two vertices come within a distance of 10 −3 m −1 µ . A more sophisticated procedure should be used for higher lepton masses. relevant for this project, choosing specifically m /m µ = 20/3. We proceed by repeating the calculation on increasingly fine lattices and perform a continuum extrapolation using a quadratic polynomial in am µ . Here, two discretisations of the vector current at the external vertex z were used, and the resulting contributions to a LbL µ were extrapolated simultaneously to the continuum, enforcing a common continuum value. The deviation of the continuumextrapolated result from the known exact result of 16.395..×10 −11 depends somewhat on the choice of the extrapolation range, but is in all cases within 2.5%. This successfully passed test gives us confidence that the setup used for the lattice-QCD calculation presented in the main text is robust for fermion masses up to at least 700 MeV.  In this appendix we collect details of the fit results to the local-conserved connected data obtained with the fit ansätze of Eq. (8). The values obtained for a HLbL,c µ ,conn. , as well as the corresponding χ 2 /dof and p values are given in Table IX.