Electromagnetic finite-size effects to the hadronic vacuum polarization

TL;DR

This paper derives an analytic expression for electromagnetic finite-size effects on the hadronic vacuum polarization in lattice QED, crucial for improving muon g-2 predictions by reducing uncertainties.

Contribution

It provides the first next-to-leading order analytical formula for finite-volume electromagnetic corrections to the hadronic vacuum polarization in scalar QED, confirming universality and agreement with numerical simulations.

Findings

Leading finite-size correction is proportional to 1/L^3.

No 1/L^2 term due to charge neutrality.

Analytical results match numerical loop integrals and lattice simulations.

Abstract

In order to reduce the current hadronic uncertainties in the theory prediction for the anomalous magnetic moment of the muon, lattice calculations need to reach sub-percent accuracy on the hadronic-vacuum-polarization contribution. This requires the inclusion of $\mathcal{O}(\alpha)$ electromagnetic corrections. The inclusion of electromagnetic interactions in lattice simulations is known to generate potentially large finite-size effects suppressed only by powers of the inverse spatial extent. In this paper we derive an analytic expression for the $\mathrm{QED}_{\mathrm{L}}$ finite-volume corrections to the two-pion contribution to the hadronic vacuum polarization at next-to-leading order in the electromagnetic coupling in scalar QED. The leading term is found to be of order $1/L^{3}$ where $L$ is the spatial extent. A $1/L^{2}$ term is absent since the current is neutral and a photon far away thus sees no charge and we show that this result is universal. Our analytical results agree with results from the numerical evaluation of loop integrals as well as simulations of lattice scalar $U(1)$ gauge theory with stochastically generated photon fields. In the latter case the agreement is up to exponentially suppressed finite-volume effects. For completeness we also calculate the hadronic vacuum polarization in infinite volume using a basis of 2-loop master integrals.

Figures

Equations

$${\Pi }_{\mu \nu }(q)=\int d^{4}x{e}^{iq\cdot x}⟨0|\mathrm{T}[j_{\mu }(x)j_{\nu }^{†}(0)]|0⟩,$$

$$\hat{\Pi }\left(q^2\right)=\Pi \left(q^2\right)-\Pi \left(0\right).$$

$$\mathcal{L}=\left({\partial }_{\mu }{\varphi }^{\ast }+ie{A}_{\mu }{\varphi }^{\ast }\right)\left({\partial }_{\mu }\varphi -ie{A}_{\mu }\varphi \right)+m^2{\varphi }^{\ast }\varphi +\frac{1}{4}{F}_{\mu \nu }{F}_{\mu \nu },$$

$$(F):$$

$$(G):$$

$$(A):$$

$$(B):$$

$$(E):$$

$$(C):$$

$$(T):$$

$$(S):$$

$$(X):$$

$$\hat{\Pi }\left(q^2\right)=\frac{1}{3q_0^2}\sum _{j=1}^3\left[{\Pi }_{jj}\left(q_0,0\right)-{\Pi }_{jj}\left(0\right)\right],$$

$${\Pi }^{(1)}\left(q^2\right)=$$

$${\rho }_U\left(\mathbf{k},\mathbf{\ell },q_0\right)=\int \frac{d{k}_0}{2\pi }\frac{d{\ell }_0}{2\pi }{\pi }_U\left(k,\ell ,q_0\right).$$

$$\Delta {\hat{\Pi }}_U\left(q_0^2\right)=\left(\frac{1}{L^6}{\sum _{\mathbf{k}}}^{\prime }\sum _{\mathbf{\ell }}-\int \frac{d^3\mathbf{k}}{{\left(2\pi \right)}^3}\frac{d^3\mathbf{\ell }}{{\left(2\pi \right)}^3}\right){\hat{\rho }}_U\left(\mathbf{k},\mathbf{\ell },q_0\right),$$

$$\Delta {\hat{\Pi }}_U\left(q_0^2\right)=\left(\frac{1}{L^3}{\sum _{\mathbf{k}}}^{\prime }-\int \frac{d^3\mathbf{k}}{{\left(2\pi \right)}^3}\right)\int \frac{d^3\mathbf{\ell }}{{\left(2\pi \right)}^3}{\hat{\rho }}_U\left(\mathbf{k},\mathbf{\ell },q_0\right)+\cdots ,$$

$${\rho }_U\left(\mathbf{k},\mathbf{\ell },q_0\right)=\sum _{j=0}^{n_U}{\left(\frac{2\pi }{\left|\mathbf{k}\right|}\right)}^j{u}_j\left(\hat{\mathbf{k}},\mathbf{\ell },q_0\right)+{\overline{\rho }}_U\left(\mathbf{k},\mathbf{\ell },q_0\right),$$

$$\Delta {\hat{\Pi }}_U\left(q_0^2\right)=\frac{{\xi }_1^U\left(q_0^2\right)}{L^2}+\frac{{\xi }_0^U\left(q_0^2\right)}{L^3}+\mathcal{O}\left(\frac{1}{L^4},e^{-mL}\right).$$

$${\xi }_j^U\left(q_0^2\right)={\Delta }_{\mathbf{n}}^{\prime }\left[\frac{1}{{\left|\mathbf{n}\right|}^j}\int \frac{d^3\mathbf{\ell }}{{\left(2\pi \right)}^3}{u}_j\left(\hat{\mathbf{n}},\mathbf{\ell },q_0\right)\right].$$

$${\Delta }_{\mathbf{n}}^{\prime }={\sum _{\mathbf{n}}}^{\prime }-\int d^3\mathbf{n}.$$

$${\xi }_j^U(q_0^2)=c_j{\varphi }_j(q_0^2),$$

$${\Omega }_{\alpha ,\beta }(z)=\frac{1}{2{\pi }^2}{\int }_0^{\infty }dx{x}^2{\omega }_{\alpha ,\beta }(x,z),$$

$${\omega }_{\alpha ,\beta }(x,z)=\frac{1}{(x^2+1{)}^{\frac{\alpha }{2}}[z+4(x^2+1){]}^{\beta }}.$$

$$\Delta {\hat{\Pi }}_E(z)=$$

$$-\frac{c_0}{m^3{L}^3}(-\frac{8}{3}{\Omega }_{0,3}+\frac{32}{3}{\Omega }_{0,4}+\frac{1}{16}{\Omega }_{2,2}+\frac{10}{3}{\Omega }_{2,3}$$

$$-\frac{32}{3}{\Omega }_{2,4}-\frac{23}{128}{\Omega }_{4,1}+\frac{5}{16}{\Omega }_{4,2}-\frac{2}{3}{\Omega }_{4,3}),$$

$$\Delta {\hat{\Pi }}_C(z)=$$

$$\Delta {\hat{\Pi }}_T(z)=$$

$$\Delta {\hat{\Pi }}_S(z)=$$

$$\Delta {\hat{\Pi }}_X(z)=$$

$$-\frac{c_0}{m^3{L}^3}(-\frac{128}{3}{\Omega }_{-2,4}-\frac{16}{3}{\Omega }_{0,3}+64{\Omega }_{0,4}-\frac{11}{24}{\Omega }_{2,2}+\frac{20}{3}{\Omega }_{2,3}$$

$$-\frac{64}{3}{\Omega }_{2,4}-\frac{17}{64}{\Omega }_{4,1}+\frac{29}{24}{\Omega }_{4,2}-\frac{4}{3}{\Omega }_{4,3}),$$

