Forward light-by-light scattering and electromagnetic correction to hadronic vacuum polarization

Lattice QCD calculations of the hadronic vacuum polarization (HVP) have reached a precision where the electromagnetic (e.m.) correction can no longer be neglected. This correction is both computationally challenging and hard to validate, as it leads to ultraviolet (UV) divergences and to sizeable infrared (IR) effects associated with the massless photon. While we precisely determine the UV divergence using the operator-product expansion, we propose to introduce a separation scale $\Lambda\sim400\;$MeV into the internal photon propagator, whereby the calculation splits into a short-distance part, regulated in the UV by the lattice and in the IR by the scale $\Lambda$, and a UV-finite long-distance part to be treated with coordinate-space methods, thereby avoiding power-law finite-size effects altogether. In order to predict the long-distance part, we express the UV-regulated e.m. correction to the HVP via the forward hadronic light-by-light (HLbL) scattering amplitude and relate the latter via a dispersive sum rule to $\gamma^*\gamma^*$ fusion cross-sections. Having tested the relation by reproducing the two-loop QED vacuum polarization (VP) from the tree-level $\gamma^*\gamma^*\to e^+e^-$ cross-section, we predict the expected lattice-QCD integrand resulting from the $\gamma^*\gamma^*\to\pi^0$ process.


Introduction
The long-standing discrepancy between theory [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] and experiment [21,22] for the muon g −2 has recently been challenged by several precision lattice QCD calculations of the HVP contribution (cf. Fig. 1) from intermediate hadronic distance scales [23][24][25][26][27][28]. One of the lattice-QCD based calculations has already reached a subpercent-level of precision for the full leading-order HVP contribution [24], thus becoming competitive with the data-driven dispersive method, which has traditionally been used to evaluate this contribution. At this level of precision, care must also be taken of the leading isospin-breaking corrections, both the strong isospin-breaking effect stemming from the unequal u and d quark masses, and the e.m. effect arising from the quarks carrying electric charges, as shown by the second diagram in Fig. 1. These effects are taken into account by lattice collaborations (e.g., Ref. [23,24,29]); however, few stringent cross-checks are possible at present. First of all, these effects depend on the precise point in the parameter space of isospinsymmetric QCD, which is not exactly the same in different calculations. Furthermore, it has not been possible to rigorously compare the size of these effects to phenomenological predictions, partly because the (QED) radiative correction to the HVP is divergent if one does not account for the counterterms associated with the quark masses and the strong-coupling, whose finite parts depend on the conventional choice of the 'physical point' in isospin-symmetric QCD [30]. Here we propose a computational strategy that enables much more direct comparisons between lattice QCD and hadron phenomenology. In its simplest incarnation, the idea is to add and subtract a Pauli-Villars term [31] to the photon propagator 1 , where k is the Euclidean four-momentum and Λ is a typical hadronic scale. The first term leads to a UV-finite effect on the HVP and is sensitive to long-distance contributions such as the π 0 γ and ηγ channels; it can be treated analogously to the HLbL contribution to g −2 by using coordinate-space methods [33][34][35], whereby power-law effects due to the internal photon propagators are avoided. The second can be treated entirely in lattice regularization by having the photon field defined on the same lattice as the QCD fields [30]. However, since a photon mass is now present, no issue with the photon zero-mode arises, nor do power-law finite-size effects occur. We return to this aspect in section 6 but remark here that a number of different methods have been used in the extensive literature on incorporating the coupling of quarks to photons into lattice QCD calculations (see [36][37][38][39][40][41][42][43][44][45][46][47][48][49] for a representative set of publications). How then can one predict the leading QED correction to the HVP with a UV-regularized photon propagator in place? We shall express it through the forward HLbL amplitude [50][51][52], as shown in Fig. 2. As has been noted [52], the connection between the forward HLbL amplitude and the e.m. correction to the HVP bears a strong resemblance with the Cottingham formula [53][54][55][56][57], which expresses e.m. mass splittings in terms of the forward Compton scattering amplitude. The analogy becomes apparent if one views light-by-light (LbL) scattering as Compton scattering off a photon. However, not much practical use has so far emerged from this connection. Perhaps the main reason is that the insertion of two standardqqγ vertices leads to a divergence, requiring the insertion of O(e 2 ) counterterms to cancel it. From the standpoint of the HLbL amplitude, the divergence appears due to the forward HLbL falling off too slowly (as 1/k 2 ) for one of the incoming photon momenta becoming large. However, the first term on the right-hand side of Eq. (1) amounts to a UV-regularization of the photon propagator, and in this case the integral over the forward HLbL amplitude yielding an e.m. correction to the HVP becomes finite. Therefore, with sufficient knowledge of the forward HLbL amplitude, obtained either by using the dispersive sum rules [58] or direct lattice calculations [59], one can make a definite prediction for this correction. The comparison can even be done in a more differential way, at the integrand level, as we shall illustrate. Also, accumulated knowledge on the HLbL amplitudes, for example concerning the relative importance of the different quark Wick-contraction topologies [60,61], can be usefully applied to the leading QED corrections to the HVP.
The rest of this paper is organized as follows. We start in the continuum, deriving in section 2 the relation between the forward HLbL amplitude and the e.m. correction to the HVP. In section 3, this relation is tested against known results in a pure QED setting. Section 4 contains a derivation of the divergence that develops when the UV-cutoff Λ is sent to infinity and indicates in which flavor combinations the divergence partly cancels. In section 5, we then come to the prediction for the π 0 exchange in the forward HLbL amplitude, and thereby to its contribution in the e.m. correction to the HVP and ultimately in the muon (g − 2). We then formulate a computational strategy in section 6 for computing the leading isospin-breaking effects to the HVP in lattice QCD. The section also discusses aspects of the coordinate-space method and presents the integrand corresponding to the π 0 exchange contribution. Finally, section 7 summarizes our findings and offers an outlook into further possible applications of this work.

