Stochastic computation of $g-2$ in QED

We perform a numerical computation of the anomalous magnetic moment ($g-2$) of the electron in QED by using the stochastic perturbation theory. Formulating QED on the lattice, we develop a method to calculate the coefficients of the perturbative series of $g-2$ without the use of the Feynman diagrams. We demonstrate the feasibility of the method by performing a computation up to the $\alpha^3$ order and compare with the known results. This program provides us with a totally independent check of the results obtained by the Feynman diagrams and will be useful for the estimations of not-yet-calculated higher order values. This work provides an example of the application of the numerical stochastic perturbation theory to physical quantities, for which the external states have to be taken on-shell.


Introduction
Feynman diagrams depict correlation functions in quantum field theories as an expansion in terms of small coupling constants in the Lagrangian. They provide us with the prime tool to calculate physical quantities once the expansion parameters are identified. Indeed, the most successful prediction of the quantum field theory, the anomalous magnetic moment (g − 2) in QED, has been calculated up to O(α 5 ) by the method of the Feynman diagrams (see Ref. [1] for a review), and the value is confirmed at the experiments with precisions of O(10 −12 ) [2]! The computation based on Feynman diagrams, however, becomes pretty much complicated at higher orders. Already at O(α 2 ), it took about a decade to converge on the correct value [3,4] since the famous leading-order formula, g −2 = α/π, was obtained by Schwinger in 1948 [5]. Even today, only experts can go beyond two-loop. (See Ref. [6] for the analytic expression at the three-loop level, which is obtained after a number of numerical works [7][8][9][10][11][12].) Recently, Aoyama, Hayakawa, Kinoshita, and Nio have achieved the calculation up to the five-loop level, which requires evaluations of more than ten thousand diagrams [1,[13][14][15][16][17][18][19][20][21][22][23], for which a numerical method has been used for the integration over Feynman parameters. The current numerical uncertainties of the highest order coefficient are of order a few percent. An independent numerical calculation has been performed for a subset of diagrams in Ref. [24].
There is a discrepancy between two groups beyond the statistical uncertainties although it is not a very serious problem in the current experimental precisions. It would not be feasible to go beyond five loops in near future due to the explosion of the number of diagrams and also to the increase of the dimension of the integrals. It is desired to have another method, which does not use the Feynman diagrams, for higher order calculations as well as for the check of the results up to five loops. Not only limited to the g − 2 computation, new methods to calculate the perturbative series will be quite welcome in general.
Correlation functions in the stochastic quantization are evaluated as long-time averages of fictitious time evolutions governed by the Langevin equation, and they are formally equivalent to the results of the path integral quantization. In the Langevin equation, one can expand the fields in terms of a coupling constant, λ, such as φ = φ (0) + λφ (1) + · · · , and write down the evolution equations for each φ (i) . The correlation functions at a fixed order O(λ n ) can be calculated as a sum of time averages, such as n i=0 φ (i) φ (n−i) in the case of two-point functions. There is no need of the Feynman diagrams in this method. All the diagrams are automatically summed in solving the set of the Langevin equations. If one tries to obtain analytic results, however, we simply recover calculations similar to the Feynman rules [30].
The formalism is more suited for a numerical simulation. By putting the theory on the lattice to make the degrees of freedom finite, the direct integration of the Langevin equation becomes possible. This provides us with an independent numerical method for the evaluation of the coefficients of perturbation series. Simplicity of the method is especially beneficial for higher-order computations.
There have been recent works on the applications of the stochastic perturbation theory in significantly reduces the cost of the fermion part [27,28]. This makes it possible to evaluate the coefficients at very high orders in QCD with dynamical fermions. The behavior of the expectation value at high orders has been used to extract the renormalon contributions [32,33].
In this paper, we perform a QED computation of the electron g − 2 in the numerical stochastic perturbation theory. Since the g − 2 factor in perturbative QED is one of the simplest physical quantities (and thus scheme independent) in particle physics, it is a good testing ground for the formalism. Also, if it is successful, one can provide a non-trivial check of the diagrammatic calculations as well as the predictions of higher order coefficients. We construct the action of lattice QED suitable for the stochastic calculation, and perform the g−2 evaluation up to O(α 3 ) on small lattices. The results are encouraging and we expect that simulations in currently available large scale computers and/or more optimized simulation codes can confirm the O(α 5 ) results of the diagrammatic methods.
Matching to the continuum theory requires various extrapolations such as the large volume limit, the continuum limit, as well as the limit to the on-shell fermion momentum. Among them, the on-shell extrapolation needs some indirect method as we work in the Euclidean space. In this work, we scan the fermion energies in the Euclidean region, and infer on-shell values by using analytic continuation. We discuss how such extrapolations cause uncertainties in the estimates of g − 2. We also discuss infrared (IR) problems in QED in a finite volume.

Lattice QED action
We first define the action of lattice QED which is to be used in the simulation. Unlike the non-perturbative lattice simulation, the fermion mass, m f , is the only dimensionful parameter in question as the Landau pole does not show up in the perturbative fixed-order calculations.
The continuum limit is, therefore, obtained by taking the dimensionless combination m f a to be zero while fixing the pole mass, m f , to be finite, where a is the lattice spacing. For this purpose, it is desirable to maintain the chiral symmetry on the lattice so that the fermion mass parameter m in the Lagrangian is not additively renormalized. If m ∝ m f , one can simply take ma → 0 as the continuum limit.
In the practical calculation, we need an IR regulator to eliminate the zero mode of the gauge field to avoid its numerically unstable random-walk evolution under the Langevin equation. We also need a gauge fixing term to eliminate the gauge degrees of freedom.
Although the ultraviolet (UV) part is already regulated by the lattice, we introduce a UV cut-off to avoid large logarithmic corrections, which makes the continuum limit far away.
Considering those requirements, we define the action in the lattice unit (a = 1) as follows: where and The derivatives are defined by and ∇ 2 = ∇ * µ ∇ µ . The Dirac operator D is that of the naive fermion, with m the fermion mass parameter. The chiral symmetry is maintained on the lattice because we use the naive Dirac operator, with which the fermion propagator has sixteen poles. Those doublers are taken care by the factor of 1/16 in Eq. (5). By this treatment, the action describes the system with a single fermion at least for perturbative calculations of quantities which do not involve γ 5 in the internal line. For the external fermion lines, one can choose the external momenta such that the doubler poles are not important.
The terms S g and S f are invariant under the gauge transformation: The gauge fixing term S gf eliminates the gauge degrees of freedom [34]. We take the noncompact QED action S g so that there is no self interactions of photons. This simplifies the perturbative expansion of the Langevin equation compared to the case of QCD. The QED coupling constant e only appears in Eq. (7).
A UV regulator is introduced through the factor e −∇ 2 /Λ 2 UV in the gauge action S g + S gf + S mass . This factor reduces the strength of the gauge interaction at high photon virtualityk 2 in a gauge invariant way. By an appropriate choice of Λ UV , one can reduce large logarithmic corrections which appear in the renormalization factors. In the continuum limit, we should send the UV cut-off to infinity while the physical fermion mass m f ∝ m fixed. This means we need to take the dimensionless quantity ma sufficiently small, while Λ UV a is fixed. Therefore, in the actual simulation, we can simply use a fixed value of Λ UV in the a = 1 unit. In the numerical calculation, this UV regularization is quite cheaply implemented by the use of the fast Fourier transform: going to the momentum space, multiply ek 2 /Λ 2 UV and coming back to the position space.
We introduce the photon mass m γ as the IR regulator as we usually do in the perturbative calculations in the on-shell scheme. We should take the limit of m γ /m → 0 later. The gauge invariance is broken by this procedure, but there are modified Ward-Takahashi identities that restrict the form of the quantum corrections. For detail, see Appendix A. Since the zero mode is regulated by the mass term, one can simply take the periodic boundary conditions for the gauge field.
In the literature, different IR regularizations in the lattice perturbation theories have been discussed. One is to remove the zero mode from the summation of the momentum [35,36].
Although this would be the simplest implementation, the estimation of the finite volume effects is quite non-trivial after the removal of the zero mode as discussed in Ref. [37]. When the zero mode is just removed, the theory is no longer a local field theory. Since we will heavily use the analyticity of the correlation functions in later discussions, non-locality of the theory may introduce an additional complication to the argument.
Imposing the twisted boundary condition is another way to remove the zero mode from the gauge field in the case of SU (N ) gauge theories [38]. In U (1) gauge theories, a similar regularization may be possible by introducing a background magnetic flux on one of the two tori, T 2 dA = 2π/e. However, since we use the non-compact QED action, where A µ (n) is treated as the fundamental field, there cannot be a magnetic monopole background to realize the flux.
For these reasons, we employ the simple IR regularization by the photon mass term. The effects due to the mass term are studied in Refs. [39,40], and it is found that the zero mode does modify the analytic structure, but in a controlled way. We will discuss it in Section 9.

Langevin equation
In the stochastic quantization, one can obtain the expectation values of operator products as a "time" average of the field products. The gauge field A µ (n) is promoted to obtain another "time" coordinate τ , A µ (n, τ ), and the evolution in the τ direction is governed by the Langevin equation: The Gaussian noise η µ (n, τ ) satisfies: where · · · η denotes an ensemble average over the random noise. The ensemble average of field products converges to the correlation functions obtained by the path integral quantization: as τ → ∞. This can be understood by the fact that the probability distribution of the path integral, e −S , is the fixed point of the Fokker-Planck equation derived from the Langevin equation in Eq. (9). The ensemble average in the left-hand-side can be evaluated as the "time" average of the Langevin trajectories of the field product, i.e., Therefore, by following long enough Langevin trajectories of fields, one can calculate the correlation functions. For a numerical method to integrate the Langevin equation, see Appendix B.
For our QED action, the only non-trivial part in the right-hand-side of Eq. (9) is the fermion part: The inverse of the Dirac operator, D −1 , is evaluated perturbatively as we see later. The trace of the space-time and spinor indices are taken by using a noise spinor ζ as: where the noise spinor satisfies The indices α and β denote those of the spinor. In practice, we generate only one noise spinor ζ for each discretized time step of the numerical Langevin simulation. The average is, therefore, taken during the large time evolution.

Perturbative expansion
We decompose the Langevin equation into those for each order in the perturbative expansion.
We follow the formalism of Refs. [27,29] for the treatment of the fermions.
We expand the field A µ (n, τ ) as and solve the Langevin equation for each A (p) µ . The equation for gauge fields in each order p is given by where | (p) denotes the p-th order term. The noise is applied only for the lowest order, A µ . Higher order fields get fluctuated through interactions with lower order fields.
The expansion of the force term, δS/δA µ , is straightforward except for the fermion part, for which The vertex factor δD/δA µ is evaluated as where The inverse of the Dirac operator, D −1 , can be obtained recursively [27] as Here, the propagator at the lowest order, D −1 0 , is given explicitly by where The momentum integration should be understood as dp 2π with p = 2πk/L for the periodic boundary conditions.
The propagator D −1 0 is diagonal in the momentum space, while D is local in the position space. Therefore, the multiplication of the Dirac operators in Eq. (22) can be effectively performed by the use of the fast Fourier transform in the numerical calculation. This significantly reduces the computational cost of the fermion implementation in the stochastic perturbation theory.
In the calculation, all the fields as well as their correlation functions are expanded in terms of the coupling constant e as in Eq. (16). The operations of the expanded quantities, such as the multiplication and the exponentiation are performed perturbatively as where O p is defined by repeatedly using Eq. (26). General functions such as the inverse and the square root of some operator are defined by the Taylor expansion in e, where the lowest order value, O (0) , is formally assumed to be non-zero.

Correlation functions
For the calculation of g − 2, we first compute the following set of correlation functions by using the Langevin evolutions. In the following discussion, all the correlation functions are implicitly expanded in terms of e. The parameters in the Lagrangian, ma, m γ a, Λ UV a, and ξ on the other hand, are kept as constants of O(e 0 ). In practice, we store the gauge field and the correlation functions as arrays of perturbative coefficients, and their operations such as multiplication are performed by applying the rule in Eq. (26).

Photon two-point function
The two-point function of the gauge field in the momentum space has the form wherek µ = 2 sin(k µ /2), andk 2 =k µkµ . The projection operators, g µν −k µkν /k 2 andk µkν /k 2 , are ambiguous atk 2 = 0, but it is not important in the following discussion, because we only use the two-point function at finitek 2 . The Fourier modes are defined as where x n +μ/2 stands for the mid point of the link specified by the site x n and the direction µ. The vacuum polarization, Π(k 2 ), is vanishing at the tree level, and is implicitly expanded as the perturbative series through the expansion of the two-point function in Eq. (28). The longitudinal part has no quantum correction as it is guaranteed by the Ward-Takahashi identity (see Appendix A): This fact is useful for the check of the gauge invariance.
We define the renormalization factor Z 3 as By this Z 3 factor, one can obtain the physical coupling constant, e P , as a power series of e as follows: where the limit of m 2 γ → 0 should be taken beforek 2 → 0. This coupling constant, e P , characterizes the coefficient of the long-range Coulomb force, and thus it is identified as the fine-structure constant in the quantum mechanics.
We define the full photon propagator as where V is the space-time volume, L 3 × T .

Fermion two-point function
The two-point function of the fermions is computed as whereD(p, q) is the Dirac operator in the momentum space. The propagator should have poles at the on-shell energy: where Z 2 denotes the renormalization factor for the fermion wave function. Doubler poles are included as we do not try to eliminate them in the action. We define the full fermion propagator as

Three-point function
We are interested in the three-point function in the momentum space, The photon momentum is denoted as k, and two external fermions have the momentum p and p + k, respectively. The three-point function includes diagrams where the external photon is not attached to the line of the external fermions (in the language of Feynman diagrams).
Among such diagrams, one-particle irreducible ones (the light-by-light scattering) start to appear from the e 7 order, i.e., at the three-loop level. Since we are interested in the threeloops and higher coefficients as well, we use the three-point function defined above rather than D −1 γ µD −1 , which is commonly used in the lattice calculations of flavor non-singlet form factors.
The effective vertex function can be defined by amputating the external legs as The prefactor κ is a product of the appropriate wave-function renormalization factors, which are not needed to be specified in the following discussion. For on-shell fermions, s(p) 2 +m 2 f = 0 and s(p + k) 2 + m 2 f = 0, the vertex function can be expressed by form factors as where u(p) and u(p + k) are fermion wave functions, and σ µν = (i/2)[γ µ , γ ν ]. We used the equation of motion to obtain this formula. On the Euclidean lattice, one cannot directly go to the on-shell momentum because of s(p) 2 ≥ 0. We will discuss how on-shell values can be extracted in the next section.
In the above expression, the O(a 2 ) correction contains the violation of the Lorentz invariance on the lattice. At the tree level, the form factor F 1 has the form F tree in the unit of a = 1 and thus depends on the direction. In the continuum limit, it is reduced to F tree 1 (0) = 1 and thus the Lorentz invariance is recovered.
The physical QED coupling, e P , is defined such that in the continuum limit to all orders in perturbation theory. This should coincide with the definition in terms of the long-range Coulomb force in Eq. (32).
The g-factor to represent the magnetic moment is given by Our goal is to express this quantity as a power series of the renormalized coupling e P . At the tree level, F 2 (0) = 0, and thus it seems automatic to obtain g = 2. However, as we will see later, we evaluate the numerator and the denominator separately by projecting Γ µ in the spatial and the time directions, respectively. Since F 1 (0) depends on the direction as already noted, we do not obtain g = 2 exactly at finite lattice spacings in our calculation.
The Ward identity:

Magnetic moments
We discuss how the ratio of the form factors at the on-shell fermion momenta are extracted from the correlation functions calculated on the Euclidean lattice. In the following discussion, we take a particular configuration of the external momenta p = −k/2 and k 4 = 0. The fermion energy p 4 is left as a variable.
Since we take the ratio in Eq. (42), we do not need to evaluate the normalization factor, κ, in Eq. (39). We first remove the photon propagator from the three-point function G µ in Eq. (38) byĜ We also defineĜ by using the fermion two-point functions. Since we choose p 4 = p 4 + k 4 and p = −(p + k), the two-point functions have a relation, by parity transformation. Therefore, in the numerical calculation, we only need to evaluatẽ We scan all possible values of p 4 on the lattice, and calculate the Fourier transform of appropriately projected functions: and F M (t) = dp 4 2π We note that the above integrals are understood as the abbreviation of the loop sums [cf. Eq. (25)]. We repeat the same calculation forĜ By defining a function g(t) as the g factor is obtained from g(t) in the limit of t → ∞ andk → 0.
One can understand this formula as follows. First, there are useful formulae, where the projection operators, P E and P M , are defined as and Here E and m f are the fermion on-shell energy and the pole mass, respectively, that satisfy Under the projection operators, the equation of motion can be used, and thus we obtain where is the location of the physical pole, and the pole at −z * represents that of the doubler degree of freedom. The function in the numerator, w(z), is a regular function at z = ±z * . Each of the fermion propagators has the simple poles at the same locations, z = ±z * , because of the relation, p 2 = (p + k) 2 . Since the matrices and in the projection operators are residues of the external fermion propagators, we see that the integrals in Eqs. (47) and (48) where · · · are terms which are subleading at large t. There are contributions from the cut, The Fourier transform of the three-point function is related to the position-space quantity as dp 4 2π where The photon fieldÃ µ is kept in the momentum space. The right-hand-side of Eq. (62) is independent of t 2 by the translational invariance. The phase factor e ik 4 t 2 can be ignored for  (t) are, on the other hand, vanishing for odd and even t, respectively. This is problematic when we take the ratio for each t in Eq. (49). In order to avoid the problem, we replace the vanishing Fourier coefficients with the geometric means of the next and previous coefficients before taking the ratio. The factor of two caused by the doubler properly cancels in the quantity (49).

Parameter choices
There are several requirements for the parameters in order to approximate the continuum theory. It is summarized as With the first condition 1/T m γ a, we can get rid of the fermion-photon intermediate states Each of them represents the size of errors from discretization, the finite photon momentum, and the finite photon mass, respectively. We also need to take a large enough T such that

Numerical simulations
We perform a simulation by keeping up to the e 7 order, which corresponds to the three-loop level in the Feynman diagrams. A set of the Langevin equations are solved by the Runge-Kutta improved algorithm [41] explained in Appendix B. We also employ the momentum dependent step size to avoid the critical slowing down in the evolution of the low-momentum modes [42]. The Langevin time τ is discretized with an interval times the momentum dependent factor, 3 e 2 + · · · . In that calculation, we use the Z 3 factor evaluated atk 2 =k 2 , wherek is the external photon momentum of the three-point function. The treatment is justified since we sendk → 0 in the L → ∞ limit.
We show in Fig. 1 the function g(t) in Eq. (49) as a function of 1/t for the 24 3 × 48 lattice. The points are shown for 3 ≤ t < L/2. At each order, coefficients, g (2n) (t)/2, of the expansion, are plotted. We see the expected behavior discussed in the previous section, where deviation can be seen at 1/t 1/12 ∼ 0.08. We extrapolate to the t → ∞ limit, i.e., a (2n) , as follows. We identify the linearly scaling region; the data at too small t are affected by subleading t dependence, whereas those at too large t are contaminated by contributions from backward propagations. We then select the points to be used in the extrapolation to 1/t → 0 as t = T /4 − 3, T /4 − 2, T /4 − 1, and T /4. We perform two linear extrapolations; one line connects the two even points among these four points, and the other does the odd points.
We then take an average of the two extrapolated values. As we discussed before, even and odd points are calculated in a different manner due to vanishing Fourier coefficients by the doubler contributions. We see in the figure that the even and odd points separately behave as linear functions but eventually merge at large t. We take the difference between the even and odd extrapolations as a systematic error. By this procedure, one obtains the g factor at We see that our parameter choice at L = 24 is not too far from the continuum limit.
Repeating the same calculations for other sizes of lattices, one can see the L dependence of the g factors. The results are listed in Table 2 Table 1. The first and second parentheses represent statistical and systematic errors, respectively. We estimate the total errors as the sum of two numbers. The values in the continuum theory are listed in the last row.
actual discretization error is O(m 2 f a 2 ) rather than O(m 2 a 2 ) ∼ O(π/L). The size of the logarithmic correction is controlled by the UV cut-off parameter Λ UV in the action. We take (Λ UV a) 2 = 2.0 which makes the logarithmic correction to be of O(1). For Λ UV a → ∞, the continuum limit gets far away as we discuss later.
We show in Fig. 2 the L dependence of the g factor coefficients, a (2n) , for each order of perturbations. The horizontal lines are again the known analytic results in the continuum theory. The grey lines represent the typical sizes of the discretization errors, ±(ma) 2 . We see a good agreement with the analytic results within the expected uncertainties. At the tree level, there is no statistical error even though we use A 3 (k 2 ) = e −2k 2 /Λ 2 UV , for the tree-level value. The deviation from the continuum theory seen in the figure is due to the O(a 2 ) correction in the fermion-fermion-photon vertex in the lattice Dirac operator. We do not try to obtain the extrapolated values in the 1/L → 0 limit as a reliable estimation seems to require points with larger lattice sizes to control the logarithmic correction of the form (1/L) ln L.
To confirm the validity of our calculation, e.g., the validity of the discretized Langevin steps, we also perform one-loop calculations with the method using Feynman diagrams. First, we compare the wave-function renormalization factor of the photon, Z 3 . We obtain Z −1 3 = 1.30728 + 0.563(3)(e 2 /(4π 2 )) + · · · with the stochastic method, whereas the diagrammatic calculation gives Z −1 3 = 1.30728 + 0.555891(e 2 /(4π 2 )) + · · · . While discrepancy beyond the statistical error is seen in the one-loop coefficient, agreement at the level of a percent is confirmed. The vacuum polarization is obtained by the expansion of the fermion determinant in the stochastic calculation. One may need to increase the number of Langevin steps and/or to generate more noise spinors in Eq. (14) if better precision is necessary. The function g(t) in the 12 3 × 24 lattice with a diagrammatic method is also compared with the stochastic result in Fig. 3 at the one-loop level. We find good consistency within the statistical uncertainties.
The a (2n) factors depend on the UV cutoff parameter Λ UV at finite lattice spacings. We show in Fig. 4 the Λ UV dependence of the one-loop coefficient, a (2) , on the 16 3 × 32 lattice.
We take (Λ UV a) 2 to be 2.0, 8.0, and 32.0, respectively, while keeping other parameters to be the same as the ones listed in Table 1. The same numbers of configurations (6,400 confs) are analyzed for each (Λ UV a) 2 . The grey lines represent the typical discretization errors, ±(ma) 2 , as before. In principle, one can take any value of Λ UV a = 0 as long as we take the continuum limit in the end, but too large Λ UV a makes the continuum limit far away as we can see from the figure. This demonstrates that the introduction of a UV regulator is essential in this calculation, otherwise we would need much larger lattice volumes to see the convergence to the continuum values.
One may also employ the gradient flow [43] to reduce the UV fluctuations rather than modifying the action. Also, one can suppress the UV modes by modifying the noise term in the Langevin equation to be momentum dependent as discussed in Ref. [44]. Those two treatments are equivalent to the method of the modified action at the tree level, and may be useful in QCD or in general theories where a gauge invariant UV regularization cannot be simply implemented. An advantage of the modified action is that one can manifestly maintain the equivalence with the path integral formalism, and thus the way to going back to the original theory is theoretically clear.

IR and finite volume effects
As we modify the QED action to include the finite photon mass, the gauge invariance is broken. In particular, the effects of the photon zero mode should be carefully treated. We also have finite volume effects in the calculation of the function g(t). We discuss those issues in this section.
We first discuss the effects of the photon zero mode. Although we assumed thatĜ µ and G (norm) µ in Eqs. (44) and (45)  However, these seemingly problematic contributions are cancelled in the function g(t) of Eq. (49). This can be understood from the fact that the contributions from the photon zero mode is factorized as an overall factor ∝ e −xt 2 with x = e 2 /(2L 3 T m 2 γ ) [39,40] in the fermion two-point functions and the three-point function we consider.
The effects of the backward propagation in the correlation functions are more serious in the extraction of the g factor we developed. For instance in Eq. (47), related to the symmetry of the loop sum under t → T − t, the leading t dependent part should be accompanied with the backward propagating contribution: By ignoring the backward propagation, one obtains F E (t) ∼ te −Et up to the e −xt 2 factor which is to be canceled in the ratio. The extra terms disturb this behavior. In principle, one can fit the function g(t) by using the above expected behavior. Rather than trying such an analysis, for simplicity, we choose a strategy to look for a region of the extrapolation where the extra term is not important. In future works, it is worth trying such an analysis.
The extra contributions are more and more important for higher orders. In the last line, we express the correction term as a perturbative series. The fermion energy as well as the factor of the photon zero mode are expanded as E = E 0 + e 2 E 1 + · · · and e −xT (T −2t) = 1 − (e 2 /(2L 3 m 2 γ ))(T − 2t) + · · · , respectively. One can see that the powers of (T − 2t) grow as (T − 2t) n at the e 2n order. Since there is an overall suppression factor e −E 0 (T −2t) , one can avoid this problem for sufficiently large T while keeping the region of extrapolation to be t T /2. Nevertheless, for a fixed T , there is a certain perturbative order beyond which the extrapolation by Eq. (70) is not reliable.
Such a limitation is already observed in Fig. 1, where the deviation from the straight line starts to show up at smaller t for higher orders. In order to visualize the deviation, we show in Fig. 5  with the estimated error given in Table 2. At the three-loop order, the plateau seems to start disappearing. A simulation with larger T will be necessary to perform a more reliable  Table 2.
estimate, especially for computation beyond three loops.

Summary
We apply the numerical stochastic perturbation theory to QED and evaluate the perturbative coefficients of the g factor of the electron up to O(α 3 ). We demonstrate that the method reproduces the known results from the diagrammatic method, and expect that it may be able to confirm or predict the higher order perturbative coefficients. Since the method does not rely on the Feynman diagrams, going to higher orders is a matter of computational time and the memory footprint. For example, the increase of the computational cost to accumulate the same number of configurations at the five-loop level is about a factor of two compared to the present work. Although we need more statistics as well as more data points at larger lattice volumes to approach the continuum limit with controlled systematic effects, the results obtained on a small-scale computer in this work are quite encouraging.
The numerical stochastic perturbation theory has been successfully applied to QCD to compute the quantities such as the vacuum energy density and the heavy-quark self-energy, which can be naturally defined on the Euclidean space-time. The challenge we faced in this work is to obtain the physical quantity, g − 2, defined for an electron on its mass-shell.
Thus, an analytic continuation from the Euclidean momentum to the on-shell electron is inevitable. We introduce the Cauchy integral of the vertex function in terms of the imaginary temporal momentum p 4 and design an appropriate combination, from which the g-factor can be obtained at large Euclidean time separations.
The main practical limitation of the method is to separate various energy scales including the electron mass, the external photon momentum and the photon mass. The latter two of these have to be sent to zero in order to obtain the physical result. The photon mass is introduced to avoid an infrared divergence due to finite volume, and it also plays the role to lift the energy of multi-particle states made of an electron and photon so that the single electron state can be easily isolated. Therefore, the extrapolation to the massless photon limit has to be done with great care. In this work, we set the parameters such that these energy scales are separated equally in the logarithmic scale. The systematic error then scales rather slowly, ∼ 1/L, as a function of the lattice extent L. The largest lattice volume in this work is L = 24; we expect that the simulation at L ∼ 128 would be possible on the large-scale machines currently used for modern QCD simulations.
For the complete evaluation of the g − 2 of the electron or the muon, we need to in-clude fermions with different masses to incorporate the mass-dependent contributions. It is straightforward to include those in the Lengiven equation, but having the mass hierarchy as well as the separation of different energy scales poses an extra complication to the parameter setting and extrapolation as described above. We leave such investigations for future works.
Of course, the stochastic method would not replace the diagrammatic computations due to limited accuracy. Especially, for higher orders, in addition to growing statistical uncertainties, the logarithmically enhanced discretization may become an important limiting factor, because they cannot be easily eliminated by numerical extrapolations. Nevertheless, even with, for example, of order of ten percent accuracy, it would provide useful estimate of the higher order coefficients for which the diagrammatic calculation is difficult to reach.
This work provides an evidence that the stochastic quantization method is practically useful in the perturbative calculation of physical quantities, and it may be more effective than the Feynman diagrams at higher loops. We expect that the range of the applications to be much wider than QCD and QED. In general, once we encounter a problem which is sufficiently complicated in the diagrammatic approach, the stochastic method may provide an efficient means to approach. Here λ is an anticommuting parameter. The fields B, c, andc are the Nakanishi-Lautrup field, the ghost field, and the anti-ghost field, respectively, with the action: where the symbol "ˆ" represents the fields multiplied by the UV regulator,Φ(n) = e −∇ 2 /Λ 2 UV Φ(n). The gauge fixing term is obtained by integrating out the B field. The equation of motion holds as the operator equation, i.e., Here D µν , S(p), and G µ (p, k) are defined in Eqs. (34), (37), and (38), respectively.
The first identity (A.13) implies that the longitudinal part of the propagator is exact at tree level, and thusk µ D −1 µν (k) =k Together with the second identity, one finds that Eq. (A.16) is maintained even with the photon mass term.

Appendix B Numerical integration of the Langevin equation
We are going to solve the following type of equation, where g(n, τ i ) are independent Gaussian noises with the variance σ 2 = 1. By taking the limit of small ∆τ , one can integrate the Langevin equation.
A simple improvement is possible by using the Runge-Kutta type algorithm [41].