The charm-quark contribution to light-by-light scattering in the muon \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(g-2)$$\end{document}(g-2) from lattice QCD

We compute the hadronic light-by-light scattering contribution to the muon \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g-2$$\end{document}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)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_\mathrm{f}$$\end{document}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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _c$$\end{document}ηc meson to define the physical charm-mass point and obtain a final value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_\mu ^\mathrm{HLbL,c}= (2.8\pm 0.5) \times 10^{-11}$$\end{document}aμHLbL,c=(2.8±0.5)×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.


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 a e-mail: meyerh@uni-mainz.de (corresponding author) 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. In the past decade, the HLbL contribution, too, has been evaluated using either dispersive methods [5], assisted by short-distance constraints and hadron structure input [6][7][8][9][10][11][12][13][14][15][16][17][18], or lattice QCD. In this case, good agreement is found between the dispersive [3] and the two lattice evaluations [19,20] within the quoted uncertainties.
One missing ingredient in the otherwise complete HLbL calculation of [20] 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 [21][22][23] 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 [24]. We quote the value and uncertainty from the 2020 White Paper (WP) [3], a HLbL,c μ (WP) = (3 ± 1) × 10 −11 . (1) While the central value comes from evaluating the QEDlike quark-loop, the uncertainty has been estimated conservatively based on the size of the η c pole contribution [12].
Since the WP appeared, the leading radiative correction to a massless quark loop has also been computed [25]. 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 [26] and Appendix A of Ref. [20] 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 [27]. 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 Sect. 2. Section 3 provides some basic theory expectations concerning the connected and leading disconnected contributions involving a charm quark. Section 4 presents our lattice results on the connected contribution for a sequence of increasing charm-quark masses, and Sect. 5 contains our results at the target charm mass for the leading topology of disconnected diagrams. We provide our final result and conclude in Sect. 6. 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.

Lattice setup
We have performed lattice-QCD calculations on gauge ensembles provided by the Coordinated Lattice Simulations (CLS) initiative [28], which have been generated using three flavours of non-perturbatively O(a)-improved Wilson-clover fermions and with the tree-level-improved Lüscher-Weisz gauge action. As in Ref. [27], 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 1.

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 1 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. [30], 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 operator c γ 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. [30] 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 [31].
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, [31] 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 whenc γ 5 c acts on the vacuum to be unity. The meson corre- Table 1 The SU(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. [29], apart from A653 and J500, where the lattice spacings were esti-mated 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 lators were computed using Z 2 × Z 2 stochastic wall sources [32,33]. 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 [30]. 1 2.2 Computing the charm contribution to a HLbL μ We apply the formalism described and used in [20,27] and therefore only recall the main aspects. The starting point of our calculation is the master formula 2 Here e 2 /(4π) = α QED is the fine-structure constant and m μ the muon mass. The QED kernelL has been computed in the continuum [34] 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 [27] with = 0.40 throughout. 3 The tensor i is a Euclidean hadronic four-point function with one of its vertices weighted linearly in one of its coordinates, 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. 3 We have also investigated the effect of choosing = 0 and found it to make very little difference.
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 four-point 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. [20] 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. [20] 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.

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 [35,36]. 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 [36], to be roughly doubled in order to include the D s loop. Taking into account the charge factor of [20], Appendix A 4 ), 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. 5 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,36]. All in all, these considerations finally lead us to expect an order of magnitude of (−0.3 ± 0.3) × 10 −11 .
In addition to the short-distance effect estimated above, the charm-light disconnected diagrams also involve a longerdistance 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, 6 acting effectively like a local gauge-invariant gluonic operator from the point of view 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 [37]. This effective Lagrangian L (c) 2γ 2g has been calculated long ago [38]; 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 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. [39], while the TFF normalisation amounts to |F η γ γ (0, 0)| 0.34 GeV −1 (see for instance [40]), 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 [41], 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 7 since it is weighted by 7 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 .
, 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.

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. [27]: in Appendix A we investigate the applicability of our QED-kernel implementation at particularly heavy scales by comparing the lepton-loop contribution to a LbL μ to the known analytical expression [35]. While the agreement 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.

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 symmetrised-conserved 8 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 8 A definition of the local and the symmetrised-conserved current can be found for instance in Ref. [42]. coefficients the same) we will only use f (|y|) along the lattice direction (1,1,1,1) for all ensembles. Figure 2 shows an example of our data for the A653 ensemble. The integrand is steeply-peaked 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 (trapezoidrule) 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.

Mass-dependence of the connected contribution
The results are given in Tables 3, 4, 5, 6, 7 and 8 of Appendix B and summarized in Fig. 4.
Expectations are that a μ scales with m 2 μ /m 2 heavy [21][22][23], 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 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.

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 Sect. 3. Indeed, this can be seen by the saturation of the tail of the lattice integrand (Figs. 2, 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 Table 1; 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 [43], where we promote each lattice spacing to a fitparameter 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 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 [21][22][23]. It is natural to first consider the η cmeson 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 ηc is the estimate from Ref. [3] 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 M S-mass of the charm quark [44], 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 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 fulfilled (see Table 9 in Appendix C and Fig. 5): 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 finitelattice-spacing 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 Table 9). 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 Our error on this quantity is entirely dominated by the systematic error from our modeling of its dependence on a and M η c .

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 charmquark-mass independent mass-shift within a given window of M η c . We estimate the mass-shift using the M S charm-quark 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.

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 [20], and partly on the arguments laid out in Sect. 3. 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 charmlight 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. [30]. This tuning is suboptimal for our present purposes, an aspect we return to below. We will also use the Z c V values from [30], 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. [27].
A plot of the partially-integrated charm-light and charmcharm 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 partiallyintegrated 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 Sect. 3, we start this fit between 0.4 and 0.5 fm. Table 2 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 charmlight 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 [29]. 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 valencequark 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 D takes the value (M D ) phys . Assuming that the disconnected contribution is roughly proportional to 1/M 2 D , we multiply our continuum-extrapolated result obtained at We neglect the charm-charm contribution as its contribution is far smaller than our final error for the charm-light.

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 Sect. 5 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. This estimate also covers the (u, d, s) valencequark mass dependence present in the charm-light disconnected contribution: as discussed in Sect. 3, for that topology we expect a short-range contribution (with a reach of order (2m c ) −1 ), plus a longer-range contribution with a very suppressed amplitude coming from the exchange of (u, d, s) isoscalar mesons. The short-range part is again expected to depend on the light-quark masses mainly via the D meson masses as estimated above, while the longer-ranger part is simply too small to affect our estimate.
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 [20] obtained with dynamical (u, d, s) quarks yields a fully non-perturbative determination of a HLbL This concludes our first-generation calculation of hadronic light-by-light scattering in the muon (g − 2). 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 [45], as well as the QDP++ library [46]. 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-bylight contributions in the continuum [27,47,48] at the onepercent 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 [34], it is important for the goal of this paper to test how robust this implementation is for computing contributions from more massive particles (Fig. 8).
As an example of a calculation performed entirely in the continuum, we quote the result we obtain with the kernel L ( =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. 9 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 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 continuum-extrapolated 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.

Appendix C: Fit results for the connected contribution
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 9.