A Cottingham-like formula for the radiative correction to the HVP
The first correction ∆Π(Q 2 ) to the leading HVP 2 Π e 2 (Q 2 ) can be written in the form where Λ is a UV-regularization parameter. We begin by establishing a formula for Π 4pt (Q 2 , Λ), named in this way because it involves the four-point function of the e.m. current, which at the same time provides the quantum field-theoretic definition of the LbL scattering amplitude. The second term, Π ct (Q 2 , Λ), consists of the required counterterms and the strong-isospin breaking contribution. While its precise form is not needed here, more details will be given in section 6. The LbL scattering amplitude M µ 1 µ 2 µ 3 µ 4 depends on the four-momenta of the incoming (q 1 , q 2 ) and outgoing (q 3 , q 4 ) photons. The forward kinematics correspond to q 1 = q 3 ≡ k and q 2 = q 4 ≡ q, see Fig. 2. Contracting the photon-line 1 with 3, we obtain a contribution to the VP tensor: where the factor of one-half is the symmetry factor; in square brackets is the Feynman-gauge photon propagator, regulated at the scale Λ, for instanceà la Pauli-Villars, Due to gauge invariance, the VP tensor has the following general form, and hence its scalar part can be expressed as: where is the traced LbL amplitude. The latter is a scalar function of three invariants: k 2 , q 2 , and ν ≡ k ·q. It is even in ν and symmetric under the interchange of k and q. We shall write it as M(ν, K 2 , Q 2 ), where K 2 = −k 2 and Q 2 = −q 2 will further be assumed to be positive, i.e., the photons are spacelike.
Introducing the helicity LbL amplitudes as with ε µ λ (q) the photon polarization vectors, the traced amplitude can be written as [64]: where For spacelike photon virtualities, the optical theorem relates the imaginary part of these amplitudes to a γ * γ * -fusion cross section [58,Eq. (16)], so that where σ = 4σ T T − 2σ T L − 2σ LT + σ LL , and X = ν 2 − q 2 k 2 . Furthermore, the analytic properties of the ν-dependence warrant a dispersive representation. Since all relevant LbL amplitudes are even in ν and require one subtraction [58, Section II. C], the dispersion relation takes the form where we are free to choose any subtraction pointν, and ν thr. is the lowest particle production threshold. For example, in QED, ν thr. = 1 /2(K 2 + Q 2 ) + 2m 2 e is the threshold for e + e − production; see Appendix A.1 for further details.
The dispersive representation justifies the Wick rotation in the evaluation of Eq. (5) and we obtain the Cottingham formula analogue: We refer the reader to Eq. (A.2) in appendix for the dispersive form, which, up to one subtraction, expresses this contribution in terms of γ * γ * -fusion cross sections. As indicated in Eq. (2), this contribution must be combined with the appropriate counterterm, to which we return in section 6, in order to obtain the first correction ∆Π(Q 2 ) to the HVP. At this point, let us briefly comment on the flavor structure of the HLbL amplitude M, particularly regarding to isospin, which plays an important role at low energies. The e.m. current carried by the quarks contains both an isovector and an isoscalar component. The LbL amplitude can be written as the sum of the three partial contributions where (i) all four currents are isovector; (ii) all four currents are isoscalar; and (iii) in one pair of currents, both are isovector, while in the complementary pair, both are isoscalar, and one sums over all six possible pairings. Pole contributions of isovector mesons such as the pion only occur in the third contribution, while isoscalar-meson exchanges appear in all three contributions.