$$\Delta\hat{\Pi}(q^{2})=\frac{c_{0}}{m^{3}L^{3}}\Bigg{(}\frac{16}{3}\Omega_{0,3}+\frac{5}{3}\Omega_{2,2}-\frac{40}{9}\Omega_{2,3}+\frac{3}{8}\Omega_{4,1}-\frac{7}{6}\Omega_{4,2}-\frac{8}{9}\Omega_{4,3}\Bigg{)}\,,$$

$${\Pi }_{\mu \nu }^{\text{charged}}(q)=\int d^{4}x{e}^{iq\cdot x}⟨0|\mathrm{T}[j_{+\mu }(x)j_{-\nu }(0)]|0⟩={\delta }_{\mu \nu }{q}^2{\Pi }_1\left(q^2\right)-q_{\mu }{q}_{\nu }{\Pi }_2(q^2),$$

$${\hat{\Pi }}^{\text{charged}}(q^2)=\frac{1}{3q^2}\sum _{j=1}^3\left[{\Pi }_{jj}^{\text{charged}}(q_0,0)-{\Pi }_{jj}^{\text{charged}}(0)\right]={\Pi }_1(q_0^2)-{\Pi }_1(0).$$

$$M=\left(\begin{array}[]{cc}\frac{1}{\sqrt{2}}\pi^{0}&\pi^{+}\\ \pi^{-}&-\frac{1}{\sqrt{2}}\pi^{0}\end{array}\right),\;\;\;\;J_{\mu}=\left(\begin{array}[]{cc}\frac{2}{3}j_{\mu}&j_{+\,\mu}\\ j_{-\,\mu}&-\frac{1}{3}j_{\mu}\end{array}\right).$$

$$D_{\mu}M=\partial_{\mu}M-i\left[J_{\mu},M\right]=\left(\begin{array}[]{cc}\frac{1}{\sqrt{2}}\partial_{\mu}\pi^{0}+i\pi^{+}j_{-\,\mu}-i\pi^{-}j_{+\,\mu}&\partial_{\mu}\pi^{+}-i\pi^{+}j_{\mu}+i\sqrt{2}\pi^{0}j_{+\,\mu}\\ \partial_{\mu}\pi^{-}+i\pi^{-}j_{\mu}+i\sqrt{2}\pi^{0}j_{-\,\mu}&-\frac{1}{\sqrt{2}}\partial_{\mu}\pi^{0}-i\pi^{+}j_{-\,\mu}+i\pi^{-}j_{+\,\mu}\end{array}\right),$$

Full text

1. I Introduction

2. II The hadronic vacuum polarization

3. III Finite-size effects to the scalar vacuum polarization

4. III.1 The full finite-size effects

5. III.2 Charged currents

6. III.3 Universality of the finite-size corrections

7. IV Numerical validation

8. IV.1 Scalar QED on a lattice

9. IV.2 Lattice perturbation theory Monte-Carlo strategy

10. IV.3 Lattice scalar QED simulations

11. IV.4 Numerical results

12. V Conclusions

13. A In continuous infinite volume

14. A.1 Diagrams in Minkowski space

15. A.2 Master integrals

16. B The scalar loop integrals $\Omega_{\alpha,\beta}$

17. C Lattice scalar QED renormalisation scheme

18. D Explicit forms of energy-integrated diagrams

19. D.1 Diagram (S)

20. D.2 Diagram (T)

21. D.3 Diagram (C)

22. D.4 Diagram (E)

23. D.5 Diagram (X)

I Introduction

One of the most precisely measured quantities in particle physics is the anomalous magnetic moment of the muon $a_{\mu}=\frac{g_{\mu}-2}{2}$, where $g_{\mu}$ describes the ratio of couplings of the muon spin and orbital angular momentum to an external magnetic field. Historically, Dirac’s original tree-level prediction $g=2$ was in good agreement with experimental results, but the discrepancy which eventually arose became very strong evidence in support of quantum electrodynamics (QED). Both experimental measurement and Standard-Model predictions for $a_{\mu}$ have by now reached a precision of about 0.5ppm where a tension of 3.5-4$\sigma$ is observed (Bennett et al., 2006; Jegerlehner, 2018; Davier et al., 2017; Keshavarzi et al., 2018). Efforts are therefore under way to increase the accuracy of both measurements and theoretical predictions. To address the former two new experiments have been planned, E989 at Fermilab (Logashenko et al., 2015) and E34 at J-PARC (Otani, 2015). The Fermilab experiment is expected to lead to first new results in 2019 and the J-PARC experiment is expected to begin in 2020. Both experiments aim at increasing the experimental accuracy by a factor of 4 to 0.14ppm. To address the latter, we note that the main challenge on the theoretical side comes from non-perturbative contributions, namely the hadronic vacuum polarization (HVP) and hadronic light-by-light scattering (HLbL), of which the HVP constitutes the dominant contribution to the theoretical uncertainty. The traditional approach to estimating the HVP uses dispersion relations together with the optical theorem to relate it to the measured cross section of $e^{+}e^{-}$ to hadrons (Davier et al., 2017; Keshavarzi et al., 2018; Davier et al., 2011; Hagiwara et al., 2011). More recently, there has been a significant progress in calculation of the HVP from first principles using lattice QCD (Boyle et al., 2012; Della Morte et al., 2012; Burger et al., 2014; Chakraborty et al., 2014, 2015; Bali and Endrődi, 2015; Chakraborty et al., 2016; Blum et al., 2016a, b; Chakraborty et al., 2017; Della Morte et al., 2017; Borsanyi et al., 2018; Blum et al., 2018; Giusti et al., 2018a; Davies et al., 2019).

Based on simple power counting we expect strong and electromagnetic isospin breaking effects to contribute at the percent level. Given that this corresponds to the level of precision state-of-the-art lattice simulations are able to achieve, these effects need to be included in future calculations. Here we concentrate on electromagnetic effects which can be computed in the lattice-discretized finite-volume theory in several ways. Common to all approaches is the difficulty of defining charged states in a finite volume with periodic boundary conditions and the resulting singularities from photon zero-modes which need to be dealt with. In QEDTL (Duncan et al., 1996, 1997; Borsanyi et al., 2013; Ishikawa et al., 2012; Aoki et al., 2012; de Divitiis et al., 2013; Fodor et al., 2016a) the global zero mode is removed by hand while in QEDL (Hayakawa and Uno, 2008; Davoudi and Savage, 2014; Fodor et al., 2016b; Blum et al., 2007, 2010; Ishikawa et al., 2012; Giusti et al., 2018b, 2017a, 2017b; Blum et al., 2018; Boyle et al., 2017; Lee and Tiburzi, 2016; Matzelle and Tiburzi, 2017; Davoudi et al., 2018) the photon zero-mode is subtracted individually on every time slice. An alternative avenue is to perform simulations with a massive photon (Endres et al., 2016) followed by an extrapolation to zero photon mass to obtain physical results (Endres et al., 2016; Bussone et al., 2018). In yet another approach one introduces charge-conjugation boundary conditions (Polley, 1993; Wiese, 1992; Kronfeld and Wiese, 1993, 1991; Lucini et al., 2016; Hansen et al., 2018) which allow for constructing gauge-invariant charged states in a finite volume. QED corrections have been performed in (Blum et al., 2007, 2010; Borsanyi et al., 2013, 2015; Horsley et al., 2016a, b; Basak et al., 2016; Fodor et al., 2016a; Blum et al., 2018). Isospin-breaking corrections to the HVP have been explicitly considered in (Giusti et al., 2017b; Boyle et al., 2017; Chakraborty et al., 2018a, b; Blum et al., 2018).

