Three-loop Euler-Heisenberg Lagrangian in 1+1 QED, part 1: single fermion-loop part

We study the three-loop Euler-Heisenberg Lagrangian in spinor quantum electrodynamics in 1+1 dimensions. In this first part we calculate the one-fermion-loop contribution, applying both standard Feynman diagrams and the worldline formalism which leads to two different representations in terms of fourfold Schwinger-parameter integrals. Unlike the diagram calculation, the worldline approach allows one to combine the planar and the non-planar contributions to the Lagrangian. Our main interest is in the asymptotic behaviour of the weak-field expansion coefficients of this Lagrangian, for which a non-perturbative prediction has been obtained in previous work using worldline instantons and Borel analysis. We develop algorithms for the calculation of the weak-field expansions coefficients that, in principle, allow their calculation to arbitrary order. Here for the non-planar contribution we make essential use of the polynomial invariants of the dihedral group D4 in Schwinger parameter space to keep the expressions manageable. As expected on general grounds, the coefficients are of the form r1+r2*zeta(3) with rational numbers r1, r2. We compute the first two coefficients analytically, and four more by numerical integration.


Introduction
In one of the first serious computations of early QED, Heisenberg and Euler in 1936 obtained their famous effective Lagrangian induced for a constant electromagnetic field by an electron loop [1]. For this effective Lagrangian they obtained the following one-parameter integral representation: (1.1) Here m and T are the mass and proper-time of the electron and a, b the Maxwell invariants, related to the electric and magnetic fields by a 2 − b 2 = B 2 − E 2 , ab = E · B. The second and third term in brackets implement the renormalization of vacuum energy and charge. The superscript '(1)' stands for one-loop. This Euler-Heisenberg Lagrangian ('EHL' in the following) is to be added to the classical Maxwell Lagrangian, and contains a wealth of information on nonlinear quantum effects such as the field-dependence of the speed of light, vacuum birefringence, and dichroism (see [2,3] for reviews). Moreover, it contains this information in a form which is ready for use with standard methods of nonlinear optics. From a particle theory point of view, the EHL holds the information on the one-loop N -photon amplitudes for arbitrary N in the low-energy limit (see, e.g., [3,4]). In terms of Feynman diagrams, it is thus equivalent to the sum of graphs shown in Fig. 1, where all photon energies are small compared to the electron mass, ω i m. Figure 1: Sum of diagrams equivalent to the one-loop EHL.