Reproducing the two-loop QED vacuum polarization
In order to test our Cottingham analogue, Eq. (12), we apply it to the QED VP: we expect the one-loop LbL amplitude to provide the two-loop VP, see  Substituting the one-loop LbL amplitude (cf. Appendix A.1) into Eq. (12) yields a complicated expression, which we only show here in the expanded form: up to terms that vanish for Λ → ∞. Hereafter m stands for the lepton mass appearing in the loops. In this calculation we have in fact adopted a simpler, momentum-cutoff regularization: In the present context, this form of regularization is equivalent to Pauli-Villars regularization, up to terms suppressed by 1/Λ 2 .
The counterterm, Π ct (Q 2 , Λ), can be obtained by applying the standard rules of renormalized perturbation theory, for which we use the Pauli-Villars regularization, see Appendix A.3: where κ = Q/m . Altogether [cf. Eq. (2)], we obtain the following result for the small-Q 2 expansion of the two-loop VP: Note that: with Π 4pt (0, Λ) given in Eq. (A.13). To check this result, we use the well-known dispersion relation: together with the well-known O(e 4 ) imaginary part, first computed in 1955 [65] as well as in [66], given for instance in [67] (Eq. (5-4.200) on p. 109) and in [68]. Using this well-known result, we reproduce Eq. (15) obtained via the Cottingham formula. We have also checked this numerically for arbitrary Q 2 .

Evaluation of the fourth-order vacuum polarization contribution to the muon (g − 2)
We now test the Cottingham-like formula beyond the small-Q 2 expansion and compute the fourth-order VP contribution to the muon (g − 2) using numerical integration. In order to avoid numerical instabilities emerging in loop integrals, that in the case of the virtual LbL amplitude we found out to persist in all familiar packages for one-loop numerical integration, we choose the more elegant way and calculate the LbL amplitude M(ν, K 2 , Q 2 ) itself via the dispersive approach; see Eq. (11). The imaginary part Im M(ν, K 2 , Q 2 ) at the tree level in QED, as well as expressions for the subtraction term M(ν, K 2 , Q 2 ) for the two choices of subtraction pointν 1 = 0 andν 2 = KQ are provided in Appendix A.1.
We extracted the 'Thomson limit' Π 4pt (0, Λ) from the first term in the series expansion of M, and after performing the integration over x, we arrive at Eq. (A.13). In particular, this quantity is finite for a fixed value of Λ. From here, the O(e 6 ) contribution of the regulated fourth-order QED VP to the anomalous magnetic moment can be computed using the general formula [69][70][71]: where the kernel function is given by with m µ the muon mass. The integral over the kernel function alone gives the Schwinger term, ∆a µ = α/2π. Therefore, substituting the renormalized VP, we obtain: This means that the subtraction term in Π 4pt (Q 2 , Λ) yields the contribution − α 2π Π 4pt (0, Λ) to a VP µ . Results obtained with a Pauli-Villars regulator and Λ = 3m µ are provided as a function of the lepton mass m appearing in the VP in Fig. 4. Furthermore, adding the contribution of the appropriate counterterm and subsequently taking the limit Λ → ∞, we obtain the full O(e 4 ) QED VP to a µ . For a muon in the VP loop, we obtain: The result numerically agrees with Ref. [72].