A recurring systematic in QCD+QED calculations is the presence of large finite volume (FV) effects, which scale as $1/L^{n}$ with the box size $L$ for some exponent $n$. This is the result of the photon being a massless particle and the long-ranged nature of electromagnetic interactions. The finite-volume corrections have been studied in effective field theories for the meson masses in (Hayakawa and Uno, 2008; Borsanyi et al., 2015; Davoudi and Savage, 2014; Endres et al., 2016; Lucini et al., 2016; Davoudi et al., 2018) and decay rates in (Lubicz et al., 2017). The finite-volume correction to the HVP at order $\alpha$ however has not been previously calculated and this is the subject of this paper. We will extend the methodology for computing finite volume effects described in (Davoudi et al., 2018) to the electromagnetic correction to the HVP.

This paper is organized as follows. section II describes the preliminaries of the HVP function, which are relevant in both finite and infinite volumes. section III describes the analytic derivation of the finite-volume correction of the HVP with $\mathcal{O}(\alpha)$ electromagnetic correction. This is the main result of this paper. Finally, section IV contains numerical tests of the analytic expressions derived in section III. In appendix A we present a calculation of the NLO HVP in continuum Minkowski space.

II The hadronic vacuum polarization

The main object of interest is the Euclidean 2-point function

[TABLE]

where $j_{\mu}(x)$ is a charged or neutral vector current and $q^{2}$ is the external, Euclidean photon momentum. We start by presenting the calculation for neutral currents relevant for the HVP and then, since the calculation is equivalent up to numerical factors in the summation of diagrams, briefly present the result for charged ones. Note that for neutral currents Ward-Takahashi identities imply that $\Pi_{\mu\nu}(q)=\left(q_{\mu}q_{\nu}-q^{2}\delta_{\mu\nu}\right)\Pi(q^{2})$. The quantity $\Pi(q^{2})$ is ultraviolet divergent and it is conventional to calculate the finite, subtracted quantity

[TABLE]

This may be expanded in powers of the electric charge as $\hat{\Pi}\left(q^{2}\right)=\hat{\Pi}^{(0)}\left(q^{2}\right)+\hat{\Pi}^{(1)}\left(q^{2}\right)+\ldots$, where $\hat{\Pi}^{(0)}\left(q^{2}\right)$ and $\hat{\Pi}^{(1)}\left(q^{2}\right)$ are the leading order (LO) and next-to leading order (NLO) terms, i.e. $\mathcal{O}\left(1\right)$ and $\mathcal{O}\left(\alpha\right)$, respectively.

Although the two-pion contribution to the HVP is small compared to the vector resonance one, it is the lightest contribution and dominates finite-volume effects (Aubin et al., 2016). Here we consider the electromagnetic corrections to this contribution, and use scalar QED as an effective theory of elementary pions. The scalar QED Lagrangian in Euclidean space is given by

[TABLE]

for a scalar field $\phi$, a photon field $A_{\mu}$ and the electromagnetic tensor $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$. We only consider the leading order scalar interactions, higher-order $\mathcal{O}\left(\lambda\alpha^{2}\right)$ contributions, where $\lambda$ is the four-scalar vertex coupling, enter at three loop order. The connected diagrams needed at NLO for the HVP are therefore those in fig. 1. Seeing that some of these diagrams are equal up to relabelling of momenta in the loops, we only need to calculate seven topologies, namely (A), (B), (E1), (C1), (T1), (S) and (X). Diagrams (A) and (B) do not depend on the external momentum and thus cancel in the subtraction, so only (E), (C), (T), (S) and (X) contribute. The topology subscripts have here been suppressed and will remain so in the rest of the paper. The labelling refers to embedded sunrise (E), contact (C), embedded tadpole (T), sunset (S) and photon exchange (X). It should be noted that also the diagrams in (D) in fig. 2 are in general needed at NLO, but they are excluded for $\mathbf{q}=0$ in $\mathrm{QED}_{\mathrm{L}}$ and are in infinite volume simply related to the square of the LO contribution, given by diagrams (F) and (G), as will be discussed in more detail later. The specific choice of the photon rest frame is elaborated on in section III.

For completeness, we also calculated the NLO HVP in infinite volume. We considered both a Euclidean lattice using lattice perturbation theory, as well as continuum Minkowski space. The Minkowski space calculation is presented in appendix A.

Using the Feynman rules from the continuum Euclidean space Lagrangian yields the momentum-space integrands of the diagrams

[TABLE]

where $k$ is the photon loop momentum, $\ell$ is the pion-loop momentum and $d$ is the number of dimensions. Similar expressions for Minkowski space are given in section A.1.

III Finite-size effects to the scalar vacuum polarization

We can express the renormalized HVP function $\hat{\Pi}(q^{2})$ through the following trace of the subtracted vector two-point function

[TABLE]

where the photon rest frame has been specifically chosen. There are two main reasons for choosing this frame. First and foremost, this is typically the frame used in current lattice calculations. Moreover, it simplifies the finite-volume calculation immensely, in particular as the coefficients $c_{j}$ defined below then are independent of the photon momentum.

Note that diagrams (A) and (B) automatically vanish in the subtraction in eq. 13, as they are independent of the external momentum. Moreover, the disconnected contribution (D) is zero in $\text{QED}_{\mathrm{L}}$ because the photon propagator vanishes in the rest frame. We are thus left with diagrams (E), (C), (T), (S) and (X), so that, including all permutations of the diagrams, the $\mathcal{O}(\alpha)$ contribution to the HVP can be written as

[TABLE]

where $\Pi_{U}(q^{2})$ denotes the contribution from diagram (U). Next we define the corresponding integrand (excluding the factors of $2\pi$ in the measure) as $\pi_{U}\left(k,\ell,q_{0}\right)$.

In finite-volume, we assume space to be periodic with spatial extent $L$ and time to remain infinite. We now present the procedure followed to determine the finite-size effects to $\Pi^{(1)}(q^{2})$ which decay like powers of $1/L$. This strategy is a direct generalization of the procedure for one-loop integrals in Ref. (Davoudi et al., 2018). The remaining part of this section is a formal description of our approach to compute the large volume expansion. Although the final result presented in section III.1 is quite compact, intermediate expressions can be quite cumbersome. It is therefore desirable to implement the whole strategy in a computer algebra system. The calculations presented here were performed using FORM (Vermaseren, 2000) and Mathematica (Inc., ), and the associated Mathematica notebook is provided as a supplement of this paper under the General Public License version 3. Intermediate products of the derivation are provided for future reference in appendix D.

