The electron self-energy in QED at two loops revisited

We reconsider the two-loop electron self-energy in quantum electrodynamics. We present a modern calculation, where all relevant two-loop integrals are expressed in terms of iterated integrals of modular forms. As boundary points of the iterated integrals we consider the four cases p 2 = 0, p 2 = m 2 , p 2 = 9 m 2 and p 2 = ∞ . The iterated integrals have q -expansions, which can be used for the numerical evaluation. We show that a truncation of the q -series to order O ( q 30 ) gives numerically for the ﬁnite part of the self-energy a relative precision better than 10 − 20 for all real values p 2 / m 2 .


Introduction
The two-loop contribution to the electron self-energy in quantum electrodynamics (QED) has been computed for the first time by A. Sabry [1] in 1962.Already at that time it was noticed that the calculation involves certain elliptic integrals.For lack of better techniques at that time the integrands have been approximated by series expansions, and the analysis has been restricted to the region above the threshold.We are now in a better position: Feynman integrals associated to elliptic curves are now a current topic of research interest and our techniques to compute these integrals have evolved .It is therefore natural to revisit the two-loop electron self-energy in QED and to present the result in a modern language.The two-loop electron self-energy is of course the central piece for the determination of the α 2 -term of the electron mass renormalisation constant Z m and the electron field renormalisation constant Z 2 .Let us note that in the MS-scheme only the pole terms are relevant whereas in the on-shell scheme all loop integrals are evaluated on-shell at p 2 = m 2 .In both cases, we do not encounter elliptic integrals for the determination of the renormalisation constants.In this paper we are interested in the finite part of the twoloop electron self-energy for generic values p 2 /m 2 where elliptic integrals do occur.We view the electron self-energy as a (off-shell and gauge-dependent) building block entering two-loop scattering amplitudes in QED.We are interested in an analytic expression for the two-loop selfenergy and in efficient methods for the numerical evaluation of the occurring transcendental functions.As a new result our final answer allows a numerical evaluation with arbitrary precision for all real values of p 2 /m 2 .Other methods for the numerical evaluation of some of the relevant Feynman integrals have been discussed in [36][37][38][39][40].
The electron self-energy is a gauge-dependent object.We perform the calculation in a covariant gauge with gauge parameter ξ.
The renormalised electron self-energy is a renormalisation-scheme dependent quantity.For this reason we present the bare electron self-energy independently of the contributions from the counterterms.In QED the on-shell scheme is the conventional choice.We give the contributions from the counterterms in the on-shell scheme.Renormalisation removes ultraviolet divergences, infrared divergences remain or are introduced through infrared poles in the renormalisation constants.The latter case already occurs at one-loop in the on-shell scheme in QED.
The bare electron self-energy is not a pure function.We may associate a weight to the iterated integrals and the transcendental constants appearing in our calculation.A function is pure if each term in the ε-expansion has uniform weight, where ε denotes the dimensional regularisation parameter.Although we may choose our master integrals as pure functions, the bare electron self-energy is not pure (and is not expected to be pure).
We follow the standard approach for loop calculations: The two-loop electron self-energy is expressed as a linear combination of master integrals.The coefficients of the master integrals are given in appendix E (for Feynman gauge ξ = 1) and in a supplementary electronic file attached to this article (for a general covariant gauge ξ = 1).The master integrals are computed from their differential equations [17,[41][42][43][44][45][46][47][48][49][50][51].For the case at hand, they can be expressed as iterated integrals of modular forms [18].In order to give the reader some orientation on the level of complexity of the calculation, let us classify loop integrals along two criteria: (i) the number of the involved dimensionless variables and (ii) whether the loop integrals may be expressed in terms of multiple polylogarithms or not.The simplest case is given by loop integrals depending on a single dimensionless variable and expressible in terms of (multiple) polylogarithms.Loop integrals associated to two-loop 2 → 2-scattering amplitudes in massless theories are a typical example [52][53][54][55][56][57][58][59], and the class of functions reduces to the sub-class of harmonic polylogarithms.The next difficult case is given by loop integrals depending on several dimensionless variables and expressible in terms of multiple polylogarithms.Loop integrals associated to the two-loop scattering amplitudes for the process e + e − → qg q in massless QCD are an example [60][61][62].In both cases, the standard approach is to transform the system of differential equations to an ε-form [47] through an algebraic change of the kinematic variables and an algebraic change of the basis of the master integrals.If we leave the class of multiple polylogarithms, we first have the case of Feynman integrals beyond the class of multiple polylogarithms, but depending on a single dimensionless variable.This is the case discussed in this paper.Of course, there is also the case of Feynman integrals beyond the class of multiple polylogarithms and depending on several dimensionless variables.An example for the latter would be given by the two-loop integrals associated to the process gg → t t [17,28,29,63,64].Let us emphasize that for the case of interest of this article we are nevertheless able to transform the system of differential equations to an ε-form [27], albeit through a non-algebraic change of the basis of the master integrals.
For the numerical evaluation of the master integrals we switch from the variable x = p 2 /m 2 to a variable q, related to the nome squared of an elliptic curve.We may choose q such that q vanishes at one of the cusps x ∈ {0, 1, 9, ∞}.Let us call this point j and the set of the remaining points S j = {0, 1, 9, ∞}\{ j}.The master integrals have a series expansion in q, which converges for all values x ∈ R\S j , i.e. everywhere on the real line except for three points.The q-expansion is a highly efficient method to evaluate numerically the master integrals for q close to zero.We give q-expansion for the master integrals for all possible choices j ∈ {0, 1, 9, ∞}, thus covering the full kinematic range x ∈ R with efficient numerical evaluation routines and thereby generalising the results of [19].
Efficient numerical evaluations are often based on fast convergent series expansions.The difference between our result and the original work of Sabry lies in the variables used to expand in.With our choice the series show a significant faster convergence (and are well motivated from the underlying mathematics).Quite recently, a purely numerical approach of solving the system of differential equations has been advocated in [65,66], based on expansions around the singular points of the differential equations.The differences and the similarities with our approach are as follows: Within our approach we first obtain an analytic result in terms of iterated integrals of modular forms (with a notion of weight for these iterated integrals, such that the ε j -term of each master integrals has uniform weight j), and only in a second step we use efficient numerical methods for the numerical evaluation of these iterated integrals.The approach of [65,66] is numerical from the beginning.We use a variable q as expansion parameter, whereas the method of [65,66] uses p 2 /m 2 or a rational function of this variable.Common to both methods is the expansion around all singular points of the differential equations, and -on a technical levelthe determination of the boundary constants for the expansion around the second, third and any further singular point from the boundary constants for the expansion around the first (or any other already known) singular point with the help of high-precision numerics and the PSLQalgorithm [67].This paper is organised as follows: Section 2 briefly reviews the Lagrange density of quantum electrodynamics and the (known) renormalisation constants to two-loop order.In section 3 we introduce the master integrals for the two-loop electron self-energy.Section 4 is a brief introduction to iterated integrals.In section 5 we express the bare two-loop electron self-energy in terms of master integrals.In addition we give the counterterms from renormalisation.In section 6 we evaluate the master integrals in terms of iterated integrals of modular forms.The iterated integrals start at a boundary point j and we do this for the four choices j ∈ {0, 1, 9, ∞}.The iterated integrals have a q-series expansion, which may be used for the numerical evaluation.In section 7 we present numerical results.Our conclusions are given in section 8.The appendix contains useful information on the QED Feynman rules (appendix A), the numerical computation of the complete elliptic integral with the help of the arithmetic-geometric mean (appendix B), the Eisenstein series of modular weight 1 appearing in the calculation (appendix C), the one-loop electron self-energy (appendix D), the coefficients appearing in the expression of the two-loop electron self-energy in terms of the master integrals (appendix E), the boundary constants for the four cases j ∈ {0, 1, 9, ∞} (appendix F) and a description of the content of the supplementary electronic file attached to this article (appendix G).

Lagrange density and renormalisation
The bare gauge-fixed Lagrangian for QED reads in a covariant gauge We use dimensional regularisation and set D = 4 − 2ε.Under renormalisation one redefines the fields and the parameters The arbitrary scale µ is introduced to keep the renormalised coupling e dimensionless.The factor S ε = (4π) ε exp(−εγ E ) absorbs artefacts of dimensional regularisation (logarithms of 4π and Euler's constant γ E ).For convenience we set Substituting the above relations into the Lagrange density we obtain where L renorm is given by L QED where all bare quantities are replaced by renormalised ones (the bare coupling e 0 is replaced by e (D) ).The counterterms are given by The renormalisation constants are known to very high loop order [68][69][70][71][72][73][74][75][76][77][78][79].Here, we only need them to order O(α 2 ).We write In the on-shell scheme we have [71] Z For the two-loop contribution to the electron self-energy we need Z 3 to order α and Z 2 and Z m to order α 2 .

Definitions and notation for the master integrals
There are three Feynman diagrams contributing to the two-loop electron self-energy in QED.These diagrams are shown in fig. 1.We label these diagrams "rainbow diagram", "kite diagram" and "fermion insertion diagram", respectively.The second diagram (the kite diagram) can be drawn equivalently as shown in fig.2, motivating the name "kite diagram".The master integrals for the rainbow diagram and the fermion insertion diagram will be a subset of the master integrals for the kite diagram.To see this, let us first introduce the kite integral.We set  with the propagators 10) and ν 12345 = ν 1 + ν 2 + ν 3 + ν 4 + ν 5 .The internal momenta are denoted by k 1 and k 2 , the internal mass by m, the external momentum by p and the dimension of space-time by D = 4 − 2ε.The arbitrary scale µ renders the integral dimensionless.In the following we set µ = m.We further define The five propagators D 1 -D 5 are indicated by the numbers 1-5 in fig. 2. We note that all propagators of the rainbow diagram and the fermion loop insertion diagram are a subset of these.In order to show this, we labelled all propagators in fig. 1 with the appropriate numbers.Therefore it is sufficient to consider only the master integrals of the kite integral.In order to present these master integrals let us first denote by ψ 1 and ψ 2 two independent solutions of the second-order differential equation [2,36] x Of course, this does not fully specify ψ 1 nor ψ 2 , but for the moment this is all what we would like to assume about ψ 1 and ψ 2 .The exact definitions of ψ 1 and ψ 2 will be given in section 6.We denote the Wronskian by We will normalise ψ 1 and ψ 2 such that We will see later that eq. ( 12) is the Picard-Fuchs equation for the periods of an elliptic curve and ψ 1 and ψ 2 will be taken as periods of an elliptic curve.We denote the ratio of the two periods and the nome squared by We further set Let us now return to the master integrals.There are eight master integrals, which we take as [13,27] In the master integral I 6 the sunrise integral in D = 2 − 2ε space-time dimensions appears.Using dimensional-shift relations we may express this integral in terms of integrals in D = 4 −2ε spacetime dimensions.We have The master integral I 7 involves the derivative of I 6 .We have and where the 8 × 8-matrix A is independent of ε and given by The entries of the matrix A are as follows: We first define and set then All entries may be expressed as polynomials in We have Let us stress that all formulae in this section are valid for any choice of ψ 1 and ψ 2 , as long as these are two independent solutions of eq. ( 12) and normalised such that eq. ( 14) holds.We will use this freedom to define four sets of master integrals, which we denote by I i, j (with 1 ≤ i ≤ 8 and j ∈ {0, 1, 9, ∞}), corresponding to the different q-expansions around the four cusps x ∈ {0, 1, 9, ∞}.Where it is not relevant, we will drop the additional index j.The first five master integrals and the last one are identical in all four sets and differ only in the variables they depend on.However, in the definition of the sixth and the seventh master integral the period ψ 1 appears explicitly and these integrals are not identical in the four sets.Of course, the dependence on our choice of ψ 1 is also reflected in the coefficients expressing the two-loop self-energy as a linear combination of the master integrals, such that the final result is independent of the choice of ψ 1 and ψ 2 .To cut the story short: This setup allows us to use in any region a choice of variables with the best numerical convergence.

Iterated integrals
We may easily solve the differential equation in eq. ( 21) order-by-order in ε.The solution is expressed in terms of iterated integrals.Let us briefly review iterated integrals [80].For differential 1-forms ω 1 , ..., ω k on a manifold M and a path γ : The iterated integral is defined by Harmonic polylogarithms are a special case of iterated integrals [81].We consider two integration kernels and define the harmonic polylogarithms by and The last equation defines harmonic polylogarithms for trailing zeros.
A second special case are iterated integrals of modular forms.Let f 1 (τ), f 2 (τ), ..., f k (τ) be modular forms of a congruence subgroup.and assume that f k (τ) vanishes at the cusp τ = i∞.We define the k-fold iterated integral by The case where f k (τ) does not vanishes at the cusp τ = i∞ is discussed in [18,82] and is similar to trailing zeros in the case of harmonic polylogarithms.If we change the integration variables from τ to q we obtain with Given the q-expansion of the modular forms f 1 , ..., f k , we may easily obtain the q-expansion of the iterated integral F ( f 1 , f 2 , ..., f k ; q) by integrating term-by-term and multiplication of power series.
5 The two-loop self-energy

The bare two-loop self-energy
We first compute the bare two-loop electron self-energy in QED.We write with α = e 2 /(4π) for the bare self-energy separating the part proportional to / p and the part proportional to m.The quantities Σ bare,V and Σ (2) bare,S are expressed as linear combinations of the master integrals I 1 -I 8 : bare,S = We work in a general covariant gauge with gauge parameter ξ.The coefficients c V j and c S j are rational functions in x, ε, ξ, ψ 1 /π and 1/π • dψ 1 /dx.For ξ = 1 (Feynman gauge) they are listed in appendix E. For a general covariant gauge (ξ = 1) they are given in a supplementary electronic file attached to this article.The master integrals I 1 -I 8 satisfy a differential equation in ε-form, therefore the ε-expansion of Σ bare,S is easily obtained from eq. ( 36) by expanding in ε to the desired order.The pole terms are rather simple.Let us write We have The finite parts Σ (2,0) bare,V and Σ (2,0) bare,S will be given in section 6.

The counterterms from renormalisation
We write for the counterterms from renormalisation at order O(α 2 )
The relevant diagrams are shown in fig. 3. Σ CT,V and Σ CT,S are given by CT,S = Z (2) Figure 3: The Feynman graphs corresponding to the counterterms.In the first graph we take the O(α 2 )-term.
J 1 and J 2 are the two master integrals of the one-loop calculation.They are given in appendix D.
The renormalisation constants in the on-shell scheme are given in eq. ( 8).In the on-shell scheme we have CT,S = (2 − 6ξ) 6 Evaluation of the master integrals

Transcendental constants
There are a few transcendental constants appearing in the boundary constants.The transcendental constants relevant to our calculation are weight 1 : π, ln (2) , ln (3) , and products of those.Our convention for the notation of the multiple polylogarithms is Li ...
The Clausen function is defined by thus where r 3 = exp(2πi/3) denotes the third root of unity.We list the relevant boundary constants for the boundary points x ∈ {0, 1, 9, ∞} in appendix F.

Integrals expressible in terms of harmonic polylogarithms
The master integrals I 1 -I 5 may be expressed (to all orders in ε) in terms of harmonic polylogarithms [13,14].The first few orders of the ε-expansion of the integrals I 1 -I 5 may be expressed in terms of classical polylogarithms.We have The analytic continuation and the numerical evaluation of these master integrals are well understood.The analytic continuation is dictated by Feynman's iδ-prescription: x → x + iδ.There are packages, which allow the numerical evaluation of harmonic polylogarithms H m 1 ...m k (x) in double precision and arbitrary precision [83][84][85][86].
Let us also note the following alternative: The harmonic polylogarithms with the letters f 0 and f 1 can be written as iterated integrals of modular forms.Since the latter are required anyhow for the problem at hand, we may as well treat the analytic continuation and the numerical evaluation within the context of iterated integrals of modular forms.We will follow this approach in this paper.

The elliptic master integrals
The master integrals I 6 -I 8 depend on elliptic topologies and may be expressed as iterated integrals of modular forms.We remark that also the master integrals I 1 -I 5 may be written as iterated integrals of modular forms.This follows from the relations From the maximal cut of the sunrise integral we obtain the elliptic curve We denote the roots of the quartic polynomial in eq. ( 48) by There is an isomorphism between the elliptic curve and C/Λ, where Λ is a lattice generated by the periods of the elliptic curve: where ψ 1,0 and ψ 2,0 are two periods of the elliptic curve generating the lattice Λ.It can be shown that ψ 1,0 and ψ 2,0 are two independent solutions of the homogeneous second-order differential equation given in eq.(12).Any other pair ψ 1, j and ψ 2, j of periods related to the first one by generates the same lattice.ψ 1, j and ψ 2, j are again two independent solutions of the differential equation given in eq. ( 12).If ψ 1,0 and ψ 2,0 are normalised according to eq. ( 14), then so are ψ 1, j and ψ 2, j .We see that we have some freedom in choosing a pair of periods ψ 1, j and ψ 2, j as independent solutions of eq. ( 12).We will label different choices by the subscript j.The definition of the master integrals I 6 and I 7 involves the period ψ 1, j and hence depends on our choice of ψ 1, j and ψ 2, j .Let us now discuss the dependence on our choice of ψ 1, j and ψ 2, j in more detail: As already mentioned, the master integrals I 1 -I 5 and I 8 do not depend at all on our choice of ψ 1, j and ψ 2, j : The dependence of I 6 on our choice is rather simple: The dependence of I 7 on our choice is more tricky, due to the appearance of the derivative of dψ 1 /dx.One finds This relation is most easily derived by relating both the basis I i, j ′ and the basis I i, j to a basis independent of our choice of periods.This intermediate basis does not need to be an ε-basis.A possible intermediate basis is Eqs. ( 52)-( 54) allow us to match the master integrals with different choices of periods in regions where both choices lead to convergent series expansions.This can be used to determine the boundary conditions for one choice from the known boundary conditions of the other choice.
In practice we proceed as follows: Suppose we already know the boundary constants for the choice j and we would like to obtain the boundary constants for the choice j ′ .Suppose further that there is a region where both choices lead to convergent series expansions.Evaluating both expressions to high precision gives us numerical values to high precision for the boundary constants of the choice j ′ .We may then use the PSLQ-algorithm [67] to match these values to a Q-linear combination of the transcendental constants from section 6.1.
In the following we will discuss four choices for the pair of periods.We label them For each choice we set The values of n j will be Each of the four choices has the property that q n j , j = 0 for x = j, j ∈ {0, 1, 9, ∞}, i.e. q n j , j vanishes at the cusp x = j.We write b i, j if ψ 1, j is substituted for ψ 1 in the generic definition of b i in eq. ( 25).For the ε-expansion of the master integrals we write We will need for the two-loop self-energy up to the finite part the integral I 6, j to order ε 2 , the integral I 7, j to order ε and the integral I 8, j to order ε 3 .The integral I 6, j starts at order ε 2 , the integral I 7, j starts at order ε, while the integral I 8, j starts at order ε 3 .Therefore we need in all three cases the first non-vanishing order.

6.3.1
The cusp p 2 = 0 We start with the cusp x = 0. We introduce the modulus k and the complementary modulus k ′ through Explicitly we have where Feynman's iδ-prescription (x → x + iδ) is understood.Our choice of periods for this case (which agrees with the choice made in ref. [27]) is given by where K(k 0 ) denotes the complete elliptic integral of the first kind.The complete elliptic integral is efficiently computed with the help of the arithmetic-geometric mean, reviewed in appendix B. The 2 × 2-matrix γ is given by x The values of the variables τ 1,0 , τ 6,1 , τ 2,9 and τ 3,∞ at the cusps x ∈ {0, 1, 9, ∞}.
The matrix γ ensures that the periods ψ 1,0 and ψ 2,0 vary smoothly as x varies smoothly in x ∈ R + iδ [19].The complete elliptic integral K(k) can be viewed as a function of k 2 : We set K(k 2 ) = K(k).The function K(k 2 ) has a branch cut at [1, ∞[ in the complex k 2 -plane.The matrix γ compensates for the discontinuity when we cross this branch cut.It is relatively easy to see that k 2 as a function of x crosses this branch cut at the point x = 3 − 2 √ 3 ≈ −0.46, the corresponding value in the k 2 -plane is k 2 = 2.The point x = 1 is a little bit more subtle.Let us parametrise a small path around x = 1 by then and the path in k 2 -space winds around the point k 2 = 1 by an angle 3π as the path in x-space winds around the point x = 1 by the angle π.Note that eq. ( 63) defines the periods ψ 1,0 and ψ 2,0 for all values x ∈ R + iδ.The periods take values in C ∪ {∞}.
One easily verifies that ψ 1,0 and ψ 2,0 are normalised according to eq. ( 14).We set The values of τ 1,0 at the points x ∈ {0, 1, 9, ∞} are tabulated in table 1.In fig. 4 we plot the values of the variable q 1,0 as x ranges over R. We see that all values of q 1,0 are inside the unit disc with the exception of the three points x ∈ {1, 9, ∞}, where the corresponding q 1,0 -values are on the boundary of the unit disc.

The cusp
For the expansion around the cusp x = 9 we set We further set The values of τ 2,9 at the points x ∈ {0, 1, 9, ∞} are tabulated in table 1.In order to simplify the notation we write in the remaining part of this section and e 1 = E 1 (τ 2,9 ; χ 0 , χ 1 ), e 2 = E 1 (2τ 2,9 ; χ 0 , χ 1 ).We also use the notation b 1 for b 1,9 and b 2 for b 2,9 .The Hauptmodul is given by The q-expansions of b 1 and b 2 are given by For the integration kernels we have The integral I (3) 8,9 is finite at x = 9.We denote its value at x = 9 by We obtain for the first non-vanishing order of the integrals I 6,9 , I 7,9 and I 8,9

The cusp p 2 = ∞
For the expansion around the cusp x = ∞ we set We further set The values of τ 3,∞ at the points x ∈ {0, 1, 9, ∞} are tabulated in table 1.In order to simplify the notation we write in the remaining part of this section and e 1 = E 1 (τ 3,∞ ; χ 0 , χ 1 ), e 2 = E 1 (2τ 3,∞ ; χ 0 , χ 1 ).We also use the notation b 1 for b 1,∞ and b 2 for b 2,∞ .The Hauptmodul is given by The q-expansions of b 1 and b 2 are given by For the integration kernels we have We obtain for the first non-vanishing order of the integrals I 6,∞ , I 7,∞ and I 8,∞

Numerical results
With the four expansions around the cusps j ∈ {0, 1, 9, ∞} at hand, we now address the question, which expansion to use for a given x.The four expansion parameters are q 1,0 , q 6,1 , q 2,9 , q 3,∞ .
For the absolute values of these we always have where the value 1 is only attained for |q n j , j | at the three points S j = {0, 1, 9, ∞}\{ j}.For a fast convergence we would like to choose j such that |q n j , j | has a small absolute value.In fig.6 we plot the absolute values of q n j , j for the four choices j ∈ {0, 1, 9, ∞}.An appropriate choice is q 3,∞ q 2,9 q 6,1 The absolute values of the variables q n j , j for j ∈ {0, 1, 9, ∞} as a function of x.There is always a choice such that |q n j , j | 0.163.This value is indicated by the dashed black line.
We denote the chosen variable simply by q.In this way we can ensure that for all The value |q| ≈ 0.163 is indicated by a dashed line in fig.6.Let us now discuss the precision, which can be reached by truncating the q-series to a certain order O(q N ).For a quantity O and an approximation O approx to this quantity we define the relative precision δ of the approximation by We consider the ε 0 -terms of the bare two-loop self-energy in Feynman gauge, i.e. the terms Σ (2,0) bare,V and Σ (2,0) bare,S .We may express these two quantities as linear combinations of iterated integrals of modular forms with coefficients, which are rational functions in x, ψ 1 /π and 1/π •  bare,S (right) obtained by truncating the q-series of the iterated integrals at order O(q 30 ) for the various choices q ∈ {q 1,0 , q 6,1 , q 2,9 , q 3,∞ }.
There is always a choice such that the relative precision is below 1.5 × 10 −21 .This value is indicated by the dashed black line.Within the first possibility (method A) we compute the coefficients from the known values of x ψ 1 /π and 1/π • dψ 1 /dx, whereas we approximate the iterated integrals by their q-series, truncated to order O(q N ).Within the second possibility (method B) we compute Σ (2,0)   bare,V and Σ (2,0) bare,S from their q-series, truncated to order O(q N ), i.e. we use the analogue of eqs.( 76), ( 85), (95) or (104) to the appropriate order O(q N ).
Let us consider a truncation of the q-series to order O(q 30 ).In fig.7 we show for method A the relative precision for Σ (2,0) bare,V and Σ (2,0) bare,S by truncating the iterated integrals to order O(q 30 ) for the various choices q ∈ {q 1,0 , q 6,1 , q 2,9 , q 3,∞ } as a function of x.We see that for all values x ∈ R the O(q 30 )-approximation of the iterated integrals gives us a relative precision on Σ (2,0)   bare,V and Σ (2,0) bare,S better than 1.5 × 10 −21 .Fig. 8 shows the corresponding plot for method B. We observe, that in some regions of x the relative precision is only below 1.3 × 10 −11 .The bad regions are the ones where we switch from one choice of q to another, i.e. the regions where the optimal expansion parameter q is close to its maximum |q| ≈ 0.163.These regions are regions away from the singular points x ∈ {0, 1, 9, ∞} of the differential equation.
The q-series for Σ (2,0) bare,V and Σ (2,0) bare,S have their virtues close to the singular points x ∈ {0, 1, 9, ∞}.An inspection of the formulae from appendix E shows, that the coefficients of the iterated integrals have poles like These poles don't show up in the final result for Σ  inside the unit disc.Let us now study how fast or good this series converges.In fig. 9 we plot the relative precision at the point x = −10 if we truncate in eq. ( 112) the q 6,1 -expansion at order N as a function of N. We see that this series converges rather slowly.Only after including roughly 150 terms in the q-expansion we reach a precision of 10%.With a truncation after 250 terms we reach a relative precision of 10 −20 .Let us note that a truncation below 150 terms gives unreliable results.If one only considers truncations up to N = 30 one might be led to the false conclusion that the series is divergent.As already mentioned, this is a constructed worst case scenario.The actual relative precisions for the finite part of the two-loop self-energy for method A and method B have been given in fig.7 and fig.8, respectively.
In fig. 10 we plot the final result for Σ numerically for the finite part of the self-energy a relative precision better than 10 −20 for all real values x.We expect the methods discussed here to be useful also in other precision calculations involving elliptic integrals.
The coupling e (D) is defined in eq. ( 4).The Feynman rules for the counterterms are

B The arithmetic-geometric mean
In this appendix we review the numerical evaluation of the complete elliptic integral of the first kind with the help of the arithmetic-geometric mean.Let a 0 and b 0 be two complex numbers.
For n ∈ N 0 one sets The sign of the square root is chosen such that [87] and in case of equality one demands in addition The sequences (a n ) and (b n ) converge to a common limit known as the arithmetic-geometric mean.The complete elliptic integral of the first kind is given by

D The one-loop self-energy
In this appendix we consider the one-loop electron self-energy.The relevant family of (one-loop) Feynman integrals is given by where the propagators D 1 and D 4 have been defined in eq. ( 10).As before we set µ = m.As a basis of master integrals we use The differential equation for J = (J 1 , J 2 ) T reads ,

F Boundary constants
In this appendix we list the relevant boundary constants when we integrate the differential equation eq. ( 21) from one of the points x ∈ {0, 1, 9, ∞}.For the finite part of the two-loop electron self-energy we need the master integrals I 1 -I 6 to order ε 2 , the master integral I 7 to order ε 1 and the master integral I 8 to order ε 3 .We denote by C k i, j the boundary constant for the ε k -th term of the master integral I i at boundary point j.
c_V and c_S are vectors, giving the coefficients c V j and c S j appearing in eq. ( 36) in a general covariant gauge with gauge parameter ξ.I_case_0 is a vector, giving the results for the eight master integrals I 1 -I 8 with boundary point j = 0 as an expansion in ε to the relevant order (I 1 -I 6 to order ε 2 , I 7 to order ε and I 8 to order ε 3 ) and an expansion in q to order O(q 50 ).I_case_1, I_case_9 and I_case_inf are similar, but for the boundary points j = 1, j = 9 and j = ∞, respectively.Sigma_V_case_0 and Sigma_S_case_0 give the ε 0 -parts Σ (2,0) bare,V and Σ (2,0) bare,S of the bare two-loop self-energy as an expansion in q to order O(q 50 ) in a general covariant gauge.The quantities Sigma_V_case_1, Sigma_S_case_1, Sigma_V_case_9, Sigma_S_case_9, Sigma_V_case_inf, and Sigma_S_case_inf are similar, but for the boundary points j = 1, j = 9 and j = ∞, respectively.

Figure 1 :
Figure 1: The Feynman graphs contributing to the two-loop electron self-energy.

Figure 2 :
Figure 2: The kite graph.This graph is equivalent to the second graph in fig. 1.
which allows us to express all integrals in D = 4 − 2ε dimensions.Let us set I = (I 1 , I 2 , I 3 , I 4 , I 5 , I 6 , I 7 , I 8 ) T .The differential equation for I with respect to τ n reads 1 2πi [a, b] → M let us write for the pull-back of ω j to the interval [a, b]

Figure 6 :
Figure6: The absolute values of the variables q n j , j for j ∈ {0, 1, 9, ∞} as a function of x.There is always a choice such that |q n j , j | 0.163.This value is indicated by the dashed black line.

dψ 1 /
dx.There are now two possibilities to compute numerically the values of Σ

Figure 9 :Figure 10 :
Figure 9: The relative precision of the order q 6,1 N -truncation of the q 6,1 -series of the function f (x) = 1/x 3 at the point x = −10 as a function of N. The dashed line indicates a relative precision of 10%.

Figure 11 :
Figure 11: The real part (left) and the imaginary part (right) of the ε 0 -term of the bare quantity Σ (2) bare,S in Feynman gauge.