Forward LbL amplitude at large virtuality from the Operator Product Expansion
Since the forward HLbL amplitude M(k, q) is finite, the divergence of Eq. (5) as Λ → ∞ can only arise from performing the k integral. The question is then, what is the large-k behaviour of M(k, q) for fixed q. This is a typical application for the Operator Product Expansion (OPE). In this section we work in Euclidean space and our starting point is the e.m. current carried by the quarks, in units of the positron charge, being given by V em . We recall that the HLbL amplitude is directly related to Π µνσλ [59]. In particular, for the forward amplitude Eq. (6), the connection reads where the scalar products on the left-hand side are Euclidean. The large momentum k "forces" the two vertices x and y to come close together. From a power-counting perspective, it is the dimension-four operators that can cause a logarithmically divergent behaviour in Eq. (5), since they contribute as O (4) /k 2 . It is then only necessary to know their Wilson coefficients to order α s included, since α s (k 2 ) 2 O (4) /k 2 multiplied by a photon propagator already yields a UV-finite integral.
Note that the two indices of the vector currents are contracted with each other (thus cancelling the axial current contribution in the OPE), and that we may average the result over the direction of k, given that we are interested in subsequently integrating over k in Eq. (5). The result of the OPE can then only contain operators with vacuum quantum numbers.
From a different perspective, the divergence resulting from the integral over the photon momentum k must be removable by the available counterterms of the theory. Moreover, since vector currents do not renormalize in QCD, the relevant counterterms are only those associated with the parameters of the theory, which are the gauge coupling and the quark masses. These parameters are respectively associated with the operators G a αβ G a αβ and m fψf ψ f for each quark flavor f . The Wilson coefficients of scalar operators appearing in the OPE of QCD currents have been calculated a long time ago [73]. We report only the result up to the order required for our purposes 3 , Interpreting the gluonic operator in terms of the (renormalization group invariant) trace anomaly Thereby we arrive at the following prediction for the asymptotic large-k 2 behaviour of the fourpoint amplitude . .} are the quark electric charges. At this point, we keep the currents V λ (0) unspecified, in particular in their flavour structure. The effect of inserting the mass operator into a correlation function is to differentiate the latter with respect to the quark mass, while the effect of the trace anomaly on a renormalization-group invariant correlation function of mass-dimension n is to differentiate with respect to all scales on which the correlation function depends, where y and z are space-time coordinates. We now set V . Let H λσ (z) be the kernel yielding the leading-order subtracted VP when integrated over with the correlator e 2 V em σ (z)V em λ (0) (see Eqs. 56 and 57 below). Acting with the linear operator on both sides of Eq. (26), we conclude that the asymptotic large-Λ behaviour of the four-point amplitude contribution to the fourth-order VP is given by 4 In these conventions, b0 = 1 (4π) 2 (11 − 2 3 n f ) and b1 = 1 (4π) 4 (102 − 38 3 n f ).

Explicit OPE calculation at leading order
Consider then the OPE of two vector currents at leading order, with S(x) the position-space fermion propagator. Eventually one finds We have already noted the effect of the mass-operator insertion in terms of differentiating the correlation function with respect to the quark mass. To understand the effect of inserting the 'equation of motion' (EOM) operator appearing in brackets in Eq. (32), imagine multiplying the Euclidean quark action by λ, S E (λ) = λψ(D + m)ψ. In the Euclidean path integral, we can take expectation values A λ using exp(−S E (λ)) as weight: it simply means that each quark propagator contains an additional 1/λ factor. Thus, if computing A λ involves n p propagators, On the other hand, the same derivative can be expressed as Thus the insertion of the EOM operator simply multiplies the observable with the number of propagators n p needed to compute it. Thus the leading contribution of x and y being close together in the vector four-point function (V µ ≡ψγ µ ψ) is We now move to the case of x and y simultaneously being in close vicinity of a third current at position z. In terms of Wick contractions, the relevant case is where the connecting point is z, i.e. there are propagators (y → z → x) or (x → z → y). The cases where the connecting point is x or y do not contribute to the O(1/k 2 ) behaviour. A straightforward if somewhat tedious calculation then gives The same contribution appears when (x, y) are close to the origin, thus doubling this contribution in the four-point function of the vector current. Altogether, from Eqs. (35) and (36), we then find in leading order of the OPE The terms not leading to mass-derivatives cancel, and only the mass-derivative of the vector twopoint function determines the large-k 2 asymptotics of the forward LbL amplitude. Thus we have reproduced the very first term in Eq. (26). Acting on both sides of Eq. (37) with Eq. (29) leads to the first term in Eq. (30) (with n f = 1 and Q f = 1). The leading-order calculation above is equally valid for the QED as for the QCD four-point function. In the pure QED context, one easily verifies with the help of Eqs. (A.14) and (A.15) that the textbook O(e 2 ) mass counterterm removes the log(Λ) term in Π 4pt predicted by the OPE. We have thus verified Eq. (30) in the pure QED case: given that Π e 2 (Q 2 ) given in Eq. (A.16) only depends on Q 2 /m 2 , the effect of inserting the trace anomaly, proportional to (2Q 2 ∂ ∂Q 2 + m ∂ ∂m ), would cancel.