For a given diagram, we start by computing the two energy integrals in $k_{0}$ and $\ell_{0}$ using contour integration. Feynman integrands are rational functions and this integration is systematically feasible analytically. We thus obtain

[TABLE]

In analogy with eq. 13, we also define the subtracted quantities $\hat{\rho}_{U}$. The finite volume effects on $\hat{\Pi}(q^{2})$ for diagram (U) in $\text{QED}_{L}$ can then be written as

[TABLE]

where finite-volume sums are on quantized momenta of the form $\mathbf{k}=\frac{2\pi}{L}\mathbf{n}$ with $\mathbf{n}$ a vector with integer components, and a primed sum means that the origin is excluded, which here comes from the $\mathrm{QED}_{\mathrm{L}}$ prescription. One important aspect here is that we are only considering the $q^{2}>0$ case. This means that pions in diagrams are purely virtual and cannot generate power-like finite-volume effects through on-shell singularities. Using the Poisson summation formula for the pion part yields

[TABLE]

where the omitted terms denoted by ellipsis are the exponentially suppressed contributions from the virtual pions.

To determine power-like finite-size effects in the five diagrams (E), (C), (T), (S) and (X), we closely follow the strategy laid out in Davoudi et al. (2018). One starts by isolating the singularities in the photon momentum $\mathbf{k}$ in $\rho_{U}(\mathbf{k},\boldsymbol{\ell},q_{0})$,

[TABLE]

where $n_{U}$ is an integer that depends on the diagram in question, $\hat{\mathbf{k}}=\mathbf{k}/\left|\mathbf{k}\right|$, and $\overline{\rho}_{U}(\mathbf{k},\boldsymbol{\ell},q_{0})$ is an analytic function in the norm $\left|\mathbf{k}\right|$ such that $\overline{\rho}_{U}(\mathbf{0},\boldsymbol{\ell},q_{0})=0$. The analytical structure of all five diagrams is such that $n_{U}\leq 1$. If we now substitute $\mathbf{k}=\frac{2\pi}{L}\mathbf{n}$ and expand in $1/L$, the finite volume effects for diagram $U$ can be written as a power series in $1/L$ (up to exponentially small corrections),

[TABLE]

The coefficients $\xi_{j}^{U}\left(q_{0}^{2}\right)$ are given by

[TABLE]

where $\Delta_{\mathbf{n}}^{\prime}$ is, as in Ref. (Davoudi et al., 2018), the $\mathrm{QED}_{\mathrm{L}}$ sum-integral difference operator

[TABLE]

Although $\overline{\rho}_{U}(\mathbf{k},\boldsymbol{\ell},q_{0})$ is an analytic function in the norm $\left|\mathbf{k}\right|$, the norm itself in not analytic in the components of $\mathbf{k}$ at the origin, which generates the $\mathcal{O}(1/L^{4})$ effects in eq. 19. We will now present the full expressions for the finite-size effects to each of the five diagram topologies (S), (T), (C), (E) and (X).

III.1 The full finite-size effects

To find the finite-size effect to a diagram (U), the last step to perform is the calculation of the $\xi_{j}^{U}$ coefficients in eq. 20. For the specific kinematics chosen here, i.e., spatial momenta equal to zero (cf. eq. 13), the integrand $u_{j}\left(\hat{\mathbf{n}},\boldsymbol{\ell},q_{0}\right)$ is independent of the photon momentum direction $\hat{\mathbf{n}}$. The function $\xi_{j}^{U}(q_{0}^{2})$ then has the form

[TABLE]

where identically to Ref. (Davoudi et al., 2018), we define the coefficients $c_{j}=\Delta_{\mathbf{n}}^{\prime}|\mathbf{n}|^{-j}$. These can be calculated numerically in several ways, and one possibility is presented in (Davoudi et al., 2018). The first three coefficients are $c_{0}=-1$, $c_{1}=-2.83729748\ldots$ and $c_{2}=\pi c_{1}$.

The functions $\phi_{j}(q_{0}^{2})$ in eq. 22 can be written as linear combinations of integrals of the form

[TABLE]

where $z=q_{0}^{2}/m^{2}$, and

[TABLE]

They arise after integrating the angular dependence of the integrals over momentum $\boldsymbol{\ell}$ for which $\left|\boldsymbol{\ell}\right|=mx$. There are several useful recursion relations and properties for these integrals, which we summarise in appendix B. For instance, for $\alpha+2\beta>3$ in $d=4$, it is possible to write any $\Omega_{\alpha,\beta}$ as a linear combination of the six simple functions $\Omega_{2,1}$, $\Omega_{3,1}$, $\Omega_{4,1}$, $\Omega_{5,1}$, $\Omega_{0,2}$ and $\Omega_{1,2}$ as well as their respective derivatives. The complete list of expressions that lead to these integrals, in particular, the expansion eq. 18, are given explicitly for all diagram topologies in appendix D.

Finally, we summarize here the final expressions for the finite-volume effects to each diagram, where every $\Omega_{\alpha,\beta}$ term implicitly depends on $z=q_{0}^{2}/m^{2}$

[TABLE]

and where all the expressions are given up to $\mathcal{O}(\frac{1}{L^{4}},e^{-mL})$ corrections. We can sum these terms according to eq. 14. The resulting series in $1/L$ for the HVP at NLO is

[TABLE]

where one notices the important cancellation of the $1/L^{2}$ terms. This result can be understood from the underlying physics since the current is neutral and a photon far away thus sees no charge. This cancellation has potentially important consequences regarding the prediction of the contribution $a_{\mu}^{\mathrm{HVP}}$ from the HVP to $a_{\mu}$ using lattice simulations. Indeed, for typical physical simulations with $mL>4$, one has $1/(mL)^{3}<1.5\%$. Under the safe assumption that the QED corrections to $a_{\mu}^{\mathrm{HVP}}$ are $\mathcal{O}(1\%)$, the electromagnetic finite-size effects discussed here would represent a contribution smaller than $0.02\%$, well below the $0.1\%$ level required to reduce by a factor of $4$ the current theoretical uncertainties on $a_{\mu}^{\mathrm{HVP}}$. Finally, for $mL>4$ one has $e^{-mL}<1.8\%$, which means that in this regime the new, power-like finite-size corrections introduced by QED are in principle not dominant compared to the exponential QCD effects. In the following sections, we demonstrate that this cancellation does not occur for charged currents, and that it is universal in full QCD+QED and therefore directly applicable to lattice results.

III.2 Charged currents

For the neutral currents only charged pions are considered. If also $\pi^{0}$ is included, the current $j_{\mu}$ in eq. 1 can be charged and the current-current correlator can therefore be rewritten as

[TABLE]

for two functions $\Pi_{1}(q^{2})$ and $\Pi_{2}(q^{2})$ that are equal for neutral currents. We are again interested in the case $q_{\mu}=(q_{0},\mathbf{0})$ and calculate the subtracted quantity

[TABLE]

