Three-loop Euler-Heisenberg Lagrangian in 1+1 QED. Part I. 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 expansion 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ζ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:
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]: (훽 = 푒퐸/푚 2 ). In the following we will focus on the weakfield limit 훽 ≪ 1 where only the first of these "Schwingerexponentials" is relevant.
The nonperturbative dependence of the Schwingerexponentials on the field supports the interpretation of field-induced pair creation as a vacuum tunneling effect, as proposed by Sauter as early as 1931 [18].
As usual in quantum field theory, the real and imaginary parts of the EHL are related by a dispersion relation. For the 푁-photon amplitudes at full momentum, this would be a standard dispersion relation performed diagram-by-diagram, relating the diagrams of Figure 1 to the "cut diagrams" shown in Figure 2, involving on-shell electrons.
However, in the zero-energy limit the cut diagrams all vanish, since a finite number of zero-energy photons cannot create a pair on-shell. Thus what counts here is only the asymptotic behaviour for a large number of photons, and instead of an ordinary dispersion relation we have to use a Borel dispersion relation. This works in the following way [7]. Consider the purely magnet powers of the field yields with an effective expansion par coefficients 푐 (1) 푛 that can be writt numbers 퐵 푛 : Here 푐 (1) 푛 holds information on tudes. The asymptotic behaviou easily studied using well-known numbers. One finds Thanks to the factor (−1) 푛 , the in hand side of (16) all give conv is one (rather roundabout) way magnetic EHL has no imaginary to pair creation. The analogous expansion for is almost the same: where now 푔 = (푒퐸/푚 2 ) 2 , but w the additional factor (−1) 푛 make which is crucial, because now th (16) leads to divergent Borel integrals do, however, all poss parts, by a (now ordinary) dis 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 figure 1, where all photon energies are small compared to the electron mass, ω i m. 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 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 the important difference that no particular ordering of the photon legs along the loop needs to be fixed. Thus the term 푆 ext alone upon expansion yields the diagrams of Figure 1 (where each leg now stands for an interaction with the arbitrary field 퐴(푥)).
The "worldline instanton" of Affleck et al. [26] is an extremal trajectory of the worldline path integral for a stationary phase approximation electric field in the 푧 direction th is given by a circle in the (Euclid Thus the contribution of 푆 0 + 푆 e (one-loop) Schwinger-exponent the e 훼휋 factor. Thus Affleck et al. arrive, wit exponentiation for Scalar QED t Spinor QED: Their argument assumes the fiel restriction on the strength of t following: (i) Formula (30), if true, c all-loop summation of of arbitrary loop order. Table 1.
(iii) Perhaps, most surprising in (30) is already the implying that the worldl matically takes all mass grams into account. Th that the determination o for the EHL becomes a r at two loops [20,22,23].
Thus, according to Affleck et has produced the factor e 훼휋 , wh simple but also perfectly anal constant 훼!. According to wha 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 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, 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 figure 2 (here it is understood that internal photon corrections are put in all possible ways).
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].

JHEP03(2019)167
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 expansion [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 mass-shift,

JHEP03(2019)167
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 electronpositron 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 -5 -

JHEP03(2019)167
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 would be a formidable challenge 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], 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 appendix 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 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 first 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 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 figure 3, and at the threeloop level by the three diagrams shown in figure 4. The electron propagators are the full ones in the external field.
Previously this type of diagram was believed to vanish, because it contains the one-loop 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 -8 -

JHEP03(2019)167
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 figure 4. These are also the ones most relevant to our main objective, the verification of the exponentiation conjecture at the three-loop level. In the analysis of [19] only quenched diagrams contribute to the exponentiation formula; diagrams with more than one fermion loop do not survive in the weak-field/large N limit. For our three-loop diagram C in D = 2, we have verified already in previous work that this is indeed the case (see figure 8 in [30]). For the "new" tadpole diagrams, the same can be easily seen from the Gies-Karbstein factorization formulas [38,39], as we will show in part 2.
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 superrenormalizability 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,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 of diagrams A and B manifest, and moreover it 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 -9 -