Cottingham-like formula for the isovector contribution to HVP
Starting from Eq. (26) with the isovector component of the e.m. current, the same steps lead to the analogue of Eq. (30), with Π e 2 now replaced by the isovector contribution Π 33 e 2 (Q 2 ) to the leading HVP, and Π γγ33 4pt (Q 2 , Λ) the corresponding QED correction to that contribution. The same steps once again could be taken with the charged isovector currents which leads to the quantities Π −+ e 2 (Q 2 ) and its QED correction Π γγ−+ 4pt (Q 2 , Λ). Taking the difference of Eq. (26) obtained once with the choice of currents from Eq. (38) and once with the choice of Eq. (39), one finds that all terms on the right-hand side cancel. This is clear for the contribution from the insertion of isoscalar operators, whose correlation function is identical with two members of the same isospin multiplet. For the isovector mass insertion (ūu −dd), G-parity ensures that this insertion vanishes separately in the neutral and in the charged isovector channel. In other words, all counterterms from the action cancel in this difference 5 , as has been noted in Ref. [75], which contains an exploratory lattice-QCD calculation of this quantity. However, care must be taken of the fact that, unlike the cases considered so far, the currents of Eq. (39) are not gauge invariant with respect to QED. Therefore this case requires further study. We note that, in the radiative corrections to the leptonic decay of a charged pion [76], the e.m. correction to a charged-current correlator represents one of several contributions.
Other flavor cases may be of interest, in particular the correlator between the isovector and the isoscalar components of the photon, which vanishes in isospin-symmetric QCD in the absence of quark electric charges. In this case, the action counterterms do not vanish altogether, though only the isovector mass insertion (ūu −dd) contributes in the analogue of Eq. (30).

The π 0 -exchange contribution
In this section, our goal is to present the form of the π 0 -exchange contribution to the forward HLbL amplitude M, cf. Fig. 5, since it is the longest-range contribution. As in the previous section, we work in Euclidean space. We define the Fourier transform of the four-point function of this current as in Eq. (22) and the forward amplitude is obtained as in Eq. (23). The O(e 4 ) contribution to the polarization tensor with a regularized internal photon propagator then reads For the π 0 -exchange contribution, proportional to the square of its transition form factor F, we have We perform the angular integration by using the Gegenbauer polynomial expansion of propagators (see for instance [77]), and the final expression is Z m |q|,|k| = 1 2|q||k| This expression, once inserted into Eq. (18), can be viewed as the VP analogue of the Jegerlehner-Nyffeler relation for the π 0 contribution to HLbL scattering in the muon (g − 2) [78,79]. In the present case, the kinematics are simpler, and correspondingly ∆a π 0 µ takes the form of a two-rather than three-dimensional integral for a yet to be specified transition form factor F. The unsubtracted VP is UV-finite for fixed Λ, while the subtracted VP Π 4pt (q 2 , Λ) remains UV-finite for Λ → ∞, unlike in the full QCD case, as we have seen in section 4, when short-distance contributions from quarks are taken into account.
As an example, for the VMD parameterization of the transition form factor, one obtains, near the chiral limit, the singular behaviour For a pion mass which is still heavy relative to the muon mass, the contribution reads ∆a µ αm 2 µ 3π Π 4pt (0). Parametrically, this contribution behaves similarly to the π 0 contribution to the HLbL contribution to a µ [51,80], except that in the latter case the chiral logarithm enters quadratically. Numerically, with F(0, 0) = (4π 2 f π ) −1 and f π = 92.4 MeV, m V = 0.77549 GeV, and the physical π 0 mass one obtains from Eq. (42) with the QED kernel (18) the following contribution to a µ , ∆a π 0 µ = 0.370 × 10 −10 .
This result agrees with the value given in [51]. We note that the result is more than an order of magnitude smaller than the contribution of the e + e − → π 0 γ channel in the dispersive representation of a VP µ 6 , however the quantity ∆a π 0 µ computed here is not precisely the same. We finally remark that the master relation Eq. (42) applies equally well to the other pseudoscalar mesons, notably the η and η .