The function $\Pi_{1}(q^{2})$ can be expanded in the electromagnetic coupling just as before, and we here denote the NLO contribution by $\Pi_{1}^{(1)}(q^{2})$. The possible topologies of the NLO diagrams are the same also here, but having charged currents implies that some of them may be forbidden and the overall numerical factors can be different compared to the netural case for those that are not. In order to find these differences, we include the neutral pion by defining the meson matrix $M$ and the current matrix $J_{\mu}$ as

[TABLE]

The covariant derivative of $M$ can then be put in the form

[TABLE]

and the kinetic part of the Lagrangian is given by

[TABLE]

Using this, we find that diagram (X), one of the permutations of diagram (T), two of the permutations of (C) and one of the permutations of diagram (E) are forbidden. Moreover, the overall numerical factors change such that the relevant NLO contribution becomes

[TABLE]

This yields the NLO FV effects as

[TABLE]

up to $\mathcal{O}(\frac{1}{L^{4}},e^{-mL})$ corrections. The $1/L^{2}$ part does not vanish here. This is expected, since the current no longer is neutral and the physical argument used for the neutral case no longer applies.

III.3 Universality of the finite-size corrections

In the above we showed that in one-loop scalar QED the leading contribution to the FV effects on the subtracted vacuum polarization function $\hat{\Pi}(q^{2})$ is of order $1/L^{3}$ in $\textrm{QED}_{L}$. We will now show that this conclusion is independent of the effective-field-theory formulation chosen for the finite volume calculation.

The $\mathcal{O}(\alpha)$ corrections to the current-current correlator $\Pi_{\mu\nu}(q)$, which we denote $\Pi_{\mu\nu}^{(1)}(q)$, can be written as

[TABLE]

where the symbol $\int_{x,y,z}$ is an abbreviation for the integration over the $3$ space-time positions $x$, $y$, and $z$. This amplitude is identical to the amputated light-by-light scattering Green’s function with two legs contracted with the photon propagator. We will argue below that the light-by-light matrix element at our choice of kinematics is free of singularities and therefore expected to have only exponential finite volume corrections. As a consequence, the only source of power-like corrections will necessarily have to come from the photon propagator pole.

A useful form-factor decomposition of the light-by-light amplitude with form factors that are free of kinematic singularities is given in equation (3.14) of (Colangelo et al., 2015). We note that all tensor structures have at least one factor of each of the $q_{1}$, $q_{2}$, $q_{3}$ and $q_{4}$. In our formula, we need to replace two of those external momenta with $k$ and the remaining two with $q$ and $-q$, respectively. Therefore all the tensor structures will be proportional to either $k^{2}$ or $k_{\mu}k_{\nu}$ for some Lorentz indices $\mu$ and $\nu$, which in turn can be replaced with $k^{2}$ using the $d$–dimensional integral (or sum in FV) formula:

[TABLE]

This means that each of the tensor structures contribute a factor of $k^{2}$, which cancels the pole in the photon propagator. Since the light-by-light form factors $F(k^{2},q^{2},k\cdot q)$ are free of kinematic singularities, and we work with Euclidean momenta which cannot give rise to any on-shell singularities, they do not have any poles in $k$. This means that the leading FV correction will come from a term which is constant in $k$, which has the contribution proportional to $c_{0}/L^{3}=-1/L^{3}$, completing the proof.

The above argument can be simplified by considering the scalar HVP form factor

[TABLE]

The function $\Pi$ then has the form

[TABLE]

where the kernel $K_{\rho\sigma}$ satisfies the Ward identities $k_{\rho}K_{\rho\sigma}=k_{\sigma}K_{\rho\sigma}=0$. The kernel function $K_{\rho\sigma}$ is a momentum space amplitude and any singularities must correspond to physical states. In Euclidean space with external momenta being real, the function can not have any poles corresponding to physical states, which would have to satisfy $p^{2}+m^{2}=0$ where $p$ is the momentum going through any cut in the diagram. We can decompose $K_{\rho\sigma}$ in a series of tensor structures multiplying form factors, which can generally be written as

[TABLE]

Since $K_{\rho\sigma}$ is free of singularities and the tensor structures are linearly independent, the form factors must be free of singularities as well. The Ward identities $k_{\rho}K_{\rho\sigma}=k_{\sigma}K_{\rho\sigma}=0$ impose a relation on the form factors simplifying the expression to

[TABLE]

As noted before, the form factors $F_{kk}$ and $F_{kq}$ are free of singularities, however $F_{kq}$ must have a zero at $k\cdot q=0$ to cancel the pole originating from the tensor structure it is multiplying. We conclude that $F_{kq}$ must be proportional to $k\cdot q$. Finally, the form factor decomposition of $K$ consistent with Lorentz symmetry, parity, and Ward identities is

[TABLE]

As before, we note that the form factors $F_{1}=-F_{kq}/(k\cdot q)$ and $F_{2}=F_{kk}$ do not have poles in $k^{2}$ and the tensor structure has two factors of $k$ which become $k^{2}$ under the integral, which cancels the photon propagator pole. As before, this results in the leading contribution to the FV correction to be proportional to $c_{0}/L^{3}$.

IV Numerical validation

In this section we provide two different numerical checks of the finite-volume corrections derived in section III.1. Scalar QED is ideally suited for numerical simulations. Indeed, as we will now explain in detail, the theory can be written on a discrete space-time simply by replacing derivatives with finite differences. In the two following subsections, we describe two different Monte-Carlo strategies to compute the volume dependence of the scalar vacuum polarization. Firstly, the master formula eq. 17 is evaluated directly using a Monte-Carlo integrator. Secondly, the finite-volume vacuum polarization is calculated at $\mathcal{O}(\alpha)$ using lattice -scalar-QED simulations, following the strategy described in (Davoudi et al., 2018). Finally, we discuss the comparison of these results with analytical predictions.

IV.1 Scalar QED on a lattice

In this section we explain our definition of the lattice discretized theory. The conventions and notations are identical to (Davoudi et al., 2018). We consider space-time to be an Euclidean four-dimensional lattice with spatial extent $L$, time extent $T$, and lattice spacing $a$. Lattice QED is then defined by the action

[TABLE]

with scalar and gauge actions

[TABLE]

respectively, with $\Delta=m^{2}-\sum_{\mu}D^{*}_{\mu}D_{\mu}$. The summation is over all the sites of the lattice. The covariant derivative is defined in terms of the $U(1)$ gauge link $U_{\mu}(x)=e^{iqaA_{\mu}(x)}$, where $q$ is the electric charge of the scalar particle, as

[TABLE]

We also introduce the forward derivative $\delta_{\mu}A_{\mu}(x)=a^{-1}[A_{\mu}(x+a\hat{\mu})-A_{\mu}(x)]$, which appears in the Feynman gauge-fixing term. The electromagnetic tensor is defined as

[TABLE]

Expectation values in this theory are expressed in terms of the path integral

[TABLE]

where the integral measures represent integrations over the field variable at each lattice site. The subscript L indicates that we are working within the QEDL prescription where the spatial zero mode is set to zero on each time slice,

[TABLE]

Below we will expand the path integral to NLO in $\alpha$. To this end it is instructive to first integrate out the scalar fields analytically,

[TABLE]