JHEP03(2019)167
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 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 five 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. 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. Finally, appendix E provides a bibliographical overview over the rather dispersed results that exist on multiloop Euler-Heisenberg Lagrangians.

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. Further, let us define the "one-tail" T (i) and "two-tail" T (ij) by -11 -

JHEP03(2019)167
The four-photon amplitude in the constant field can then be written as [54] (omitting the global factor (2π) 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 aṡ 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 indices i and j: using an arbitrary covariant gauge with gauge parameter ξ, this is implemented by setting

JHEP03(2019)167
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 Putting things together, 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, where nowĠ B12 = 1 − 2u.

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)) Remarkably, after the gaussian k-integration all u-dependence has already cancelled out, and one is left with a one-parameter-integral: 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] 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 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 . 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

JHEP03(2019)167
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 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 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 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 threeloop level it turns out to be essential for avoiding spurious IR divergences to take both photon propagators in the "traceless" gauge ξ = −1, where 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 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

Term (1) p-integral
We decompose further A straightforward calculation gives after some work 2 Figure 7. Parametrization of diagram A.
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 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) -20 -

(3.11)
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. 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, qintegrals 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 : 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 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 With these conventions, the contribution of this diagram is written as The Dirac trace in (3.21) can again be simplified using 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 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 and work out the trace:

JHEP03(2019)167 4 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 section, we will develop algorithms for the extraction of the three-loop expansion coefficients from the integral representations obtained in the previous two sections.
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. where the coefficients β easy n can be decomposed into terms of the form 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 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

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 k can 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.

JHEP03(2019)167
Further, the symmetries of diagram B, together with the natural appearance of the variablew, motivate us to introduce a differential operatord, 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 alsõ d 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 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 (α, α ,α) .

JHEP03(2019)167
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 Λ, α(α, Λ)) . (4.38) 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

JHEP03(2019)167
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 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 P (a,w, v, j) a 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.

JHEP03(2019)167
We are now ready to tackle the n = 1 case. For β hard 1 one finds β hard 1 = u 00 1 + u 20 with coefficient functions that in our new basis read It 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  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 led to four-parameter integral representations, although of a quite different structure. The Feynman parameter calculation results in relatively compact integrands, particularly for JHEP03(2019)167 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 is 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.
We find that the coefficients are of the form r 1 + r 2 ζ 3 with rational numbers r 1 , r 2 , where the ζ 3 comes purely from the non-planar diagram B. The appearance of ζ 3 at the three-loop level is generic, so that it would have been not much short of a miracle to find our three-loop coefficients to be rational. However, cancellations of odd ζ n values abound in QED, particularly in the quenched sector (see, e.g., [63,64]), and one could have expected such a miracle to happen for two reasons: first, the fact that the two-loop result (1.25) is (up to a constant) just a second derivative of the one-loop one (1.19), which conceivably could have been the beginning of an exponentiation that would have implied rationality to all loop orders. Second, the already mentioned work by Dunne and Krasnansky [3,35,36,65] where they develop, for the case of scalar QED in a self-dual background in any (even) dimension, a "loopology" algorithm that allowed them, using only integration-by-parts and algebraic manipulations, to express the two-loop EHL in terms of the one-loop one, thereby reproducing the scalar QED analogues of the two-loop formulas (1.25) and (1.12) without having to actually confront any two-loop integrals. Vice versa, we can predict from the appearance of ζ 3 in our calculation that their algorithm must run into some obstacle at the three-loop level, since otherwise a complete reduction to one-loop and two-loop diagrams would be possible, and those can be expressed in terms of the diagamma function and its derivatives, whose large-argument expansion (which is the relevant one for the weakfield expansion) contains only even ζ n values. Presumably this obstacle is due to the first appearance of a nonplanar diagram, and it would be very interesting to pinpoint what it is (first steps at the three-loop level were already taken in [66,67]).
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 (in the sense that in the -34 -

JHEP03(2019)167
self-dual case the field strength tensor consists of two copies of (1.18) and fulfills the same identity F 2 = −f 2 1l).
The special utility of the gauge ξ = −1 in D = 2 calculations also deserves further investigation. Although this gauge is known from four dimensions to have some special properties [68], possibly related to the fact that it makes the photon propagator proportional to the conformal tensor, 1 we would like to be able to predict whether higher-loop corrections to the EHL in D = 2 are still manifestly IR finite in this gauge. For this purpose, it should be useful to study more generally the fermion self-energy diagrams at various loop orders in the Schwinger model (for the one-loop fermion self energy, which is purely IR divergent in two dimensions, it was already shown in [69] that it becomes finite in the ξ = −1 gauge).