Electromagnetic correction to the HVP in lattice QCD: a computational strategy
The general structure of Eq. (2) also applies to the calculation of the isospin-breaking contribution to the leading HVP Π e 2 (Q 2 ) (from QED and strong isospin breaking) in lattice regularization, in which case the inverse lattice spacing 1/a plays the role of the UV cutoff. To be more specific, we note that lattice QCD admits (N f + 1) bare parameters, i.e., the SU(3) gauge coupling and the N f quark masses, which we assemble into a vector b lat . When correcting the isosymmetric theory for isospin breaking, the bare parameters must be readjusted. The shifts δb lat i in the bare parameters are determined by requiring that the theory with isospin breaking reproduces (N f + 1) suitable experimental observables [30]; typically, the masses of hadrons which are stable in the absence of weak interactions, lat |h is given by the forward matrix element of the operator conjugate to parameter b lat i . Thus, given a lattice calculation of M iso , M lat 4pt (a; 0) and the matrix J lat , the vector δ b lat is obtained by solving a linear system.
We are now in a position to write the lattice-regularization analogue of Eq. (2) for the subtracted HVP as ∆Π(Q 2 ) = lim where the counterterm has the form and Π lat 4pt (Q 2 , a; M γ ) denotes the e.m. four-point function contribution to the HVP, computed (in general) with an internal photon of mass M γ . Massive QED has previously been used to control photon zero modes in finite volume, with the physical limit M γ → 0 taken after extrapolating to infinite volume [40,48,82,83]; our approach, by contrast, is to keep Λ = M γ fixed and use it for separating long-range contributions from UV-divergent ones. Based on the identity of Eq. (1) for the photon propagator with a fixed Λ ∼ 400 MeV, we propose to perform the following decompo- Up to the subtraction at Q 2 = 0, the function Π 4pt (Q 2 , Λ) is the same as in Eq. (2). Since Λ plays the role of the Pauli-Villars UV regularization scale, the continuum limit a → 0 can be taken for this quantity. Similarly, M 4pt (Λ) represents the e.m. hadron-mass corrections computed with the Pauli-Villars regulated photon propagator. The continuum limit can also be taken in this case, since the same OPE, as reviewed in section 4, determines the asymptotic behaviour of the forward Compton amplitude on hadron h [74]. Given a choice of M γ , each term on the right-hand side of Eq. (50) is affected by rather different systematics on the lattice and is meant to be evaluated separately. The same observation applies to the two terms on the right-hand side of Eq. (51). The UV-finite part Π 4pt (Q 2 , Λ) of Eq. (50) receives long-distance contributions due to the long-range photon propagator. A coordinate-space representation free of power-law finite-volume effects is presented below in subsection 6.1. The result can be compared to an evaluation based on our Cottingham-like formula, Eq. (11) and Eq. (12). For the second term of Eq. (50), one possible expression in the time-momentum representation [63] is where V em ν (y) and V em µ (x) are discretized as conserved currents and the simplest form of the photon propagator (in Feynman gauge) is and the tadpole terms ensure the transversality of the four-point function with respect to contracting it withk µ ork ν . Similarly, M lat 4pt (a; M γ ) can be determined by well-established methods where the (now massive) photon is treated as part of the finite-volume lattice field theory.
We now briefly discuss how to compute the hadronic mass shifts M 4pt (Λ) with a long-range, but Pauli-Villars regulated photon propagator. As for Π 4pt (Q 2 , Λ), it is possible to avoid powerlaw effects in the volume [84] by using coordinate-space methods and, additionally, by explicitly correcting the elastic contribution of the forward Compton amplitude for finite-volume effects. As a slight variation to the concrete proposal in [84], this correction could be done with the help of a separate calculation of the e.m. form factor(s) of the hadron whose mass correction is being computed. These methods could also be applied to the difference of M 4pt,h (Λ) between proton and neutron, a quantity that could be compared to predictions based on the original Cottingham formula.
We remark that the currently most frequently used formulation of lattice QCD coupled to photons consists in removing the photon zero-mode in every time-slice [85]. The corresponding finite-size effects on the HVP have recently been investigated and found to be parametrically of order 1/L 3 , and numerically small in the framework of scalar QED [86]. Another recent investigation provides a systematic analysis of various finite-size effects beyond the pointlike approximation of hadrons, in particular of pseudoscalar mesons masses [87] (see also references therein). As an alternative method, first results (primarily on hadron masses) based on simulating QCD+QED with C * boundary conditions have recently been presented [88].
In conclusion, while the numerical practicability of the presented method remains to be demonstrated, we have established that it is possible to avoid power-law finite-volume effects altogether in computing ∆Π(Q 2 ) on the lattice. How large the discretization errors on this quantity are at finite lattice spacing with the method proposed above will need to be explored in practice. Here we remark that the Pauli-Villars regularization of the photon propagator is only one of many posssible choices. For instance, with the decomposition which amounts to a 'double Pauli-Villars' regularization of the photon propagator, the same strategy as described above can be carried out, now with the expression in brackets in Eq. (54) falling off as fast as 1/k 6 at large k 2 .
6.1. Coordinate-space representation of Π 4pt (Q 2 , Λ) free of power-law finite-size effects A Euclidean coordinate-space expression for the subtracted HVP is where the leading contribution isΠ e 2 ;σλ (z) = e 2 V em σ (z) V em λ (0) , the relevant (Q-dependent) coordinate-space kernel H λσ (z) was derived in [89] (Sect. II.B.2) and we have abridged z ≡ d 4 z. Expanding a QCD correlation function to second order in the e.m. coupling leads to the insertion of the product of two e.m. currents, whose relative positions are weighted by the internal photon propagator. Thus, using Feynman gauge for the latter, we arrive at the expression for the regulated contribution to the subtracted HVP. The Pauli-Villars regulated photon propagator in position space reads which is only logarithmically divergent for x 2 → 0. Here K 1 is the modified Bessel function of the second kind. We note a close analogy of expression (57) with the master relation used for the HLbL contribution to the muon (g − 2) in Refs. [35,90,91]. Similarly, a VP µ can be obtained from Eq. (57) by replacing the kernel H λσ by the appropriate one given in Sect. II.B.3 of Ref. [89]. The main feature of our proposal is that no IR-regularization of the photon propagator is needed. Thus finite-size effects are expected to be on the order of exp(−m π L/2), as in the case of the HLbL contribution [90].