where $O_{\rm Wick}$ represents the observable after the Wick contraction. The action is symmetric under $A_{\mu}\to-A_{\mu}$ and therefore, contributions odd in $q$ do not contribute to expectation values. To NLO we can therefore set ${\rm det}(\Delta)=1$.

We rewrite the operator $\Delta$ with the help of the translation operator $\tau_{\mu}f(x)=f(x+a\hat{\mu})$, as

[TABLE]

The expansion of $\Delta$ in the electric charge $q$ takes the form,

[TABLE]

where

[TABLE]

Inserting the kernel $\Delta$ expanded in $q$ into the scalar-QED action,

[TABLE]

allows us to identify the Feynman rules for the inverse free propagator, the scalar-photon-vertex and the scalar tadpole, respectively. In particular, the scalar propagator in the background field $A_{\mu}$ is then readily given by

[TABLE]

From this expansion it is a simple exercise to derive the associated Feynman rules for lattice perturbation theory.

IV.2 Lattice perturbation theory Monte-Carlo strategy

In order to numerically check the analytic results we numerically calculate the finite-size corrections $\Delta\hat{\Pi}_{U}$ for each diagram $U$ in scalar QED using lattice perturbation theory (LPT). We present below the analytic expressions for the diagrams (E), (C), (T), (S), and (X) in lattice perturbation theory, which are the discrete version of eqs. 8, 9, 10, 11 and 12.

[TABLE]

where $\hat{k}_{\mu}=\frac{2}{a}\sin(\frac{ak_{\mu}}{2})$.

On the lattice, there is potentially an infinity of new scalar-photon vertices because of the compactification of the gauge field in eq. 48. These vertices are classically discretisation effects, but at the quantum level they can generate finite contributions when multiplying power divergences, and ignoring them can potentially break Ward-Takahashi identities. If one consistently keeps contributions which do not vanish in the continuum limit, four new diagrams appear in lattice perturbation theory, represented in fig. 3. Both diagrams (L3) and (L4) are independent of the external momentum and therefore vanish in the subtracted vacuum polarization function eq. 13.

The integrand for diagrams (L1) and (L2) is given by

[TABLE]

We integrate these expressions using the VEGAS algorithm (Lepage, 1978), and more specifically its C++ implementation in the Cuba library (Hahn, 2005). This integration algorithm builds upon Monte Carlo techniques and creates histograms approximating the shape of the function which are then used as probability distributions for importance sampling. This is particularly useful for the integrals considered here, which are eight-dimensional, and have a complicated sawtooth-like structure, as we discuss now. In finite volume, the lattice momentum is discretized and the corresponding sums can be dealt with in VEGAS by realising that for a function $f(k)$

[TABLE]

where $\lfloor k\rfloor$ is the floor operator rounding $k$ down to the nearest integer. This is extendable to any number of dimensions. As for the analytic results, we assume an infinite time extent and pions are in infinite volume, cf. eq. 17, so that only three of the eight integrals are sums in the finite-volume calculation. The implementation of the calculation is distributed as a C++ source code under the General Public License v3 in the supplementary material of this paper, and it features an option to also have the pions in a finite volume as instructed in the comments. This is particularly useful when comparing to lattice data, as discussed later.

Each diagram depends on the lattice spacing $a$. We introduce a scaling parameter $\sigma$ such that the lattice spacing is varied according to $a\rightarrow a/\sigma$ and calculate the diagrams for four different values of $\sigma$, namely $\sigma\in\left\{1,1.5,2,3\right\}$, from which a continuum extrapolation is made by fitting against some polynomial in $a$. We find, using a pion mass such that $am=0.2$, that the $a$ dependence is mild. We find the best description of the data in terms of a leading $\mathcal{O}(a^{2})$ correction, as expected from the $\mathbb{Z}_{2}$ symmetry of scalar QED.

Rewriting the sums as in eq. 64 yields sawtooth-like behavior since the integrands of the finite-size effects then are of the form $f(\mathbf{k})-f(\lfloor\mathbf{k}\rfloor)$. The number of discontinuities in this function is on the order of $(\sigma L/a)^{3}$ (or $(\sigma L/a)^{6}$ if also the pions are put in finite volume) which means that it can be hard to sample the integrand efficiently and thus get reliable values and errors from Cuba. The reliability can be checked by comparing the results from calculations with a varying number of Monte Carlo evaluation points for a certain $\sigma$. We find that using $10^{11}$ points gives reliable results for $\sigma<4$. Our result are summarized at the end of this section.

IV.3 Lattice scalar QED simulations

An alternative avenue which we also explore is to evalute the lattice-discretized path integral in eq. 52 by means of a Monte Carlo integration for a series of different spatial extents $L$. This allows for mapping out the volume dependence, thereby checking our analytical predictions. Instead of numerically solving the momentum sums as in the previous section one directly samples the path integral in eq. 52. In particular, we compute the vacuum polarization tensor $\Pi_{\mu\nu}(q)$ as the discrete Fourier transform of the two-point function

[TABLE]

with the lattice conserved vector current

[TABLE]

After carrying out the Wick contractions we can write the expression for the vacuum polarization tensor in terms of the propagator in eq. 57 acting on a point source $\delta(x)$,

[TABLE]

where the expectation value represents the functional integration on the gauge potential $A_{\mu}$. We evaluate this correlation function numerically using a setup identical to the one in (Davoudi et al., 2018), in fact the data used here are a side-product of the calculation presented in this previous work. The covariant Klein-Gordon equation is solved in a stochastic background field $A_{\mu}$ to form the interacting scalar propagator $\Delta^{-1}\delta(x)$. Using the expansion eq. 57, this can be achieved using the fast Fourier transform algorithm. As a consequence, this method has a reasonable numerical cost, which is independent from the chosen scalar mass and has a quasilinear complexity in the number of lattice points. We refer the reader to (Davoudi et al., 2018) for more details on the computational aspects.

In principle, the full $\mathcal{O}(\alpha)$ correction to the scalar vacuum polarization also receives contributions from the diagrams in fig. 4, coming from the 1-loop counter-terms of scalar QED. We assume that these counterterms are determined through a set of renormalization conditions in infinite volume, and therefore are independent of the volume. Because these diagrams do not contain photon propagators, they clearly do not contribute to eq. 17. However, the same formula assumes scalar particles to be in infinite-volume, which is not the case in the lattice simulation. Although these finite-volume corrections are exponentially suppressed, they can be greatly enhanced by the ultraviolet-divergent values of the counter-terms. We therefore included these diagrams to ensure that exponential finite-volume corrections are negligible for reasonably large values of $mL$ (the typical threshold for lattice QCD simulations is $M_{\pi}L>4$). The details of the renormalisation prescription used here are given in appendix C. The cost of computing the extra counter-term diagrams is negligible, since they do not depend on the gauge field.

IV.4 Numerical results