JHEP03(2019)167 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.
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) . The number of linearly independent polynomial invariants of degree k is given by the k-th coefficient of the Molien series, around zero, of the generating function [59] 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,70], 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| [71]. 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 [71]. 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) .
We will denote such a set by (1 : a, 2 : λ, 2 : µ, 4 : σ 4 ; 3 : c), where the first entry denotes the degree of the polynomial and the semi-colon separates the primitive from the secondary invariants. This set is however not optimal for our algorithm developed in section 4.3. based on integration by parts with the operator an invariant polynomial. In fact, in this casedP is a semi-invariant polynomial. A semiinvariant polynomial is an invariant polynomial up to a phase (in our case a minus sign) for some of the transformations of D 4 . Thus by means of d one can associate to any invariant polynomial not in the kernel of d, a semi-invariant polynomial. In our case, this suggests to introduce the semi-invariants: 14) It is easy to check that the Jacobian of the transformation (w, w ,w,ŵ) → (w, g, j, j ) is nonvanishing, which implies that the semi-invariants (w, g, j, j ) are algebraically independent. Now since dw = 4 , dj = 0 , (D. 15) this suggests to consider the new set of primitive invariants of degree respective 1, 2, 2, 4: (1 : a, 2 : v, 2 :w 2 , 4 : j 2 ), where we have further introduced v = 2λ + µ, which now are all butw in the kernel of d. Again the Jacobian of the transformation (w, w ,ŵ,w) → (a, v,w 2 , j 2 ) is found not to vanish, so that this new set does not satisfy any polynomial relation and thus constitutes an alternative set of primitive invariants, well-adapted to d. Thus any invariant can now be written as I(w, w ,w,ŵ) = f (a, v,w 2 , j 2 ) + g(a, v,w 2 , j 2 )c , (D. 16) with some polynomials f and g. In this new basis the syzygy is given by Since dc = −2g we can further refine the procedure. Introduce h = ac + g 2 which is in the kernel of d and using the relationships we can eliminate first c and then g so that the polynomial in (D.16) reduces to I(w, w ,w,ŵ) = f (a, v,w 2 , j 2 ) + g(a, v,w 2 , j 2 ) 4a 2 v + j 2 − a 4 − (aw − j) 2 16a = P (a,w, v, j) a .
(D. 19) Here P is a polynomial in the variables (a,w, v, j) which are algebraically independent since the Jacobian ∂(a, v, j,w) ∂(w, w ,w,ŵ) = −32(w − w )(w −ŵ) = 0 .   Since the polynomial f and g are uniquely defined, and since there are no relation between the variables (a,w, v, j) the polynomial P is also uniquely defined. Note that in this "basis" I has a non-polynomial expression because we divide by a. However this does not matter, since a will eventually be replaced by unity anyway in the further procedure of evaluating diagram B . The benefit is that now all polynomial invariants are expressed in terms of (a,w, v, j) which are all, exceptw, in the kernel of d. This is our optimal "basis". Note that (a, v) are invariant polynomials whereas (w, j) are semi-invariants polynomials. Note also that more precisely P depends onw 2 , j 2 andwj which are invariant under the whole action of the group D 4 .

E Summary of results on multiloop Euler-Heisenberg Lagrangians
In the following table we summarize the available information on multiloop Euler-Heisenberg Lagrangians, classified according to loop order, spin, dimension, and background field. The column "Computation" points to the articles where the Lagrangians were computed, "Discussion" to articles linking them to the exponentiation conjecture. For an exhaustive discussion of Euler-Heisenberg Lagrangians at the one-loop level see [3].
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.