Computing a HVP
µ : π 0 -exchange contribution to the coordinate-space integrand By Fourier-transforming the polarization tensor associated with Eq. (42) (see Eq. (40)), one obtains the O(e 4 ) contributionΠ 4pt;σλ (z, Λ). After insertion into Eq. (56) and contraction of the indices, the integrand is a scalar function of |z|. We illustrate the integrand of this last scalar integral in the following. Alternatively to the proposed position-space approach, one can reach the widely used time-momentum representation (TMR) [63] by Fourier-transforming the polarization tensor at vanishing spatial momentum only with respect to q 0 .
Choosing the same parameters as in section 5 for the VMD parameterisation of the transition form factor, we obtain the integrands as functions of R ≡ |z| (and R ≡ z 0 for the TMR case) as shown in Fig. 6, where the results with the original [89] position-space kernel (kerO) and the TMR are compared, calculated in the continuum and infinite volume at two different Pauli-Villars masses Λ = 3 m µ and 200 m µ ; recall that the π 0 -exchange contribution by itself remains finite as Λ → ∞. Both integrands are rather long-range, an observation which implies a certain difficulty for lattice calculations if the O(e 4 ) contribution is to be computed with good relative precision. We remark that the integrand displayed in Fig. 6 corresponds to the sum of all Wick contractions contributing to the four-point function of the e.m. current, but, using the results in appendix A of Ref. [35], it would be fairly straightforward to adapt the prediction to individual Wick contractions of the quark fields.
As a consequence of the Ward-Identity of the vector current, a term ∂ λ [z σ F (|z|)] can be added to the position-space kernel H λσ (z) without changing the integrated result of Eq. (56) in infinite volume [92]. With a judiciously chosen subtraction, one can make the z-integrand in Eq. (57) more peaked in the small-|z| region. As in practice, a calculation on the lattice is limited by the degrading signal-to-noise ratio when the arguments of the correlator are far apart in position-space, the possibility of reshaping the integrand makes the position-space representation appealing. Such a technique has also been used for lattice determinations of the HLbL scattering contribution to a µ [90,93]. Figure 6 thus also shows the result of improving the kernel (kerM6) to make the integrand shorter-range; details of its construction are given in Appendix B. For the Λ = 3 m µ case, the partially-integrated a µ (R) obtained with kerM6 already reaches about 70% of its final value a µ (∞) at R = 2 fm, but only about 50% with the TMR. For Λ = 200 m µ , the benefit becomes even more apparent: 95% with kerM6 and merely about 65% with the TMR. We thus expect that the position-space method with an improved kernel should offer a good opportunity to compute the e.m. correction to a HVP µ on the lattice, with better controlled finite-volume effects.

Conclusion
We have written a Cottingham-like formula for the leading QED correction to the HVP, mainly in terms of the traced forward HLbL scattering amplitude. While the latter is a physical amplitude, when contracting an incoming with an outgoing photon line, the corresponding momentum integral diverges logarithmically in the UV. The required counterterms that remove that divergence have been worked out in section 4 with the help of the OPE: they involve the derivatives with respect to the quark masses and the virtuality of the leading HVP. The finite part of the counterterms, however, depends on the precise choice of the point in the parameter-space of isospin-symmetric QCD around which the isospin-breaking effects are computed.
At present, it appears that the most promising application of the Cottingham-like formula is to implement it with a finite regulator on the order of a few hundred MeV. The regularization amounts to replacing the internal photon line by a Pauli-Villars regulated propagator, or any other convenient form of propagator regularization. The leading QED correction to the HVP with such a regularized photon propagator in place can be computed in lattice QCD using coordinatespace techniques similar to the calculation of the HLbL contribution to the muon (g − 2), without incurring power-law finite-size effects. When working on very large lattices, as realized in masterfield simulations [94,95], these techniques are particularly natural [96]. Irrespective of whether one uses this new approach or the more established method involving the removal of the spatial zeromode of the photon, a direct comparison becomes possible between the lattice-QCD calculation and the prediction based on the Cottingham-like formula. For the latter, we recall that the traced forward HLbL amplitude can be represented dispersively in terms of the γ * γ * → hadrons fusion cross-section, up to one subtraction term. The complementary part, i.e. the second term in Eq. (1), involves a massive photon propagator in the case of the Pauli-Villars regularization choice. For this part, the lattice provides a natural UV-regularization, and a prescription for handling the photon zero mode on a finite lattice is then no longer needed, owing to the photon mass.
Finally, in the context of the original Cottingham formula, it could be interesting to compute the e.m. contribution to the proton-neutron mass difference with a regularized photon propagator, which could then be compared in detail and without scheme uncertainty to predictions based on the (rather mature) dispersive treatment of the forward Compton amplitude (see [56] and Refs. therein).