In fig. 5 we compare the analytic results to LPT and lattice data. We use $am=0.2$ and $aq_{0}=8\pi/128$, i.e. $z=q_{0}^{2}/m^{2}\approx 0.964$. The red dashed line is the $1/L^{2}$ term and the green solid line is the full expression of the form $1/L^{2}+1/L^{3}$ in the corresponding analytic expression in section III.1. The purple points are the infinite volume pion LPT points for a finite $a$, and the crossed blue points are the continuum extrapolated values. The orange box shaped points are finite volume pion LPT data. From the infinite volume pion LPT data we clearly see that the full analytic form is much better than when including only the $1/L^{2}$ term and the agreement is excellent up to $1/mL<0.3$ for all diagrams. Other values of $z$ yield a similar level of agreement. We see that the lattice data starts to deviate from the analytic curve after $1/mL>0.2$, but the finite volume pion calculation reproduces precisely this behavior. We thus attribute the discrepancies to the exponential finite-size effects that are neglected in eq. 17 for the analytic calculation as well as in the infinite volume pion LPT Monte-Carlo. Moreover, for $mL\simeq 4$ we found that the difference between the infinite-volume pion and finite-volume pion data is on the order of $10^{-6}$, an order of magnitude smaller than the naive suppression from a factor of $\alpha$ between the LO and the NLO HVP, viz. $\hat{\Pi}^{(1)}\sim\alpha\hat{\Pi}^{(0)}\sim\alpha 10^{-3}\sim 10^{-5}$.

V Conclusions

We have performed a 2-loop calculation of the $\mathcal{O}(\alpha)$ corrections to the hadronic vacuum polarization in scalar QED. We presented the infinite volume results in terms of 2-loop master integrals from which we obtained an analytic expression for the finite volume correction to the HVP at this order. We found that even though each of the individual diagrams contributes as $1/L^{2}$, these terms all cancel when combined. We argued that this cancellation is expected on physical grounds for neutral currents and show that it does not occur for charged currents. We also argued that this cancellation is universal, i.e., it occurs regardless of the effective theory used to derive this result.

All our results were tested numerically using two different approaches - direct integration of lattice perturbation theory integrals using VEGAS and lattice scalar U(1) gauge theory with stochastically generated photon fields. We find good agreement between analytic results and results from both numerical approaches. While absent from our analytical expressions, exponentially suppressed finite volume effects are visible in our results from the lattice simulation

Finally, an important consequence of this work is that for the foreseeable future finite-volume effects on the QED corrections to the hadronic vacuum polarization are likely to be negligible in lattice simulations. For instance, we expect this to hold even if lattice computations aimed at matching experimental projections of a four-fold reduction in the error on $a_{\mu}$ by Fermilab (Logashenko et al., 2015) and J-PARC (Otani, 2015) down to 0.14ppm. This assumes a typical lattice simulation where the pion mass times the spatial extent is larger than four, for which this work estimates the finite-size effects to be at the level of only a few percent of the $\operatorname{O}(\alpha)$ correction to the HVP function $\hat{\Pi}(q^{2})$. Unless these effects come with an unnaturally large coefficient in the full theory, they should be negligible compared to the per-mil accuracy required on the HVP contribution to $a_{\mu}$. Of course, large coefficients cannot be excluded considering how critical it is to properly estimate the theoretical uncertainty on $a_{\mu}$, particularly in the perspective of confirming or excluding the current discrepancy between experiment and theory on this quantity. The results of this work together with simulations of full lattice QCD+QED even with a limited number of volumes, should allow to constrain the size of these effects.

Acknowledgements.

A.P. would like to thank the High Energy Physics Group of Lund University for its warm welcome, important parts of this work were initiated during A.P.’s visit in Lund. Lattice computations presented in this work have been performed on DiRAC equipment which is part of the UK National E-Infrastructure, and on the IRIDIS High Performance Computing Facility at the University of Southampton. T.J. and A.P. are supported in part by UK STFC grants ST/L000458/1 and ST/P000630/1. A.P. also received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No 757646. J.B. and N.H.T. are supported in part by the Swedish Research Council grants contract numbers 2015-04089 and 2016-05996, and by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No 668679. J.H. was supported by the EPSRC Centre for Doctoral Training in Next Generation Computational Modelling grant EP/L015382/1. A.J. received funding from STFC consolidated grant ST/P000711/1 and from the European Research Council under the European Union’s Seventh Framework Program (FP7/2007- 2013) / ERC Grant agreement 279757.

Appendix A In continuous infinite volume

In this section the HVP is considered in continuous infinite volume Minkowski space. We calculate $\Pi^{(0)}\left(q^{2}\right)$ and $\Pi^{(1)}\left(q^{2}\right)$ in $\overline{MS}$ and numerically compare their respective sizes. The corresponding calculation in QED can be found in (Kallen and Sabry, 1955; Barbieri and Remiddi, 1973).

In Minkowski space the scalar QED Lagrangian is

[TABLE]

The relevant counterterms for the parameters above, are defined in the counterterm Lagrangian

[TABLE]

In $d=4-2\varepsilon$ dimensions these are given by

[TABLE]

Note that diagrams (A) and (T) identically vanish in dimensional regularization, so that we are left with diagrams (F), (G), (D), (B), (E), (C), (S) and (X). The two HVP contributions can thus be written (cf. the FV case in eq. 14)

[TABLE]

Using the tensor structure of $\Pi^{\mu\nu}\left(q^{2}\right)$ and the Ward identity it is easy to see that the disconnected part is given by the squared LO contribution, $(D)=\left(\Pi^{(0)}\left(q^{2}\right)\right)^{2}$. The diagrams are given in section A.1.

Using Lorentz invariance identities and integration by parts, the 2-loop integrals can be rewritten in a basis of master integrals. The program Reduze2 (von Manteuffel and Studerus, 2012) employs a Laporta algorithm in order to do this, and allows the user to define such a basis. The master integrals used here are the $\overline{MS}$ subtracted parts of

[TABLE]

All but integral $M$ are divergent and thus require expansion in $\varepsilon$ in order to isolate the divergent parts from the finite ones, something which can be done in both Euclidean and Minkowski space. For a Euclidean spacetime the analytic results can be found in (Martin, 2003). However, working in Minkowski space, the corresponding expressions are here given in section A.2.

Note that in $\overline{MS}$ the threshold shift can, and does, induce a sign change of the imaginary part of $\Pi^{(1)}\left(q^{2}\right)$ at some $q^{2}$. However, this does not occur for an on-shell scheme or with the physical mass. The physical mass $m_{ph}^{2}$ is related to $m^{2}$ through

[TABLE]

The HVP can therefore also be expanded around this mass,

[TABLE]

where in the last step the disconnected part was separated from the NLO contribution. To simplify the expressions, let us further define

[TABLE]

The HVP contributions at LO and NLO are thus

[TABLE]

where the quantities with bars are the finite parts of the integrals in appendix A. These contributions as well as the corresponding subtracted quantities are plotted in fig. 6 for $m_{ph}=139.5$ MeV, $\mu=500$ MeV and $e=0.303$. As can be seen, the NLO parts are roughly two orders of magnitude smaller than LO, this is due to the additional power of $\alpha\sim 10^{-2}$. Moreover, it can be noted that $\Pi^{\text{disc}}$ on average is significantly smaller than the other parts, and that $\Pi^{\text{disc}}$ and $\Pi^{\text{2-loop}}$ combine to give the proper non-singular threshold behavior.

A.1 Diagrams in Minkowski space

Using the Feynman rules for the Lagrangian in eq. 68, the diagrams are