Advances in High Energy Physics
low-energy limit (since a constant field can emit only zeroenergy photons). Thus diagrammatically L (1) ( ) is equivalent to the sum of the Feynman graphs shown in Figure 1, where all photon energies are small compared to the electron mass, ≪ .
In [10] it was shown how to carry out explicitly the construction of these limiting low-energy amplitudes starting from the weak-field expansion of the EHL: It turned out that if one fixes the number of photons, their momenta 1 , . . . , , and a helicity assignment for each photon, then in this limit the full dependence on the momenta and polarization vectors is carried by a unique invariant. Thus the magnitude of the amplitude can be specified by a single number, which will be essential for our whole approach. Except for the purely magnetic case, the EHL has also an imaginary part related to vacuum pair creation by the electric field component (to be called "Sauter-Schwinger pair creation" in the following) [18,19]. In the purely electric case one finds, from the poles in the -integration, the following decomposition due to Schwinger [19]: Consider the purely magnet powers of the field yields with an effective expansion par coefficients (1) that can be writt numbers : (1) = − 2 2 (2 ) (2 − Here (1) holds information on tudes. The asymptotic behaviou easily studied using well-known numbers. One finds The construction of these low-energy amplitudes from the weak-field expansion of the EHL, (1. 2) can, using standard spinor helicity techniques, be carried out explicitly for any number of photons and helicity assignments [5,6]. It turns out that, in this limit, once a helicity assignment has been fixed for all photons, the full dependence of the N -photon amplitude on the external momenta and polarization vectors can be absorbed into a single invariant χ N . Except for the purely magnetic case, the EHL has also an imaginary part related to vacuum pair creation by the electric field component [7]. In the purely electric case, which we will focus on here for simplicity, this imaginary part allows the following decomposition due to Schwinger [7]: . Physically, the kth term in this decomposition relates to the coherent production of k electron-positron pairs in one Compton volume of spacetime (see, e.g., [8]). The nonperturbative dependence of these "Schwinger exponentials" on the field supports the interpretation of field-induced pair creation as a vacuum tunneling effect, as envisioned by Sauter as early as 1931 [9].
The kth term in the decomposition (1.3) can be obtained from the kth pole in the proper-time integrand of (1.1) by a simple application of the residue theorem. Still, for higher-loop considerations it turns out to be advantageous to observe that the real and imaginary parts of the EHL can be related alternatively through the asymptotic properties of the coefficients of the weak-field expansion (1.2) [10]. In the purely electric case, this expansion is where g ≡ eE (1. 5) with B N the Bernoulli numbers. Replacing the weak-field expansion coefficients c (1) n by their leading asymptotic growth, which is one can Borel sum the series. One obtains a Borel integral that is singular, leading to an imaginary part that is precisely the leading term in the Schwinger decomposition (1.3): This procedure can be repeated with the subleading, sub-subleading etc. asymptotic growth of the coefficients c (1) n to reproduce the full Schwinger expansion in (1.3) [10]. However, -3 -in this paper we will be concerned only with the leading term in this expansion, which dominates in the weak-field limit. Our interest here is in multi-loop corrections to the EHL. The two-loop correction to the EHL, involving one internal photon exchange, corresponds to the diagrams shown in Fig. 2 (here it is understood that internal photon corrections are put in all possible ways).
Advances in High Energy Physics QED -matrix by representing the scalar particles in terms of relativistic particle path integrals, and coupling them through photons in all possible ways. Upon restriction to the purely photonic part of the -matrix (no external scalars) and to the "quenched" contribution (only one virtual scalar), this "worldline representation" can be stated very compactly in terms of the (quenched) effective action Γ[ ]: Here denotes the proper-time of the scalar particle in the loop, its mass, and ∫ ( )= (0) D ( ) a path integral over all closed loops in space-time with fixed periodicity in the proper-time. The worldline action [ ( )] has three parts: They are given by The kinetic term 0 describes the free propagation of the scalar, ext describes its interaction with the external field, and int generates the corrections due to internal photon exchanges in the loop. Expanding out the two interaction exponentials leads back to Feynman diagrams, however with stationary phase approximation. Fo electric field in the direction this e is given by a circle in the (Euclidean It can be shown that in the weak-fi the imaginary (although not the r Lagrangian can be well-approximat integral with this single trajectory: This is easily evaluated to be Their argument assumes the field to restriction on the strength of the c following: (i) Formula (30), if true, cons all-loop summation of an of arbitrary loop order. Tho Table 1.
(ii) According to [26], the co quenched diagrams gets sup limit.
(iii) Perhaps, most surprisingly, t in (30) is already the phys implying that the worldline matically takes all mass reno grams into account. This is that the determination of the This two-loop EHL was first studied by Ritus in 1975 [11]. That calculation, as well as later recalculations [12][13][14], led to a representation of L (2) in terms of rather intractable two-parameter integrals, and so far only the first few coefficients of its weak-field expansion have been computed [10,11,13,15,16].
Schwinger's formula for the imaginary part (1.3) also generalizes to the two-loop level, in the following way [8,11,17]: ). Thus at two-loop one finds the same decomposition into Schwinger exponentials, but the prefactor of the kth Schwinger-exponential now is not a constant but a function K k (β) of the field strength. All that is presently known about these functions K k (β) explicitly are their leading orders in the weak-field epxansion [11,17]. At this leading order, things become extremely simple: adding up the one-loop and two-loop EHL's one finds And in [10] it was checked, albeit only numerically and based on a calculation of fifteen coefficients, that the method of constructing the imaginary part of the EHL by Borel summation of the weak-field expansion coefficients still works fine at the two-loop level, at least in this weak-field limit. However, it works now in a more complicated manner, since at the two-loop level mass renormalization kicks in. It turns out that, before mass renormalization, the coefficients of L (2) have an asymptotic growth that is faster than what we have seen above for the one-loop coefficients; only after adding the counterterm from mass renormalization, which is of the form δm ∂ ∂m L (1) , and only if δm is the physically renormalized mass (renormalized at the one-loop level) the addition of the counterterm leads to a cancellation that reduces the leading asymptotic growth of the two-loop coefficients to be precisely the same as at one loop.
Thus the study of the two-loop EHL allows one to understand already two important aspects of the QED S-matrix: first, although the physically renormalized electron mass is usually determined through a calculation of the electron propagator, it can as well be obtained from the effective Lagrangian (although one has to go one loop higher in the perturbative expansion). This has been called "immanent renormalization principle" by Ritus (see, e.g., [8]). Second, unless the physically renormalized electron mass is used, the two-loop N -photon amplitudes dominate over the one-loop ones for sufficiently large N , i.e. perturbation theory breaks down already at the two-loop level.
In [17] it was further noted that, if one would assume that in this weak-field approximation the higher-loop order corrections just lead to an exponentiation, then the e απ factor can be absorbed into the Schwinger factor e − π β by the following massshift, This extrapolation, to be called "exponentiation conjecture" in the following, would appear far-fetched if based only on a two-loop calculation, but there is strong independent support for its correctness: first, the same mass shift had been found by Ritus already before for the crossed process of electron propagation in a constant electric field [18]. Second, in the vacuum tunneling picture of Sauter-Schwinger pair creation it can be interpreted as taking into account the leading-order Coulomb energy correction due to the fact that the electron -positron pair materializes at a finite distance [17]. Third, the analogue exponentiation has also been obtained for Scalar QED by Affleck et al. [19], and using a totally different approach based on a semi-classical approximation to Feynman's worldline path integral presentation of the multiloop QED effective action [20,21]. In 2002 G.V. Dunne and one of the authors [22][23][24] noted that the electric or magnetic backgrounds are by no means the simplest ones in the Euler-Heisenberg context. Computationally, the most favorable case is the one of a (euclidean) self-dual ('SD') field. Such a field is defined by F =F , which has the consequence that F 2 = −f 2 1l. For real f , the SD effective Lagrangian has properties similar to the magnetic EHL, for imaginary f similar to the electric one. Surprisingly, for such a background the EHL turned out to be computable in closed form even at the two-loop level: with ψ the digamma function ψ(x) = Γ (x)/Γ(x). This self-dual EHL cannot be realized with real fields in Minkowski space, but still holds physical information on the photon amplitudes, namely on their "all +" helicity components [25,26].
Moreover, in the self-dual case the study of the weak-field expansions and the construction of the imaginary parts involve only well-known properties of the digamma function, and thus could be done much more completely and explicitly than for the electric or magnetic case [24]. This made it possible to verify that the above-mentioned construction of the imaginary part by Borel summation works as well at the two-loop level, that is, certain a priori possible complications such as the appearance of additional poles or cuts in the complex g plane apparently do not occur in QED.
Having thus gained confidence that the correspondence between the leading asymptotic behaviour of the weak-field expansion coefficients (1.6) and the leading weak-field behaviour of the imaginary part (1.7) of the EHL persists to higher loop orders, it was suggestive to apply this correspondence in reverse to the exponentiation conjecture to gain information on the multiloop photon amplitudes. This was carried out by G.V. Dunne and one of the authors in [27]. There we essentially transferred the factor e απ first from the leading weakfield behaviour of the imaginary part to the leading large N -behaviour of the expansion coefficients, and from there to the N -photon low-energy amplitudes in the limit of N → ∞. For the second step it was essential that, as mentioned above, in the low-energy limit the whole kinematical dependence of the photon amplitudes can be absorbed into a single invariant, effectively reducing the amplitude to a single number.
However, this leads to an analytic dependence on α, which seems at variance with wellknown arguments [28,29] that exclude a non-vanishing radius of convergence for generic amplitudes in QED. For a further discussion of this apparent paradox see the recent [30].
Considering its relevance for the asymptotic structure of the QED perturbation theory, as well as its relation to the vacuum tunneling picture of Schwinger pair creation and the Ritus mass shift, we consider it of the utmost importance to clarify whether the exponentiation really works beyond the two-loop level. However, a calculation of the three-loop EHL seems presently technically out of reach even for the simplest case of a self-dual field. Now, so far we have focused on spinor QED in D = 4 dimensions. However, it is important to stress at this point, that spin has not played any role in the above issue. All the statements that we made above about the weak-field expansions at one and two loops hold, qualitatively, as well for Scalar QED [13,31], starting with Weisskopf's "Scalar EHL" [32], ( 1.15) Similarly, at least at the one-loop level the structure of the EHL's and associated Schwinger exponentials at one-loop is essentially independent of the space-time dimension [33,34]. In particular, in D = 1 + 1 dimensions the EHLs become where Z = ef T , κ = m 2 /(2ef ), and the constant field strength parameter f is defined by (see app. A for our conventions). We note that the integrals can also be done in closed form, e.g., for the Spinor QED case: In 2006 Dunne and Krasnansky studied the Scalar QED EHL in various dimensions at the two-loop level [35,36] and found, in particular, the following explicit formula for this EHL in the D = 1 + 1 case: Hereα ≡ 2 e 2 πm 2 is our definition of the fine structure constant in two dimensions, and This result is formally very similar to the one above for the D = 4 EHL in a self-dual background, eq. (1.12), which suggests that QED in two dimensions (scalar or spinor) may be sufficiently similar to serve as a simpler testing ground for the exponentiation conjecture. But f1irst it had to be established whether the exponentiation conjecture itself admits a generalization to the two-dimensional case. This was settled in [37] where, using the same method of worldline instantons as Affleck et al, it was found that the exponentiation formula (1.10) generalizes to 2D Scalar QED in the form (1.22) By Borel analysis, this leads to the following formula for the limits of ratios of l -loop to one -loop coefficients: where the expansion coefficients in 2D are defined as Thus, differently from the 4D case, in 2D the asymptotic growth of the weak-field expansion coefficients increases with increasing loop level, which presumably is related to the fact that the Coulomb interaction in two dimensions is confining. In [37] it was further shown that the coefficients of the two-loop Scalar QED EHL found by Dunne and Krasnansky, eq. (1.20), indeed fulfill the asymptotic prediction (1.23) with l = 2. And there also the corresponding Spinor QED calculation was performed, leading to Here ξ 2D is the same function that appeared in the Scalar QED case, eq.(1.21), and λ 0 an infrared cutoff that affects only the irrelevant vacuum term. Note that in 2D the Spinor QED EHL has a simpler structure than the Scalar QED one, since it involves the function ξ 2D only linearly. Nevertheless, it is easy to check that the resulting weak-field expansion coefficients obey the same asymptotic relation (1.23). Thus in 2D QED we have an asymptotic prediction for the spinless case, and explicit two-loop computation confirms this prediction and suggests that it does not depend on spin. Moreover, the calculations were significantly simpler than in four dimensions, in particular not requiring mass renormalization; although in 2D QED (finite) mass renormalization exists, it turns out that, differently from the 4D case, the corresponding terms in the EHL do not contribute to the leading order growth of the weak-field expansion coefficients, and thus also not to the weak-field limit of the imaginary part. All this provides strong motivation for pushing the calculation of the EHL in 2D to higher orders. In this series of papers, we will present a complete calculation of the three-loop EHL in 2D spinor QED. Until recently, it was believed that this effective Lagrangian is (in any dimension) given at the two-loop level by the diagram shown in Fig. 3, and at the threeloop level by the three diagrams shown in Fig. 4. The electron propagators are the full ones in the external field. However, in 2016 Gies and Karbstein [38] (see also [39][40][41]) showed, that also the oneparticle reducible diagrams shown in Figs. 5 and 6 have to be taken into consideration.  Previously this type of diagram was believed to vanish, because it contains the oneloop tadpole diagram, which can be shown to formally vanish using momentum conservation and gauge invariance (see, e.g., [12,42]). However, in [38] it was shown that this is not the case any more when it appears as a subdiagram, due to the IR divergence of the photon propagator connecting it to the rest of the diagram in the zero -momentum limit. Thus we will have to include these diagrams here, too.
In this first part of our three-loop analysis we will, however, restrict ourselves to the one-fermion loop (or "quenched") three-loop diagrams, that is, the diagrams A and B of Fig. 4. These are also the ones most relevant to our main object, the verification of the exponentiation conjecture at the three-loop level.
Partial results of this calculation have been published in various conference proceedings [6,30,[43][44][45]. A first attempt at this calculation [43] used the standard approach to this type of calculation, namely Feynman diagrams with the exact electron propagator in the field. The photon propagator was taken in Feynman gauge. Due to the super-renormalizability of 2D QED, at the three-loop level the effective Lagrangian is already UV finite. However, we encountered spurious IR divergences that greatly complicated an already cumbersome calculation. In a second run [44] we used the worldline formalism [47][48][49][50] along the lines of [13,51,53,53,54], and also encountered IR divergences. Those could be removed by suitable integrations by parts, but it then was found that it is also possible to avoid the appearance of IR divergences altogether by the choice of a particular covariant gauge. This is the gauge ξ = −1, to be called "traceless gauge" in the following, since it makes the photon propagator traceless in D = 2 (it had already been used in [37], but only in the worldline instanton calculation). Returning then to the Feynman diagram calculation, we found that here, too, this gauge makes the IR finiteness manifest, and moreover leads to very significant simplifications due to the identities (3.3), (3.4) below.
Thus both the Feynman diagram and the worldline calculation in this gauge yielded integral representations for diagrams A and B that are manifestly finite term by term. However those representations are quite different: the one from Feynman diagrams is much more compact, while the one from the worldline formalism has the advantage of combining the contributions of the planar and the non-planar diagram, something that would be hard to achieve in the diagrammatic approach. Thus we have chosen to present here both representations, for whatever their worth may be. Since our specific purpose will require a fairly high-order computation of the weak-field expansion coefficients of this Lagrangian, we further show how these coefficients can be obtained from the Lagrangian in an efficient way. Here we use the representation from the diagrammatic approach. As expected on general grounds, the coefficients are of the form r 1 + r 2 ζ 3 with rational numbers r 1 , r 2 , where the ζ 3 comes from the non-planar diagram only. We compute the first two coefficients analytically, and four more by numerical integration.
The organisation of this paper is as follows. In section 2 we calculate the three-loop quenched EHL in the worldline formalism, including also a recalculation of the two-loop EHL. The two-loop EHL had been obtained already in [44] using the Feynman diagram approach, but we include the worldline calculation here because it displays some interesting simplifications. In section 3 we calculate the three-loop quenched EHL again in the diagrammatic approach. Based on the resulting four-parameter integral representation, we present algorithms for the computation of the weak-field expansion coefficients in section 4. Here for the non-planar contribution we make essential use of the high symmetry of diagram B, in two ways: first to develop a certain integration-by-parts procedure, and second to rewrite the integrand in Schwinger parameter space in terms of polynomial invariants of the dihedral group D 4 , which helps greatly to keep the expressions manageable. In section -10 -5 we give the results of a calculation of the first six coefficients for both diagrams. Here for the planar diagram A we have analytic results for all six coefficients, while for the incomparably more difficult non-planar diagram B we have computed the first two coefficients analytically, the remaining four numerically. Section 6 gives our summary and outlook. There are four appendices: appendix A gives our conventions. In appendix B we list the momentum integrals appearing in the three-loop worldline calculation. Appendix C contains the corresponding list for the Feynman diagram calculation. Finally, in appendix D we present some elements of the invariant theory of the dihedral group D 4 , and sketch the derivation of the basis of invariants used in the weak-field expansion algorithm of diagram B in section 4.
2 Worldline calculation of the quenched three-loop EHL

Worldline representation of 2-photon and 4-photon amplitudes in a constant field
The starting point for our calculation of the two-loop and three-loop Euler-Heisenberg Lagrangians in the worldline formalism are the following representations of the 2-photon and 4-photon amplitudes in a constant field (see [53,54] for context and derivation of these formulas). Define the field strength tensor by and further Z ≡ ef T . Introduce the vacuum worldline Green's functionṡ and the generalized (constant field) Green's functions with the scalar, dimensionless coefficient functions We note the symmetry properties Denote further by the field strength tensor associated to photon i, and define the "super-bicycle of length n" G S (i 1 i 2 · · · i n ), involving a subset i 1 , i 2 , . . . , i n of the integration variables, bẏ Then, the two-photon amplitude in the constant field can be written as where D is the space-time dimension. We will set D = 2 from the beginning, since dimensional regularization will not become necessary in our calculation.
-12 -Further, let us define the "one-tail" T (i) and "two-tail" T (ij) by The four-photon amplitude in the constant field can then be written as [54] (omitting the global factor (2π where the upper index on a Q i denotes the "cycle content": Let us mention that each of the sixteen terms appearing in the decomposition (2.11) gives a contribution to the four photon amplitude that is separately gauge invariant [53][54][55]. The representations (2.8), (2.10) hold mutatis mutandis also in 4D QED, however in the 2D case they are particularly useful, because here all matrices appearing in the above expressions (the F i 's and all worldline Green's functions) commute with each other. This is because in two dimensions all antisymmetric matrices are multiples of each other, and the Green's functions involve only the matrix F . Thus any matrix, or product of matrices, appearing in our calculations below is of the form M = a1l + bE, which also entails the useful identities These simple observations will lead to many simplifications in the following calculations.
In particular, the commutativity implies that, for even n, the super-bicycles (2.7) factorize asĠ 14) The representations (2.8), (2.10) hold off-shell, so that the one-loop amplitudes can be used to construct (quenched) higher-loop photon amplitudes by sewing off pairs of photons, say, the photon legs with index i and j: Using an arbitrary covariant gauge with gauge parameter ξ, this is implemented by setting and adding the integration d 2 k i /(4π 2 k 2 i ) (there is also a combinatorial factor of 1 2 for each pair of legs).
Thus in the following we will construct the two-loop Euler-Heisenberg Lagrangian from the one-loop vacuum polarization tensor, and the three-loop EHL from the one-loop fourphoton amplitude. In the worldline formalism, one could construct these higher loop Euler-Heisenberg Lagrangians also directly using the concept of multi-loop worldline Green's functions [53,56,57]; however, this would obscure the decomposition (2.11), which we will find very useful in the following. This is because, first, it will allow us to substantially reduce the number of terms in the integrand using the symmetries of the problem; and second, because the gauge invariance term-by-term of this decomposition implies that for each term we can use a different gauge parameter ξ in the sewing procedure (2.15).

One-loop vacuum polarization tensor in a constant field
Let us start with calculating the vacuum polarization tensor, related to the two-photon amplitude (2.8) by (2.16) We can simplify this calculation by observing that the constant-field vacuum polarization tensor carries the transversal projector familiar from the vacuum case, We note that this is another peculiarity of the 2D case, and due to the fact that no other symmetric and transversal tensor of second degree can be built from δ µν , k µ , and F µν . In 4D QED the tensor Π µν in a constant field involves several independent transversal structures (see, e.g., [2]). Thus Π µ µ (k) = k 2 Π(k), and it is sufficient to calculate the trace of Π µν , corresponding to the replacement ε µ 1 ε ν 2 → δ µν . Making this replacement in F 1,2 , setting k = k 1 = −k 2 , and using (2.7), we are led to compute In the exponent we get (2.20) Putting things together, (2.21) Finally we use, as usual in this type of calculations [53], the translation invariance of the worldline Green's functions to put τ 2 = 0, τ 1 = τ , and rescale τ = T u. With a further change of variables from T to Z and from m 2 to κ = m 2 /(2ef ) we obtain our final form for the vacuum polarization tensor,

Two-loop EH Lagrangian
We now construct the two-loop EH Lagrangian from the polarization tensor. Using (2.22) in the sewing procedure gives (note that cosh Z ≥ cosh(Ġ B12 Z)) (2.23) Remarkably, after the gaussian k -integration all u -dependence has already cancelled out, and one is left with a one-parameter-integral: (2.24) We note that the lowest order (f -independent = vacuum) term has an UV divergence at Z = 0. Subtracting this lowest order term, that is replacing 1/ sinh 2 Z by 1/ sinh 2 Z −1/Z 2 , we can apply the standard integral formula [23] (2.25) Thus we find

Three-loop quenched EH Lagrangian
We proceed to the construction of the quenched three-loop EH Lagrangians, starting from the worldline representation (2.10) of the one-loop four-photon amplitudes in the field. We -16 -apply the sewing procedure (2.15) with k 1 = −k 2 =: k and k 3 = −k 4 =: l. The universal exponential factor in (2.10) can then be written as We will also need the determinant and the inverse of M, In the following we will often treat as numbers matrices which are proportional to the unit matrix, such asĜ B12ĜB21 .
-17 -Further, we observe that, after the sewing, the resulting total worldline integral is invariant under the exchanges τ 1 ↔ τ 2 , τ 3 ↔ τ 4 , and (τ 1 , τ 2 ) ↔ (τ 3 , τ 4 ). This implies that, of the sixteen terms appearing in the decomposition (2.11), only seven are independent and need to be computed; namely, we can now make the replacements We will now first work out the effect of the sewing replacements on these seven structures. This is simplest for the ones in Q 4 and Q 22 ; e.g., inĠ S (1324) we can use the commutativity of all matrices to write In this way, we obtain (independently of ξ), where we have further introduced the abbreviation (2.37) Proceeding to Q 3 , here we also use E 2 = 1l to rewrite After the sewing of ε 3 with ε 4 this giveṡ (the term involving ξ drops out here because it contains a factor lEl = 0).
Coming to the terms in Q 2 , here things are more delicate. It turns out that, to avoid the appearance of spurious infrared divergences in the k, l integrations, one should use the gauge ξ = −1. This leads tȯ (2.42) The momentum integrals can now be easily performed, without encountering any IR divergences, making them gaussian by exponentiating each denominator, The parameter integrals over x, y are convergent since α, β, ∆ > 0. A list of all integrals is given in appendix B. The exception is the first term in the square brackets in (2.40) which contains an IR divergence in the l integral for any gauge except for our choice ξ = −1; here care must be taken, and we therefore present its calculation in some detail. We abbreviate A =Ĝ B23 −Ĝ B13 and B =Ĝ B41 −Ĝ B42 . The 1/k 2 from the propagator here cancels against the k 2 in the prefactor of (2.40). Thus the k integral is finite, and performing it first gives (2.44) The first term in the brackets, which is the IR-critical one, vanishes since tr(AP 34 B) = tr(P 34 BA), BA is a linear combination of 1l and E, and trP 34 = tr(P 34 E) = 0. In the second term, we can use the identity (2.13) to write Thus also the l -propagator cancels, and doing the l -integral we obtain for the total integral the result (2.46) All the other integrals are straightforward. In this way we arrive at our final result for the quenched part of the three-loop Lagrangian, corresponding to the sum of the diagrams A and B of Fig. 4: Note that in the calculation of I 13 we have used the identities (2.12), (2.13), the fact that all matrices are commuting, and thatĜ Bij +Ĝ T Bij ∼ 1l. Note also that our final integrand is still given in the decomposition corresponding to (2.10),(2.33), with the subscript on a term denoting its cycle content.

Feynman diagram calculation of the quenched three-loop EHL
In this section, we will calculate the quenched part of the three-loop EHL another time, now using the standard formalism. In terms of Feynman diagrams, the quenched part is given by diagrams A and B of figure 4.

Definitions
We use the 2D electron propagator in a field with constant field strength tensor F = 0 f −f 0 , where f > 0. The electron propagator in the field can, using the proper-time representation and Fock-Schwinger gauge, be written as (see appendix A for our conventions) ). In the following we will often omit the arguments in ϕ(z, p) and g(z, p) and replace them by superscripts, e.g.,ĝ ≡ g(ẑ,p). Moreover, we will abbreviate γ ≡ 1/ cosh z, γ ≡ 1/ cosh z , etc. Greek indices run over the values 1,2 as usual and we use the shorthand x ∧ y ≡ x α y β αβ . As in our three-loop calculation in the worldline formalism of the previous section, in the Feynman diagram calculation, too, at the three-loop level it turns out to be essential for avoiding spurious IR divergences to take both photonď propagators in the "traceless" gauge ξ = −1, where In this gauge one has, in D = 2 and for any n ≥ 0, the extremely useful identities

Diagram A
We start with the simpler one, diagram A. We parametrize this diagram as in figure 7, with independent momenta p, q, k, so that the kinematic relations are The only contribution from the double trace within (5) will be (odd powers of m never appear by parity, Euclidean QED so the amplitude is real) the real and even (in m) part of: the terms in the trace that contribute are two, namely, mass dependent and mass independent: (1) = m 2 tr(e σ 3 (z−ẑ) k /p / k /q / p / q / ) (2) = −γγtr(p / k /p / k /p / q / p / q / ) Also, we write the p-dependent part of the argument of the exponential by defining ϕφφϕ := Ωe −ap 2 +p·β e −(tk 2 +t q 2 ) , where we abbreviate t ≡ tanh z/ef , t ≡ tanh z /ef , etc., defined a = t +t +t + t , Ω = γγγγ (ef ) 4 e −2κ(z+ẑ+z+z ) , and the vector β = −2(tk + t q).
The amplitude is then given by where dz is an obvious shorthand. We find both contributions separately The Feynman rules give for this diagram the following expression: (3.6) Note that this diagram has to be taken twice, and we include this factor in its definition. Applying the identities (3.3), (3.4) (with n = 0) inside the Dirac trace, we can rewrite More explicitly, the trace can be written as tr(gk /p / k /ĝq / p / q / ) = I 1 + I 2 (3.8) where I 1 = m 2 tr(e σ 3 (z−ẑ) k /p / k /q / p / q / ) Further, we rewrite the p-dependent part of the argument of the exponential by defining ϕφφϕ ≡ Ωe −ap 2 −2(tk+t q)·p e −(tk 2 +t q 2 ) , (3.10) where we abbreviated t ≡ tanh z/ef , t ≡ tanh z /ef , etc., and defined The amplitude is then given by where Dz here and in the following is a shorthand for the fourfold z -integral. Working out the traces (3.9) and performing the gaussian p -integration, one finds Here J 1,2 still correspond to I 1 , I 2 , and we have introduced the further abbreviations (3.16) Note that d A ≥ 0, e A ≥ 0, as required for convergence, and also that d A e A − g 2 A ≥ 0.
-24 - We now come to the integrals over the photon momenta k, q, and considering the factors of 1 k 4 q 4 in (3.14), (3.15), it is not obvious that they are IR finite. As was already mentioned, indeed in a general covariant gauge here one would encounter spurious IR divergences. However, this happens not to be the case in traceless gauge; the k, q integrals of J 1,2 are completely finite, and moreover elementary. In appendix C we list all the k, q -integrals needed to complete our calculation of diagram A, as well as of diagram B below.
After these integrations, considerable simplifications occur, leading to the following simple results for J 1,2 : Adding up both contributions and taking prefactors into account gives finally the integral representation: where we have introduced A ≡ tanh z + tanh z + tanhẑ + tanhz .

Diagram B
We come to diagram B, which is expected to be more difficult. See figure 8 for our parametrization. We again use k, q, and p as the independent variables. The remaining variables are expressed in terms of them as We will use k, q, and p as the independent variables. The remaining variables are expressed in terms of them as With these conventions, the contribution of this diagram is written as With these conventions, the contribution of this diagram is written as The Dirac trace in (3.21) can again be simplified using (3.3), (3.4). The result can be written as Here the term in b 4 involving sinh can be dropped since it is odd under the exchange (k ↔ q, z ↔z). For diagram B it turns out to be convenient to first perform the gaussian p integral, without even working out the trace appearing in b 0 . This integral appears in the form where only b 0 involves terms with numerators. One can then rewrite the exponent in (3.24) as with a, t,t . . . defined as for diagram A, and b = (t +t)k + (t + t )q c =tk 2 + t q 2 +t(k + q) 2

(3.26)
Any Kronecker -δ coming out of the p integration will belong to b 0 , and it will make the trace in b 0 vanish on account of the first of the identities (A.10). Therefore the effect of the p integration can also be described as a replacement Thus for the b 0 term we find (3.28) We further abbreviate and work out the trace: (3.30) The remaining k, q integrals are again finite and elementary. We rewrite the new exponent in (3.28) as The five integrals arising from (3.30) are included in appendix C. After simple manipulations, one finds the following result for the total momentum integral of b 0 : For b 4 we need only a single integral, given in (C.1a), and leading to t + t +t +t × 1 + (t + t +t +t)(tt t + tt t + ttt + t tt ) (tt − t t ) 2 ln (t + t +t +t)(tt t + tt t + ttt + t tt ) (t + t )(t +t)(t + t )(t +t) (3.34) Putting things together, our final result for diagram B becomes where B = (tanh 2 z + tanh 2ẑ )(tanh z + tanhz) + (tanh 2 z + tanh 2z )(tanh z + tanhẑ) C = tanh z tanh z tanhẑ + tanh z tanh z tanhz + tanh z tanhẑ tanhz + tanh z tanhẑ tanhz G = tanh z tanhẑ − tanh z tanhz

Calculation of the weak-field expansion coefficients
As was explained in the introduction, our main motivation for this calculation of the EHL in 2D is in the properties of the weak-field expansion. For the one-loop and two-loop contributions, the expansion coefficients have already been obtained in closed form in [37].
In this chapter, we will develop algorithms for the extraction of the three-loop expansion coefficients from the integral representations obtained in the previous two chapters.
In (1.24) we defined the weak field expansion coefficients at l loops by Thus at three loops, and now introducing, instead of κ, the expansion variable ρ ≡ 1 2κ = ef m 2 , we can write It will be convenient to replace the coefficients c and to denote the contributions of the two diagrams to Γ n by Γ A,B n , of their sum by Γ A+B n .

Diagram A
The case of diagram A is straightforward. Starting from our final representation (3.18), we rescale z = ef m 2 w = ρw etc. Then we can write where a ≡ w + w +ŵ +w . (4.4) We expand the integrand I A in powers of ρ, the coefficients α n of this expansion can be decomposed into terms of the form with k ≥ 2 and polynomials P nk . After exponentiating the denominator using a Feynman parameter, the integrand factorizes in integrals over w, w ,ŵ,w that are elementary. The integrand of the final λ -integral consists of terms of the form with integer numbers m ≥ 2 and polynomial numerators P (λ). Thus the λ -integral is elementary, too. This also makes it clear that the coefficients Γ A n are all rational numbers. For example, at lowest order we find simply and to Γ A 0 = − 1 3 . Using this method, we have found it easy to obtain the first twelve expansion coefficients for diagram A.

Easy part of diagram B
We proceed to the much more difficult case of the non-planar diagram B. Rewriting again z = ef m 2 w = ρw etc. (3.35) becomes where with a corresponding rewriting of the functions A, B, C, G defined in (3.19), (3.36), a defined as in (4.4) andw ≡ w − w +ŵ −w.
We will treat I easy with c ≡ ww ŵ + ww w + wŵw + w ŵw , (4.16) k ≥ 2 and polynomials Q kl n . The factor 1 a k is exponentiated as in (4.8). In the factor 1 c l we first rewrite -31 -and then exponentiate only the second factor: After this, the integrand factorizes in integrals over w, w ,ŵ,w, and these integrals are of the form One is left with a double integral over λ and µ, involving a product of four Bessel functions and rational factors. Remarkably, MATHEMATICA has an algorithm for the analytic evaluation of this type of integral [58]. The result is of the form with rational numbers r 1,2 . At lowest order one has leading to At the next order we find polynomials Q 31 1 , Q 41 1 , Q 32 1 , where, for example, (4.23) The full calculation is already too long to be presented here, so let us jump to the final result:

Hard part of diagram B
Finding a method of closed-form evaluation for I hard B is more difficult. One of the difficulties with its integration is that straightforward attempts will create spurious divergences at G = 0. For obtaining a first integral we will make essential use of the high symmetry of this diagram. Namely, while diagram A above has only a Z 2 × Z 2 invariance generated by the interchanges w ↔ŵ and w ↔w, diagram B additionally to those has the symmetry under pair exchange (w,ŵ) ↔ (w ,w). The group generated by these three reflections is the eight-element dihedral group D 4 . The expansion of I hard yields coefficients that can be decomposed in the following way: let us introduce, besides a and c, now also Let us further introduce the basis functions Then we can, using ac = h − g 2 , write where the u k,l n and v k n are polynomial functions of the w, w ,ŵ,w. Here in general the index kcan take both even and odd values. However, for reasons that will become clear below we shall eliminate the case of odd k, multiplying by a factor of g in the numerator and the denominator wherever necessary.
Further, the symmetries of diagram B, together with the natural appearance of the variablew, motivate us to introduce a differential operatord, -33 -which is in some sense conjugate tow. Its action onw, a, c, g, h is simple: Given that h is constant under the action ofd, as far asd is concerned U k,l and V k are functions of the single variable g. Moreover, they are functions simple enough to be integrated in closed form, for any k, l and an arbitrary number of times. Given that alsod 2 g =da = 0, we can use integration-by-parts with respect to the variable g to write any term of the form appearing in (4.28) as a total derivative with respect tod, simply by repeatedly integrating the factors U k,l and V k and differentiating the polynomials u k,l n , v k n until the latter vanish. Denoting by deg the degree as a polynomial in w, w ,ŵ,w, and by, e.g., U (−i) k,l the i -fold indefinite integral of U k,l in g, we can then explicitly write β hard n as a total derivative with respect tod: β hard n =dθ n , where Thus we have now reduced the integrand to boundary terms: After this rescaling, θ n |w =0 depends on T only through a global factor of T 2n : θ n |w =0 = T 2n ω n (α, α ,α) . We use the delta function to eliminate α , rather than α, orα, since by the symmetry of diagram B the resulting integrand will be symmetric in α,α. where we have renamed f n (α,α) ≡ ω n (α, 1 − α −α,α). We note that at this stage we have the correspondences For the computation of the final integral over α andα in (4.35) it will be advantageous to perform a change of variables fromα to Λ,  Let us carry through this algorithm for the lowest coefficient n = 0. Taking the ρ → 0 limit in (4.13) and eliminating c through c = (h − g 2 )/a yields For constructing θ 0 via (4.31) we need Going through the steps leading to (4.35), we get with g, h as in (4.36). After the change of variables (4.37) this becomes Both integrals can be done in closed form, and one finds (4.46) In principle, this algorithm after computerization can be used to calculate the weak-field expansion coefficients of I hard B to arbitrary order. The problem is that the polynomials u k,l n and v k n rapidly grow in size with increasing n. However, since the complete integrand is invariant under the action of the group D 4 , and both h and g 2 possess this invariance, all those polynomials must be invariant, too (note that this would not be quite true had we permitted odd powers of g in the denominator, since g itself is only a semi-invariant: it is invariant under w ↔ŵ and w ↔w, but changes sign under the pair exchange (w,ŵ) ↔ (w ,w)). This suggests to improve the algorithm by using the representation theory of the group D 4 to rewrite those polynomials more compactly in terms of some basic polynomial invariants of D 4 .
As we explain in detail in appendix D, the application of the general theory of polynomial representations of finite groups (see, e.g., [59]) to our case of the group D 4 , acting on polynomials of four variables, shows that all our numerator polynomials can be rewritten in the form where P (a,w, v, j) is a polynomial in the four invariants (more precisely semi-invariants, see appendix D) a,w, v, j. Of those a,w we have already seen, and the remaining two are defined by v = (w +ŵ)(w +w) + 2(wŵ + w w) , This choice of a basis of invariants is well-adapted to our integration-by-parts procedure, since but forw all get annihilated byd.
We are now ready to tackle the n = 1 case.   As a final remark, at high orders it might become useful to give a separate treatment to the denominator cosh(ρw) cosh(ρw ) cosh(ρŵ) cosh(ρw) since it has full S 4 symmetry, while the rest has only D 4 symmetry; We remark that this invariant-based approach could as well be applied to I easy B , but it would require a higher-order calculation to see whether this would be more efficient than the method described above.

Low-order results
In Table 1 we give the first six coefficients for both diagrams A and B. The coefficients for A and the first two coefficients for B were calculated analytically using the methods developed in the previous section. The remaining coefficients for B were obtained by numerical integration. Note that the coefficients for diagram A are rational, for B of the form r 1 +r 2 ζ 3 with rational numbers r 1 , r 2 .
To the best of our knowledge, there are no results available in the literature that could be used to perform some check on our results. However, let us mention that we have an internal check for the first coefficient, since for this one we had obtained an analytical value already in a previous calculation that used Feynman gauge [43]. In Feynman gauge, we had found

(5.2)
We see that the sum agrees, Γ A+B 0 = − 3 2 + 7 4 ζ 3 , and that the rational part becomes gaugeindependent only in the sum over both diagrams, while the ζ 3 -part can come only from the non-planar diagram in any gauge.
Computing a number of coefficients sufficient to address the issues related to the exponentiation conjecture will require a more substantial computational effort, which we leave for future work.

Conclusions and Outlook
We have presented here the calculation of the three-loop correction to the Euler-Heisenberg Lagrangian in 1+1 dimensional massive spinor QED. The calculation has been performed in parallel using standard Feynman diagrams and the worldline formalism, treating the constant external field non-perturbatively in both cases. In both formalisms, the use of the "traceless" gauge ξ = −1 for the internal photons turned out beneficial in making the IR finiteness of the effective Lagrangian manifest term-by-term. In the Feynman diagram approach, this gauge choice moreover led to substantive simplifications. Both methods -39 -led to four-parameter integral representations, although of a quite different structure. The Feynman parameter calculation results in relatively compact integrands, particularly for the non-planar diagram A, while the worldline formalism yields a more extensive integrand, but has the advantage that this integrand holds for both the planar and the non-planar sector. We have further used the worldline formalism for a recalculation of the two-loop EHL that is simpler than the Feynman diagram calculation of [37].
To the best of our knowledge, this constitutes the first calculation of a three-loop effective Lagrangian in quantum electrodynamics, and also the first three-loop calculation in 1+1 QED (in fact even at the two-loop level the only calculations in 1+1 QED that we have been able to find in the literature are the already mentioned studies by [35,36] for the Scalar QED case, and the recent [60] on super QED (more effort seems to have gone into exploring the expansion in the mass, see [61], [62] and refs. therein).
Based on the representation obtained in the Feynman diagram approach, we have developed computerizable algorithms for the analytic calculation of the weak-field expansion coefficients that, in principle, work to arbitrary orders. For diagram B, our algorithm make use of the invariant theory of the symmetry group of the graph, the dihedral group D 4 , in a way that may be generalizable to other multiloop graphs and thus of independent interest. As to explicit computation, here we have been satisfied with a low-order test, leaving to future work the task of obtaining a sufficient number of coefficients to confirm or refute the exponentiation conjecture, which is our main motivation for pushing this computation to the three-loop level.
Moreover, the methods developed here should also become useful in an eventual calculation of the three-loop EHL in four dimensions, particularly for the self-dual case which is to some extent a "doubling up" of the two-dimensional case.

C Momentum integrals appearing in the calculation of diagrams A and B
Here we list all the integrals over the photon momenta k, q appearing in the calculations of diagrams A and B in section 3. The parameters d, e, f should be replaced by the ones defined in (3.16) for diagram A, and defined in (3.31) for diagram B. d 2 k k 4 d 2 q q 4 e −dk 2 −eq 2 +2gk·q [2(k · q) 2 − k 2 q 2 ] = π 2 1 + de g 2 − 1 ln 1 − d 2 q q 4 k 4 q 2 k · q e −dk 2 −eq 2 +2gk·q = π 2 g d 2 e D Invariants of the dihedral group D 4 Given a finite group G and a real n−dimensional representation (the representation can in fact be also complex) ρ(G) = Γ ⊂ GL(n, R). The representation ρ may or may not be irreducible. The representation ρ induces naturally an action on the set of polynonials with n variables R[x 1 , · · · , x n ]. Denoting X = (x 1 , · · · , x n ) t an n−dimensional vector of R n , the action of g ∈ Γ on X is denoted g · X. A polynomial I ∈ R[x 1 , · · · , x n ] is said to be invariant if for any g ∈ Γ we have g · I(X) ≡ I(g · X) = I(X) . (D.1) The set of polynomial invariants is denoted R[x 1 , · · · , x n ] Γ = I ∈ R[x 1 , · · · , x n ] s.t. g · I(X) = I(X) ∀g ∈ Γ .

(D.3)
All polynomial invariants can be expressed in terms of what is usually called primitive and secondary invariants. The number of primitive invariants is equal to the dimension of the representation space, i.e., is equal to n [59,63], so we denote them by P 1 , · · · , P n . They are algebraically independent, which means that they do not satisfy any polynomial identity of the type f (P 1 , · · · , P n ) = 0. Denote now by R = R[P 1 , · · · , P n ] the subalgebra of polynomial invariants generated by the primitive invariants. Equivalently R can be seen as the set of polynomials in P 1 , · · · , P n . Denote d 1 , · · · , d n the degree of P 1 , · · · , P n respectively. The number of secondary polynomial invariants is equal to m = d 1 · · · d n /|G| [64]. Denote now S 1 , · · · , S m the set of secondary polynomial invariants. The set of all invariants (P 1 , · · · , P n , S 1 , · · · , S m ) are not any more algebraically independent and consequently do satisfy algebraic relations which are called syzygies. It turns out that the subalgebra of invariants R[x 1 , · · · , x n ] Γ is a free R−module with basis (S 1 , · · · , S m ). In particular this means that any invariant I ∈ R[x 1 , · · · , x n ] Γ can be uniquely written as where f i (P 1 , · · · , P n ), i = 1, · · · , m belongs to R, i.e., are polynomials in (P 1 , · · · , P n ).
There exist many automated way to compute primary and secondary invariants and we have used the Computer Algebra System for Polynomial Computations called SINGU-LAR [64]. Before considering our specific case recall that the polynomial invariants for the n−dimensional representation of the group of permutation Σ n acting on n elements consist of the well-known symmetric polynomials. This in particular means that the set of symmetric polynomials of degree d = 1, · · · , n is the set of primitive invariants and the only secondary invariant is σ 0 = 1. In particular taking the notations of the paper we have n = 4 and x 1 = w, x 2 =ŵ, x 3 = w , x 4 =w and the symmetric polynomials in this case are given by σ 1 = w +ŵ + w +w , σ 2 = wŵ + ww + ww +ŵw +ŵw + w w , (D.5) σ 3 = wŵw + wŵw + ww w + w ŵw , The dihedral group D 4 is eight dimensional and the three generators in the four dimensional (reducible) representation are given by their action on (w,ŵ, w ,w) g 1 : w ↔ŵ , g 2 : w ↔w , (D.6) g 3 : (w,ŵ) ↔ (w ,w) .
For instance this means that there are 16 linearly independent polynomial invariants of degree six. Since we are considering a representation of dimension four of D 4 we must find four primitive invariants. Using SINGULAR, we obtain the primitive invariants a = σ 1 = w +ŵ + w +w , λ = (w +ŵ)(w +w) , Consequently any invariant can be uniquely written in the form I = f 1 (a, λ, µ, σ 4 ) + f 2 (a, λ, µ, σ 4 )c . (D.10)