Acknowledgements
We acknowledge the support of Deutsche Forschungsgemeinschaft (DFG) through the research unit FOR 5327 "Photon-photon interactions in the Standard Model and beyond -exploiting the discovery potential from MESA to the LHC" (grant 458854507), and as part of the Cluster of Ex Here we give more details on how the two-loop QED contribution to vacuum polarization is reproduced via the Cottingham-like formula. We begin with Eq. (12), which by a variable change, ν = KQx, is cast into: where Λ is the scale regularizing the integral over K 2 in the ultraviolet. We keep the regularization implicit throughout this appendix. Substituting the dispersive representation of the LbL amplitude and the optical theorem, we obtain where X = ν 2 − Q 2 K 2 , andν is the subtraction point (below we useν = 0 andν = KQ and verify that the results are equivalent). We next provide the ingredients needed to evaluate the vacuum polarization to two-loops in QED, including the necessary counterterms.
Appendix A.1. Polarized γγ fusion cross sections We start with the tree-level cross sections for the QED process γ * γ * → ( denote spinor QED fermions) with polarized virtual photons. The conventions are the same as in Refs. [58] and [97]: where s = (k + q) 2 = 2ν − K 2 − Q 2 , ν = k · q; Q 2 = −q 2 and K 2 = −k 2 are the spacelike photon virtualities. The threshold energy in this case is given by ν thr. = 2m 2 + 1 /2(K 2 + Q 2 ). The cross sections corresponding to the fusion of two polarized photons, either transverse (T ) or longitudinal (L), read: The optical theorem connects the absorptive cross sections with the imaginary part of the corresponding forward LbL amplitudes, Im M P P = 2 √ X σ P P , P, P ∈ {L, T }.
The total cross section requied in the Cottingham formula is given by 8 The vacuum polarization at Q 2 = 0, needed for renormalization, is, in general, given by: In the case of two-loop QED, we obtain: which is finite upon any regularization set by Λ.

Appendix A.3. QED counterterm in the Pauli-Villars regularization
The counterterm can be obtained in the following way where in the on-shell renormalization scheme 9 and in Minkowski-space notation After the subtraction, the functions H i , henceforth referred to as form factors, are modified tō As the coordinate-space formulation for a µ is obtained by a simple substitution of the e.m. kernel appearing in Eqs. (56) and (57), we will not introduce new notations in the following discussion. Throughout this appendix, the e.m. kernel and form factors refer implicitly to the ones related to a µ . To avoid the long-distance region where lattice calculations perform less well, we aim at reshaping the integral representation of a µ into a shorter-ranged one by introducing a phyiscally motivated subtraction to the e.m. kernel. In the leading-order HVP calculation, the contribution of the ρ-meson is dominates over a large distance interval. The simplest way to model the ρ-meson is to represent it as a δ-function in the spectral function (see Eq. (67) of Ref. [89]). Based on the behaviour of the correlator at large separations in this simple model, we can choose a one-parameter subtraction function The subtraction term Eq. (B.4) preserves the behavior of the kernel near the origin and at long distances. This subtraction ensures that the ρ-meson contribution from our model with a mass of M times that of the muon falls at least faster then re −M mµr log(m µ r) at large r.

Appendix B.2. Approximants for the form factors
The form factors require evaluating Meijer's G-functions, which are in general computationally costly. Fortunately, as the asymptotic behaviors of the form factors are known, one can efficiently approximate them according to the size of the argument. For the typical size of the boxes that we include in lattice calculations to approach the physical point, the maximal separation between the two vector-current insertion is ∼ 6 fm. Eventually, one would also want to have good control over the wrap-around effects of light intermediate states due to periodic boundary conditions. Accounting for these considerations, useful approximants of the form factors should be accurate up to about a distance of m µ r ∼ 13 . The approximants that we choose are of the form where the coefficients a j , p j and q