[TABLE]

A.2 Master integrals

Below, each Minkowski-space master integral in appendix A has been separated into a finite and an infinite part, the finite one denoted by a bar. The analytic expressions for these finite integrals are given in (Martin, 2003), and are in terms of Riemann zeta functions as well as the polylogarithm functions $\textrm{Li}_{2}$ and $\textrm{Li}_{3}$,

[TABLE]

Appendix B The scalar loop integrals $\Omega_{\alpha,\beta}$

Consider the dimensionless function $\Omega_{\alpha,\beta}$ given by

[TABLE]

where

[TABLE]

This function converges if and only if $\alpha+2\beta>3$. The relations between $m$, the external momentum $q^{2}$, the integration variable $x$ and the variable $z$ are $z=q^{2}/m^{2}$ and $x=\sqrt{\boldsymbol{\ell}^{2}}/m=\ell/m$. It is also possible to write $\Omega_{\alpha,\beta}$ explicitly in terms of hypergeometric functions as

[TABLE]

where $\,{}_{2}F_{1}$ is the hypergeometric function defined by

[TABLE]

However, the form eq. 80 is complex and might not be the most useful in practice. The $\Omega_{\alpha,\beta}$ functions are actually related to each other and can be expressed as combinations of a smaller set of functions. One starts by noticing the relations

[TABLE]

The identity eq. 82 directly implies that

[TABLE]

i.e. the index $\beta$ is decreased by taking derivatives in $z$. One can further note that

[TABLE]

which is easily proven by multiplying and dividing the integrand a factor $z+4(x^{2}+1)$. For the case $\alpha+2\beta>3$ it is possible to find a recursion relation by using eq. 83, this by partially integrating the definition of $\Omega_{\alpha,\beta}$:

[TABLE]

where in the last step the convergence requirement $\alpha+2\beta>3$ as well as eq. 83 were used. By writing x^{4}=x^{2}\Big{(}x^{2}+1\Big{)}-x^{2} one then finds the relation

[TABLE]

or, for $\alpha>0$,

[TABLE]

Inspired by the recursion relation in eq. 84, define for $\alpha>1$ the functions

[TABLE]

Now, using eq. 84 for $\alpha>1$ and an integer $\beta\geq 1$

[TABLE]

The second identity eq. 83, combined with an integration by parts of eqs. 79 and 84 leads to

[TABLE]

where we used the prime notation for derivatives in $z$. Using this last relation and eq. 89, it is clear that for any positive integer couple $(\alpha,\beta)$ such that $\alpha+2\beta>3$, $\Omega_{\alpha\beta}(z)$ is a linear combination of the six following functions and their derivatives

[TABLE]

Appendix C Lattice scalar QED renormalisation scheme

We start by rewriting the Lagrangian of scalar QED in term of the renormalized fields and parameters defined by $\phi_{0}=\sqrt{Z_{\phi}}\phi$, $A^{\mu}_{0}=\sqrt{Z_{A}}A^{\mu}_{0}$, $m=Z_{m}m$, $e_{0}=Z_{e}e$, where a subscript 0 denotes a bare quantity. The counterterm part of the lattice Lagrangian is given by

[TABLE]

At the order $\mathcal{O}(q^{2})$ relevant here, the electric charge $q$ does not renormalize, i.e. $Z_{A}=1$. The discretized action is gauge invariant and, as it is well known in the continuum, the theory can be renormalized by removing divergences in the self-energy function and by using $\delta_{V}=\delta_{Z}$ as imposed by the Ward-Takahashi identities. By denoting $\Sigma(p)$ the self-energy function at momentum $p$, we choose the following renormalization prescription

[TABLE]

with $q_{T}=(\frac{2\pi}{T},\mathbf{0})$ where $T$ is the time extent of the lattice. This prescription allows to compute the wave function renormalization at finite time extent. In all the finite-time numerical results presented in this paper we used $T=128a$. For $T\to\infty$, this prescription gives back the more traditional conditions, where one assumes that the self-energy and its derivative vanishes at $p^{2}=0$. For $T=128a$ and $am=0.2$ we found

[TABLE]

Appendix D Explicit forms of energy-integrated diagrams

The subtracted functions $\hat{\rho}_{U}$ can be written in the form

[TABLE]

where $C^{U}$, $A^{U}_{i}$ and $a_{ij}^{U}$ are functions of $\mathbf{k}$, $\boldsymbol{\ell}$ and $q_{0}$. The above factorization is chosen such that the dependence on $\mathbf{k}$ in these functions is different from pure powers of $\left|\mathbf{k}\right|$. This means that they can depend on $\mathbf{k}$ in denominators through the energy $\omega_{p}=\sqrt{\mathbf{p}^{2}+m^{2}}$ (which often shows up in denominators) as well as the combination $\mathbf{v}(\boldsymbol{\ell})\cdot\hat{\mathbf{k}}$ for the velocity $\mathbf{v}(\boldsymbol{\ell})=\frac{{\boldsymbol{\ell}}}{\omega_{\ell}}$ and the unit vector $\hat{\mathbf{k}}=\frac{\mathbf{k}}{|\mathbf{k}|}$. This separation is useful since, for a given $j$, a large volume expansion of $C^{U}A_{i}^{U}a_{ij}^{U}$, which multiplies $\left|\mathbf{k}\right|^{j-1}$, has leading power behavior of order $\left|\mathbf{k}\right|^{j-1}$. It is therefore only the $j=0$ term in the sum over $j$ which can give a contribution to $b_{1}^{U}$ and thus a $1/L^{2}$ finite-size correction, where the coefficients $b_{1}^{U}$ and $b_{0}^{U}$ are defined through

[TABLE]

Defining the velocity is particularly useful as any term with such a factor vanishes when integrating over $\mathbf{k}$. The velocities can enter also in the small $\left|\mathbf{k}\right|$ expansion, for instance through

[TABLE]

In this appendix, we list the non-vanishing functions $C^{U}$, $A^{U}_{i}$, $a_{ij}^{U}$ and $b_{i}^{U}$ separately for each diagram (U).

D.1 Diagram (S)

First consider (S), whose integrand for the momentum assignment in fig. 7 is

[TABLE]

The non-vanishing functions entering $\hat{\rho}_{S}$ and $\hat{\rho}_{S}^{\textrm{exp}}$ are

[TABLE]

D.2 Diagram (T)

Now consider the calculation of diagram (T) with momenta as in fig. 8. The integrand is

[TABLE]

The non-vanishing functions here are

[TABLE]

Note that $b_{1}^{T}=C^{T}$ since $C^{T}$ cannot be expanded in small $\left|\mathbf{k}\right|$. Also, since $b_{0}^{T}=0$ we cannot have any contributions of order $1/L^{3}$.

D.3 Diagram (C)

For diagram (C) with momenta as in fig. 9, the integrand is

[TABLE]

This gives

[TABLE]

D.4 Diagram (E)

The integrand for diagram (E), with the momentum assignment in fig. 10, is

[TABLE]

Here we have

[TABLE]

D.5 Diagram (X)

Assigning momenta as in fig. 11, the integrand of diagram (X) is

[TABLE]

The non-vanishing functions are now

[TABLE]