QED as a many-body theory of worldlines: I. General formalism and infrared structure

We discuss a reformulation of QED in which matter and gauge fields are integrated out explicitly, resulting in a many-body Lorentz covariant theory of 0+1 dimensional worldlines describing super-pairs of spinning charges interacting through Lorentz forces. This provides a powerful, string inspired definition of amplitudes to all loop orders. In particular, one obtains a more general formulation of Wilson loops and lines, with exponentiated dynamical fields and spin precession contributions, and worldline contour averages exactly defined through first quantized path integrals. We discuss in detail the attractive features of this formalism for high order perturbative computations. We show that worldline S-matrix elements, to all loop orders in perturbation theory, can be constructed to be manifestly free of soft singularities, with infrared (IR) divergences captured and removed by endpoint photon exchanges at infinity that are equivalent to the soft coherent dressings of the Dyson S-matrix proposed by Faddeev and Kulish. We discuss these IR structures and make connections with soft theorems, the Abelian exponentiation of IR divergences and cusp anomalous dimensions. Follow-up papers will discuss the efficient computation of cusp anomalous dimensions and universal features of soft theorems.


Introduction
The infrared (IR) and ultraviolet (UV) divergences that appear in perturbative computations of amplitudes in gauge theories and gravity are essential to understanding the fundamental structure of quantum field theories. While the UV structure of amplitudes in the standard model and gravity, are the primary focus of attention due to their sensitivity to novel physics, IR divergences are also fascinating because of their importance in understanding the overall consistency of the respective quantum field theories and the possible existence of universal features and emergent dynamics that they may share.
The IR structure of amplitudes in QCD and gravity have a number of features that are not fully understood and are a subject of active research. In the Bjorken asymptotics of QCD, understanding the absorption of IR divergences into universal non-perturbative structures are key to active work on factorization theorems that are essential to the predictive power of amplitudes [1]. In the Regge limit of the theory, a "reggeization" of amplitudes is seen to occur, with the emergent Reggeon and Pomeron degrees of freedom capturing the leading IR dynamics of amplitudes [2]. A further striking feature is the phenomenon of gluon saturation, where an emergent semi-hard saturation scale characterizes the maximal close packing of soft gluons [3,4]. In gravity, as in QED, soft theorems govern the structure of perturbative scattering amplitudes [5]. However, unlike QED, the self-interactions of soft gravitons produced in trans-Planckian scattering can also lead to reggeization [6,7] and novel emergent structures, an example of which are a description of Black Holes as overoccupied self-bound gravitons [8], in close correspondence with gluon saturation [9]. Further, it has been argued that the existence of Black Holes is sufficient for the UV completion of gravity [10], directly relating the UV to the IR, and signalling the breakdown of Effective Field Theory [11,12].
In contrast, the IR structure of QED is much simpler and does not pose any practical problem for the computation of scattering amplitudes. Nevertheless, there are nagging conceptual issues, a clearer understanding of which may provide deeper insight into their more challenging counterparts in QCD and gravity. One formulation of the IR problem [13] states that an infinite number of real gauge bosons, with arbitrarily low energies, accompany the in and out states in any process. Hence it is not possible to precisely define, in the Fock sense, the interacting state of a system of charged particles. Moreover, virtual IR divergences in any scattering transition set the corresponding S-matrix elements to zero [14].
In Abelian theories, the spectrum of real and virtual IR quanta agrees with the classical expectation: it is independent of the spin dimension of matter fields and dominated by long wavelength bosons of small momenta k µ that are sensitive only to large spacetime regions x µ of the interaction, k µ x µ ∼ 1. The amplitude for the emission or absorption of IR quanta can therefore be computed without precise knowledge of the physics inside the region x µ [15]. Its momenta are peaked along the momenta of the charged incoming and outgoing particles at the boundaries of that region, and its contribution to the amplitude can be isolated as an overall soft factor, order-by-order in perturbation theory, from the amplitude without the IR bosons attached [5]. In non-Abelian theories, the self-interactions of gauge bosons makes isolating such soft factors, a much more involved problem in general.
The solution to the stated IR problem in QED is usually understood to come from the unitarity of scattering amplitudes. Bloch and Nordsieck [14] showed that transitions to final states containing an arbitrary number of soft photons are given by the exponentiation of the single photon emission cross-section, which diverges logarithmically as ∼ dω k /ω k in the IR. When virtual photons are also attached, and one sums over all possible configurations with real and virtual quanta above the IR cut-off, the IR divergences cancel 1 . Hence the total cross-section evolves depending only on the residual UV cut-off Λ to some power a(γ ij ). This function depends on the initial and final momenta of the charges through geometrical factors γ ij known as cusp anomalous dimensions.
The radiative cross-sections are thus renormalizable: for any other Λ , they scale by a factor (Λ /Λ) a , independently of the physics of the IR. For instance, if one wishes to represent the radiative cross-section in terms of the non-radiative one, one can choose as Λ the mass of the charges in the interaction [5]. Hence through renormalization, radiative cross-sections are finite, with the relevant scales set by the cross-section without real or virtual bosons. This cancellation of IR singularities was proved to all orders by Jauch and Rohrlich [20,21] and by Yennie, Frautschi and Suura [22]. The modern treatment demonstrating the classical features of soft factors in Abelian gauge theories is due to Weinberg [5]. In particular, he showed that this behavior is universal to gravity even though gravitons can themselves emit softer gravitons; this is because the effective coupling of soft gravitons is proportional to their energy, and therefore vanishes in the IR limit 2 .
The remarkable thing about this cancellation of IR singularities in the cross-section is that it does not naturally occur at the amplitude level. IR singularities perturbatively are still infinite order by order, or if summed to all orders, set all S-matrix elements to zero. So despite powerful proofs of consistency, the ubiquitous cancellation of IR divergences raises the question of why observables are always found IR finite in nature but this finiteness is not manifest at the amplitude level. In perturbative calculations, interactions are typically switched off at very late and early times relative to the timescales of the process under consideration. However because softer and softer virtual and real bosons can be attached to the process under consideration, no practical definition of the in and out states can be accomplished without taking into account these interactions at asymptotic times.
Building on the work of Dollard [23], Chung [24] and Kibble [25][26][27][28], Faddeev and Kulish (FK) [29] realized 3 that when the in and out states are suitably dressed to take into account asymptotic interactions, IR divergences can be systematically absorbed in these dressed states, order to order in perturbation theory. Specifically, in this FK picture, the wave functions of the incoming and outgoing charged particles are dressed with coherent states 1 Most generally, as argued by Kinoshita, Lee and Nauenberg, for theories with massless charged particles, this requires a summation over both initial and final degenerate energy states [16,17]. For recent discussions of this KLN theorem, see [18,19] . 2 This is not true for gluons emitting softer gluons, for which case, as Weinberg states, (the) elimination of such complicated interlocking infrared divergences would certainly be a Herculean task, and might even not be possible. This is the problem of reggeization and gluon saturation in QCD alluded to previously; while much progress has been made, the treatment of IR divergences is still by no means fully understood. 3 Followed by Zwanzinger in [30]. See also [31,32] for modern treatments of the problem.
corresponding to a cloud of infinite numbers of soft photons. The corresponding S-matrices of the dressed states are infrared finite. The handling of IR singularities in many practical circumstances should benefit from the definition of IR finite S-matrix elements. As we will discuss at length in this paper, and in follow-up work, both vacuum matrix elements of operators and scattering matrix elements can be expressed in an exact worldline formulation of QED as expectation values of so-called Wilson loops, and lines, respectively, to all orders in perturbation theory. The cusp anomalous dimensions we mentioned previously in the context of total cross-sections have a natural interpretation at the amplitude level in terms of the IR and UV properties of Wilson lines, which can be shown explicitly to be IR finite to all orders in perturbation theory. Thus scattering cross-sections can be understood in QED entirely in terms of these IR finite structures, with all non-trivial information contained in their UV behavior, as sketched previously.
Besides the intrinsic interest in such a manifestly IR finite formulation of amplitudes in QED, the Wilson line is also a fundamental object in Yang-Mills theories [33,34], where their UV and IR properties can be formulated in terms of evolution equations for the cusp anomalous dimensions [35,36]. The Wilson line defined as such can be employed to derive the DGLAP and BFKL evolution equations for parton distributions [37][38][39], the low energy behavior of form factors [1], Drell-Yan processes and Higgs production or the re-summation of large Sudakov logarithms [40,41]; for recent work on computing cusp anomalous dimensions to high orders in QED, QCD and N = 4 supersymmetric Yang-Mills theory, see [42][43][44][45][46]. There is likewise work in parallel on understanding the cusp anomalous dimensions of gravity amplitudes [7,47,48] and possible color-kinematics dualities between such infrared structures in QCD and gravity [49], which show great potential in simplifying computations of gravitational radiation amplitudes [50].
In addition to their key role in precision computations, the resurgence of interest in soft theorems and FK dressings in gauge theories arose independently motivated by developments in the context of gravity. It was shown around the same time as Weinberg's work by Bondi, Metzner, Van der Burg, and by Sachs [51,52] (BMS), that degenerate vacua in gravity separated by differing numbers of soft gravitons are connected by supertranslations corresponding to spontaneously broken symmetries of the infinite dimensional BMS group of asymptotic spacetime symmetries 4 . More recently, it was shown by Strominger and collaborators that Weinberg's soft graviton theorem could be understood as the Ward identity corresponding to BMS invariance of gravitational amplitudes [54,55]. They showed subsequently that the Weinberg formula is equivalent to that of gravitational displacement (or memory) experienced by two inertial detectors following the passage of a pulse of gravitational radiation [56]. This "triangle" of connections between asymptotic symmetries, soft theorems and memory has been conjectured by Strominger to be a universal feature of the infrared behavior of gauge theories and gravity [57].
Indeed, it was shown in [54] that scattering amplitudes in QED can be understood to satisfy an infinite dimensional symmetry group of large U (1) gauge transformations that approach angle dependent constants on the celestial sphere at null infinity. As with the BMS group in gravity, the spontaneous breaking of this symmetry by the QED vacuum results in an infinite number of degenerate vacua characterized by differing numbers of soft photons, which are Goldstone modes on the celestial sphere. The Ward identity [58][59][60] corresponding to the QED BMS-like symmetries is the soft photon theorem [5,15]. Completing the triangle is the QED memory effect [61,62]. This approach reimagines our conventional understanding of infrared divergences in the language of symmetries of the S-matrix of charged particles. Specifically, the usual S-matrix vanishes for the emission of ultrasoft photons precisely because it does not allow for vacuum transitions between the degenerate vacua which are allowed by the BMS-like symmetry. Likewise, the FK dressed states are coherent states of Fock vacua that are eigenstates of the BMS charges at null infinity. The FK cancellations of infrared divergences of the S-matrix by the soft photon dressings of the incoming and outgoing charges is simply reinterpreted as the statement that matrix elements of the S-matrix between the dressed eigenstates is finite [60,63]. This FK phase is the previously noted Goldstone mode of large gauge transformations satisfying a two-dimensional current algebra on the celestial sphere at null infinity [64]. A further intriguing prospect is the possibility of a clear dictionary between this language of asymptotic symmetries and that of Wilson lines and cusp anomalous dimensions [65].
A question that one can ask about this program of constructing and exploring infrared finite amplitudes is their validity to all orders in perturbation theory and, as a corollary, its relative efficacy compared to conventional approaches in the computation of cross-sections for physical processes. As indicated earlier, we will show here that this question can be cleanly formulated and addressed within the worldline formulation of QED, which is exact. While this framework for computing cross-sections in terms of manifestly IR safe structures is robust, its advantage remains to be established, though, as we will discuss, there are good reasons to believe that this may be the case.
The formulation of QED as a theory of first quantized 0+1-dimensional worldlines can be traced back to the original work of Feynman and Schwinger [66][67][68]. However a complete treatment of spin (and color) in terms of Grassmann worldlines was only developed in the 1970s [69][70][71][72][73][74][75]. The modern treatment of worldlines in gauge theories, following Polyakov's pioneering work [34], is due to Strassler [76], who also showed its equivalence to the powerful string amplitude formalism of Bern, Dixon, Dunbar and Kosower [77][78][79][80] for multi-leg loop amplitudes. A classic review of the subject is due to Schubert [81]; for a recent update, see [82].
In the worldline formulation of QED, matter and gauge fields can be integrated out explicitly, resulting in a many-body theory of point particle bosonic and Grassmann worldlines, with the multiplicity of the closed worldlines growing with each loop order and with the external particles described by a fixed number of open worldlines. The propertime variable controls the evolution of these first quantized worldlines; UV and IR scales are encoded respectively in the short and long time behavior of the Abelian classical Coulomb interaction, that emerges naturally exponentiated.
The UV singularities found in these worldline interactions can be mapped on to the UV divergences found order by order in perturbation theory and treated with standard renormalization methods. The striking feature of S-matrix elements, when expressed as worldlines, is that they can be constructed to be manifestly free of soft singularities, with IR divergences captured and removed by endpoint photon exchanges at infinity. As we will see, these asymptotic contributions exactly match the soft coherent dressings of the Dyson S-matrix proposed by Faddeev and Kulish [29] to define IR finite amplitudes.
It can be argued that the Faddeev-Kulish dressings only reshuffle IR divergences from the amplitudes to the states-see [32] for a nice discussion. However the same conceptual issue can be raised for any renormalization procedure eliminating singularities, be they IR or UV. Worldline amplitudes offer a covariant proper-time method to introduce IR and UV cut-offs in perturbation theory, with a transparent physical interpretation: the possibility of further interpretation in the language of symmetries makes it all the more attractive. In the worldline approach, both self-energy corrections to individual worldlines, and virtual exchanges between different worldlines, are treated on a common footing. As noted in [83], the improper treatment of these leaves the cancellation of IR divergences incomplete.
Worldline methods in QED have been employed previously to compute the Schwinger pair production rate in background fields [84], utilizing the method of worldline "instantons" (in propertime) to extend the computation to arbitrary coupling strengths when the background fields are weak [85][86][87]. Besides such semi-classical approaches, there have also been significant developments in worldline numerical Monte Carlo techniques [88][89][90]. Though these were previously applied primarily for studies of pair production in background fields, applying these techniques to compute high order contributions scattering amplitudes in the formalism we will develop is worthy of further exploration.
The formulation of QED (in particular its IR structure) as a many-body theory of worldlines may find an important application in quantum computing where such methods were previously discussed in the context of a systematic mapping the Wigner-Weyl-Moyal formalism [91,92] for the phase space distributions we discussed earlier to the full set of Clifford gates. Similar single particle digitization frameworks have since been developed for scattering in a φ 4 scalar field theory [93] and a hybrid framework for deeply inelastic scattering in Regge asymptotics [94]. Our many-body framework allows one to extend these ideas to QED; for earlier work in this direction, see [95]. Non-perturbative worldline approaches [96] are also under active study in the context of quantum computing. Relating these approaches to our discussion offers a possible path towards quantum computation of scattering amplitudes.
Another attractive feature of this formalism is its immediate generalization to Yang-Mills theories, allowing one to treat on common footing the exponentiation of spinor, helicity as well as color degrees of freedom through the use of Grassmann variables [71,97,98]. In particular, the fact that colored Wilson lines can also be expressed in terms of the Grassmannian counterparts [36,[99][100][101][102] of the spin degrees of degrees of freedom does not appear to be widely appreciated. The Grassmann framework also admits an elegant semi-classical interpretation 5 in terms of spinning, colored particles in background fields [73,74,110,111]; a canonical coordinate transformation of these variables to Grassmann bilinears representing semi-classical spin vector and color charges allows one to construct, from first principles, extended phase space distributions including these so-called 6 Darboux variables [92].
While a detailed discussion is beyond the scope of our discussion here, we note that the worldline representation of colored Wilson lines in terms of a path integral over a local onedimensional Grassmann action presents potential advantages in wide range of practical QCD computations ranging from a chiral kinetic theory of the Quark-Gluon Plasma [113,114], to deeply inelastic scattering in the Regge limit [115], to the role of the chiral anomaly in the proton's spin [116,117]. For computations of IR soft factors for the problems we mentioned earlier, unlike the QED case, one does not have an exact exponentiation in QCD [33,118,119], with IR factorization requiring a systematic eikonal and next-to-eikonal exponentiation of socalled web connected diagrams [120][121][122][123]. The web structure of non-Abelian exponentiation is reproduced in a generating functional formalism which allows for matrix exponentiation of the products of Wilson lines [124]. The representation of Wilson lines in terms of a local Grassmann action suggests that our framework can aid in simplifying such computations.
We envisage this paper as the first of several that will explore in some depth the infrared structure of QED as a many-body theory of worldlines. We will develop the formalism in this paper and outline some of the key features that are relevant for the IR issues we have noted. Paper II will discuss at length the computation of real soft infrared divergences and present a detailed and efficient computation of the two loop cusp anomalous dimension. Future papers will discuss the extensions of this formalism to higher orders and to QCD. An intriguing offshoot of this work that we noted, and will pursue separately, is in the application of non-perturbative worldline methods to quantum computing.
The paper is structured as follows. In Sections 2 and 3, we will present a formal derivation of worldline amplitudes in QED to all orders in perturbation theory. In Section 2, we show that the QED path integral, involving closed vacuum loops to all orders, can equivalently be expressed in terms of first-quantized generalized Wilson loops of "super-pairs" of pointlike boson and Grassmann variables. This framework is then extended in Section 3 to consider general scattering processes, which can be formulated exactly, to all orders, as first-quantized point-like boson and Grassmann generalized Wilson lines and loops. In Section 4, we will take these results to familiar ground by reinterpreting them in the language of conventional Feynman diagrams. We will first show in Section 4.1 how the worldine structures obtained in the previous sections lead to universal "master" expressions for generalized polarization tensors that allow one to construct compactly QED amplitudes order-by-order in perturbation theory. This allows, in principle, efficient computations of QED processes at higher loops, where the factorial growth in the number of diagrams with each loop order in conventional perturbation theory significantly hinders the quest for precision. In Section 4.2, our results are expressed equivalently in terms of a generating functional for arbitrary background fields. pendix D) and Wong's equations in QCD [104]. In gravity, their worldline counterparts [105] are the Papapetrou-Mathisson-Dixon equations [106][107][108]; for a recent discussion and development of the latter in the language of effective field theories, we refer the reader to [109]. 6 Ref. [92] established the equivalence of the semi-classical limit of the worldline formalism to the Darboux phase space approach of Alekseev, Faddeev and Shatashvili [112].
As we will discuss, this approach is useful in strong field settings and potentially, in the use of worldline Monte Carlo methods in high order computations. Section 5 examines the IR structure of the QED Dyson S-matrix in the worldline formalism. We demonstrate in Section 5.1 a formal proof of the IR finiteness of the Faddeev-Kulish S-matrix to all orders in perturbation theory. The proof follows fundamentally from the simple asymptotic classical structure of the first-quantized worldline currents. Our result is illustrated for the case of Möller scattering in Section 5.2. In particular, we show how the asymptotic real and virtual contributions of the FK S-matrix arise in this approach, and can be expressed naturally in the first-quantized language of cusps in worldline currents. In Section 5.3, we connect our discussion of IR safety in QED to recent discussions and alternative approaches to this problem. In Section 6, we recover tbe well-known result for the one loop QED cusp anomalous dimension, and discuss its computation to higher orders. As noted earlier, a further detailed discussion will be the subject of Paper II. We briefly summarize our conclusions in Section 7, and provide an outlook on extensions and applications of our work.
The paper contains five appendices. In Appendix A, we state the conventions that are used in this work. Appendix B provides a detailed derivation of the all-order dressed spin-1/2 propagator, which clarified subtle features of a previous computation by Fradkin and Gitman [125]. The compact worldline result for n th rank polarization tensors is derived in Appendix C. As a limiting case of these general results, as previously shown by Strassler [76], one recovers the Bern-Kosower expression for the one-loop polarization tensor with arbitrary numbers of external photon legs [78]. Appendix D discusses the classical worldline equations of motion which, for homogeneous backgrounds, reduce to the well-known Bargmann-Michel-Telegdi equations [103] for spinning charges in background electromagnetic fields. The semiclassical worldline instanton method for arbitrary background fields is also outlined in this appendix. Finally, Appendix E provides details of key aspects of the derivation of the IR dressings in Section 5.2.

The QED vacuum as a many-body theory of closed worldlines
Following pioneering work on the worldline description of quantum field theory [76], there have been a number of significant advances in this framework, which are nicely summarized in [81,82]. However the fact that QED can be expressed exactly as a many-body theory of open and closed first-quantized worldlines is not widely appreciated, though, as noted, this possibility was already considered in Feynman and Schwinger's classic papers [66][67][68]. In this section, we will systematically develop the many-body worldline framework for the QED vacuum-vacuum amplitude. This will be extended to scattering processes in the next section. Our approach is pedagogical, requiring for the most part only basic familiarity with quantum field theory.
Consider the vacuum-vacuum matrix element Z = 0|0 in QED, encoding the dynamics of an infinite sea of virtual particles, Here g = ±e, ζ is the gauge fixing parameter, / D = D µ γ µ with D µ = ∂ µ − igA µ (x), and [D µ , D ν ] = −igF µν (x). We will work in Euclidean spacetime, wherein this non-perturbative amplitude is well defined. Perturbative amplitudes can be Wick rotated back to Minkowski spacetime using the conventions given in Appendix A.
To express Z entirely in terms of worldlines, we must integrate out the A µ (x) and Ψ(x) fields. Integrating out the matter fields first, we get Before integrating out the gauge fields A µ (x), we will first write the fermion determinant in worldline form. To do so, we use the identity 7 , with The fermion determinant term in the path integral then becomes As is customary, to compute the trace of this operator, we introduce the Schwinger propertime parameter ε 0 and noting that (valid for Re z > 0), we integrate both sides between z = 1 and z = a, to obtain The square of the Hermitian operator / D + m is positive, and the previous heat-kernel regularization method can be applied. The functional trace in Eq. (5) becomes 7 In the first step, we introduced γ 2 5 = 1 and then used {γ 5 , γ µ } = 0.
allowing one to reduce the problem of computing the fermion determinant to one of finding the matrix elements of a first-quantized fermion traversing a loop in the background A µ (x), governed by the Hamiltonian operator where we identifiedp µ = i∂ µ . The first term in the r.h.s of Eq. (8) acts as the required UV cut-off of the effective action, renormalizing the genuine UV divergence of the zero-point energy of the vacuum. The vacuum-vacuum amplitude now reads, To perform the trace in the expression above, the matrix elements of e −Ĥ can be expressed as point particle boson and fermion path integrals [67,76,81,98,[125][126][127][128][129]. The eigenfunctions of the fermion operatorσ µν = i[γ µ , γ ν ]/2 can be treated on equal footing with the particle 4-positionx and 4-momentump bosonic operators [71,97]: spanning the space of intermediate states of worldlines propagating in propertime τ in the gauge field background A µ (x). Our starting point is the result in [76], which can be expressed in general as where we gauge-fixed the einbein ε 0 so that τ ∈ [0, 1]. The 0+1-dimensional bosonic commuting worldline x µ (τ ) encodes the particle 4-position and its fermion anticommuting worldline counterpart ψ µ (τ ) encodes spin precession. The latter inherits the anti-commutation properties of the Dirac matrices and their action on a single particle state. The two currents coupling to A µ (x) are respectively, i) of the standard form gẋ µ (τ ) in the Wilson loop/line formulation of gauge theories (denoting the Lorentz force of a scalar charge in a background gauge field), and ii) the coupling of the local spin tensor σ µν (τ ) = i[ψ µ (τ ), ψ ν (τ )]/2 to the local magnetic and (boosted) electric components in F µν (x). If we impose periodic boundary conditions (P) x µ (1) = x µ (0) and anti-periodic boundary conditions (AP) ψ µ (1) = −ψ µ (0) respectively on the boson and fermion worldline trajectories in Eq. (12), they describe the sum over all possible configurations in propertime of worldlines representing the dynamics of virtual spin-1/2 charges in an arbitrary background A µ (x). Inserting Eq. (12) in Eq. (10), and expanding in loops, one obtains To perform the remaining integration in A µ (x), we express the free Maxwell action as where we integrated by parts to find the inverse of the gauge-unfixed Euclidean photon propagator, with the matrix elements Taking its Fourier transform, and inverting the result, the propagator in d dimensions is given by Finally, performing the momentum integrals for arbitrary values of the gauge parameter, its coordinate space counterpart is We shall also express the worldline currents as 4-densities to integrate out the gauge field. Eq. (13) can be cast as where we express the sum over the currents of the virtual charges as with the charged scalar and fermion current densities, respectively, Thus with the help of Eqs. (14) and (19), the vacuum-vacuum amplitude in Eq. (13) can be reexpressed as We can now integrate over A µ (x) and write this expression (normalizing with respect to the pure gauge vacuum sea of disconnected photon loops), The -th contribution (W ( ) ) in this loop expansion is the amplitude for finding the QED vacuum in a configuration with virtual charged particles. It can be represented as the sum over all possible worldline paths x i µ (τ ) and ψ i µ (τ ) with periodic and anti-periodic boundary conditions, respectively, of virtual worldlines describing loops in propertime ε i 0 , summed over all possible propertimes, and takes the form The dynamics of these virtual fermions are encoded in the exponentiation of the photon exchanges between the -loop currents in Eq. (19), functionals of the set of super-pairs of virtual worldlines {x µ (τ ), ψ µ (τ )}. Using the r.h.s of Eq. (19) one may express the exponential as that of the sum of contributions due to the exchanges between pairs of worldline currents i and j: Represented thus, one observes that the worldline dynamics includes the exponentiation of both the self-energy subgraphs i = j and the non-diagonal exchanges i = j on equal footing, consisting of the Lorentz forces between two spinning charges i and j. This will be important to note for future reference when we discuss in Section 5 the FK dressings of worldlines. Even though we "got rid" of the gauge fields, we can reintroduce them in a purely first quantized interpretation of the th -loop contribution to the QED vacuum-vacuum transition amplitude as a many-body theory of worldlines. Specifically, with Eq. (20), the boson-boson coupling between the coordinate space currents can be rewritten as where (for ζ = 1 in d = 4) using Eqs. (16) and (20), Here A x j µ (x) the field at point x created by the coordinate space current gẋ j µ (τ ) of the j th charge moving along x j µ (τ ). Likewise, the spin precession of the i th worldline couples independently to the magnetic and (boosted) electric components of the field created by the scalar current gẋ j µ (τ ): where we used Eq. (20), integrated by parts, and antisymmetrized to define the local spin tensor of the i th virtual charge σ i µν (τ ) = i[ψ i µ (τ ), ψ i ν (τ )]/2. The field strength tensor created by the vector current x j (τ ) refers to with A x j µ (x) given in Eq. (26). The couplings of the spin of worldline j with the field emerging from the spatial current of i can be rewritten using Eq. (20) as where the field created by the spin precession of j is We used here Eq. (16) with ζ = 1 and d = 4. Lastly, the pure spin-spin interaction of the worldlines, using Eq. (20) again, can be written as where the field strength tensor created by the spin precession of worldline j at the point x is with A ψ j µ (x) given in Eq. (30). Putting together the contributions from Eqs. (25), (27), (29) and (31), the interaction functional in Eq. (23) can now be rearranged as a more general Wilson loop describing the dynamics of particle i: It is important to note that in this generalized Wilson loop 8 the background gauge fields in Eq. (26) and (30) are replaced by dynamical fields. They include the self-interaction terms i = j, the spin precession is exponentiated on the same footing as the dynamical fields, and worldline contour averages exactly defined through path integrals. These features are not accessible in the conventional Wilson loop formulation of gauge theories. Using Eq. (26) and (30), with Feynman gauge ζ = 1 and d = 4, the generalized Wilson loop takes the elegant form where the structure of the Wilson loop is that of a normalized expectation value O ( ) of a many-body functional O ( ) of the super-pairs of worldlines {x µ (τ ), ψ µ (τ )}, with the expectation value defined as This allows us to write the QED partition function (with the shorthands x i µ = x i µ (τ i ) and σ i µν = σ i µν (τ i )) in the following compact and suggestive form: Eq. (36) is the central result 9 of this section. It expresses the QED vacuum as a many- 9 We will not discuss here UV renormalization in the sense of mass and charge renormalization. Worldline perturbation theory, as we will discuss in Section 4, consists of explicit computations of generalized polarization tensors. Employing dimensional regularization in this framework, one can identify the corresponding UV divergences as the usual ones extracted from Feynman diagram computations [81]. Renormalization of generalized Wilson lines [34] will be discussed in Section 6. body theory of pointlike 0+1 dimensional particles represented by "super-pairs" of closed boson and fermion worldlines 10 (x i µ (τ ), ψ i µ (τ )) with typical virtualities 1/ε i 0 . Eq. (36) is exact and generalizes the notion of a Wilson loop to spinning, precessing particles in fully dynamical backgrounds, with the field of the ith worldline sourced by self-interactions as well as those generated by the other particles at any given order in the loop expansion. An illustrative representation of one of these amplitudes is depicted in Fig. 1.
While it is perilous to claim anything novel in such a highly developed field, we have not, to the best of our knowledge, seen the explicit and compact many-body worldline representation given in Eq. (36) elsewhere in the literature. In the next section, we will extend this framework to open worldlines, thereby enabling one to compute many-body scattering processes.

S-matrix in the many-body worldline formalism
We will consider here the scattering amplitudes for transitions between in and out charged states via the emission of virtual photons 11 . We denote n i and n o (n i andn o ) to be the number of incoming and outgoing charged particles (antiparticles) and N i = n i +n i and N o = n o +n o the total number of in and out charges. We also identify a pair of external particles with indices n and m, with n, m = 1, . . . , r, with r = (N o + N i )/2 being the number of real worldlines in the scattering process. A pair of virtual particles is identified with indices i and j with i, j = 1, . . . , , where represents the number of virtual worldlines.
The Dyson S-matrix element for the transition amplitude i → f is defined as and the limits t f → +∞ and t i → −∞ are assumed. Here |p, s and |p,s are single particle and antiparticle states of momentum p and s. The vacuum-vacuum amplitude calculated in Section 2 corresponds to the particular case r = 0 in Eq. (37) with |in = |0 and out| = 0|. Our strategy here is to continue working in Euclidean spacetime to construct the nonperturbative building blocks, which can then be Wick rotated to Minkowski spacetime, as needed, to compute scattering amplitudes to a given order in perturbation theory. Our starting point is the Euclidean QED path integral whereη(x) and η(x) are anti-commuting external sources. Proceeding along the same lines as in the previous section, and first integrating out the Dirac fields, we get The dressed Euclidean fermion Green function D A F (x, y) satisfies We shall now construct dressed N -point Green functions and normalize the amplitudes for these with respect to the infinite sea of virtual loops in Z = 0|0 ≡ Z[0, 0]. The 2-point function is given by where the r.h.s, defined as corresponds to the propagator for a single real fermion to go from position x 1 i to x 1 f while coupled to a fully dynamical gauge field. It contains arbitrary photon insertions and fermion loops, the latter arising from the fermion determinant.
Higher order N -point functions describe the interacting problem of r = (N i + N o )/2 real spinning charges. (Odd N -point functions vanish.) The 4-point function gives The relative sign between the two terms follows from functional differentiation, yielding the odd symmetry factor of the final state under the interchange of single final states. Likewise, higher N -point functions and their matrix elements can be built as functional averages of products of 2-point functions given by D A F (x f , x i ). Thus to compute N -point amplitudes in the worldline formalism, it is sufficient to express D A F (x f , x i ) in terms of worldlines, and subsequently perform the functional average over the gauge fields.
Our starting point here is a result first obtained by Fradkin and Gitman [125], which we will rederive in Appendix B to clarify some aspects of the derivation 12 . This result (Eq. (226) of the Appendix) gives the amplitude for the propagation of a spinning charge from The various terms and fields in the above expression demand detailed explanation, which is given below. (See also Appendix B for further details.) The propagation of the real fermion is described by the super-pair (x µ (τ ), ψ µ (τ )) of commuting and anti-commuting worldlines, respectively, with the dimensions of the fermion degrees of freedom specified by five Grassmannian fields labeled by λ = 1, 2, 3, 4, 5, with the fifth necessary to account for the helicity-momentum constraint 13 . Importantly, one requires open boundary conditions of their path integrals, which satisfy where θ λ are five anti-commuting auxiliary sources, which are eventually set to zero. We also defineγ µ = −iγ 5 γ µ andγ 5 = γ 5 , with γ 5 = γ 1 γ 2 γ 3 γ 4 ; they obey the identity {γ λ ,γ ξ } = 2δ λξ . These two objects anti-commute ({θ λ ,γ λ } = 0), which allows us to generate the correct tensor structure of the Green function in the fermion degrees of freedom. The normalization of the odd path integral in five dimensions is denoted by N 5 . Further in Eq. (44), A µ -independent terms are collected in the (super-gauge unfixed) free worldline action of a real spinning charge We leave manifest the time reparametrization invariance of the worldlines as a gauge symmetry, keeping the einbeins unfixed: besides the super-pair of worldlines x µ (τ ), ψ µ (τ ) , we introduced a new super-pair of commuting and anti-commuting einbeins denoted ε(τ ) and χ(τ ), respectively, with trivial dynamics. Their conjugate partners are denoted π(τ ) and ν(τ ). The first two terms in the curly brackets of Eq. (46) fix the gauge in the reparametrization invariance of the worldlines. The next two terms in the curly brackets of Eq. (46) enforce the energy-momentum constraint on the worldline, while the terms in the parenthesis in Eq. (47) impose the helicity-momentum constraint [113,133].
The zero-modes of the new boson and fermion einbeins are denoted as and correspond to conventional commuting and anti-commuting Schwinger parameters, respectively, in Feynman diagram language. The terms π(τ )ε(τ ) and ν(τ )χ(τ ) in Eq. (47) act as super-gauge fixing terms.
With all the terms specified as above, Equation (44) is exact. It sums all possible worldline paths connecting the open boundary endpoints specified in Eq. (45) to construct the dressed quantum amplitude for a spin-1/2 particle to go from x i to x f in the presence of an Abelian background gauge field, to all orders in the coupling. It also contains the semi-classical limit of the action of a charged spinning particle in an Abelian background gauge field, proposed previously in [69,74,135].
For a single charge r = 1, plugging Eq. (44) into Eq. (42) gives the dressed fermion propagator with the real particle charged current density defined (with the einbein ε(τ ) unfixed) to be Further expressing the fermion loop contributions as functional traces in worldline form, and expanding in loops, precisely as in the previous section, we then get The charged current densities of the virtual fermions are defined identically and the corresponding free worldline action given by Substituting Eq. (51) into Eq. (49), we get Note that the gauge functional average differs from the one in the vacuum-vacuum amplitude in Eq. (21) only due to the presence of the charged current of the real particle. Since currents and actions are additive for each virtual and real particle present, we defined Upon integrating over the gauge field, as previously, we finally obtain the dressed propagator for the propagation of a real spinning charge from x 1 i to x 1 f , which can be expressed as and where, in analogy with the vacuum-vacuum amplitude, we can define a more general Wilson line containing virtual fermion loops as The structure of Wilson lines and loops are that of the normalized expectation value O (r, ) of many-body functionals O (r, ) of the set of r and super-pairs of fermionic and bosonic worldlines {x µ (τ ), ψ µ (τ )}, with boundary conditions in Eq. (45) for the r real charges and P and AP for the virtual charges, Eq. (56) then expresses the functional average over fields in the l.h.s., as a 0+1 dimensional field theory in the r.h.s., with normalized expectation values W (1, ) (x f , x i , θ) given by the path integrals in Eq. (58). Note that in the above expression the sub-indices i and n denote different and independent (virtual and real, respectively) worldlines.
We can now write down the precise form of the 2r point function, encoding the amplitude of propagation of a system of r real particles 14 where No...1 is the totally anti-symmetric symbol and the sum runs over all permutations of the final points of the dressed Euclidean Green functions, encoding the parity of the out many-body wavefunction under permutations of final single particle states. Following the same steps as for the r = 1 case we get where This functional contains the exponentiation of all tree-level photon exchanges in the selfinteractions of the charged fermions as well as of the virtual photon exchanges between them. It further includes the fermion loops arising from the couplings of the real particles with the polarized virtual fermions of the sea and photon exchanges between the disconnected loops of the sea. These last are removed order to order in perturbation theory by (Z/Z MW ) −1 in Eq. (60), calculated in Section 2. With the input from the above expression, Eq. (60) is exact just as was the case for Eq. (44). It provides the 2 r-point correlator of r spin-1/2 fields between x n i to x n f to all orders in perturbation theory. It can be interpreted as a sum over worldline contours of the interactions of r pointlike spinning charges in the field created by themselves and those induced by virtual particles at each loop order. Note that the tree-level = 0 case of a single field r = 1 in Eq. (56) should be understood as the generalization (with full spin corrections, back-reaction, helicity-momentum and energy-momentum constraints, and well-defined worldline expectation values) of a Wilson loop operator between x 1 i to x 1 f . For > 1, we note that the the exponentiation of genuine fermion loops, as shown, cannot be obtained in the conventional Wilson loop picture.
We will now provide the explicit form of the Euclidean interaction functional of QED in general d dimensions, arbitrary gauge parameter ζ and general einbein ε(τ ). Using Eq. (55), we can reexpress the argument of the exponential in Eq. (61) as where S ab denotes the contribution to the action of the pair of worldlines a and b in the many-body system comprising the complete set of interactions in the amplitude. We can write it as where the labels (B) and (F) denote the boson and fermion currents, respectively, of a and b attached to the photon. Dimensionally regularizing Eq. (17), absorbing the mass dimensions of the coupling in µ and further using Eq. (20), we obtain the following individual contributions: which is the classical interaction appearing in a Wilson loop between two scalar charges. We use as a shorthanded x a µ ≡ x a µ (τ a ), ψ a µ ≡ ψ a µ (τ a ) and a ≡ a (τ a ). We observe that the non-diagonal gauge-dependent terms in the second line of Eq. (64) can be cast as a total double derivative with respect to particle a and b proper times τ a and τ b . Thus in the worldline formalism, gauge transformations correspond to the addition of boundary terms to the action [81], but only in the scalar interactions, as we will see next.
The coupling of the scalar current of a charge a, in the magnetic and boosted electric field created by the spin precession of b, generates and its mirror image (setting a ↔ b), coming from the coupling of the spin precession of a in the field created by the scalar current of b, .
The final term is the pure fermion-fermion contribution corresponding to the interaction induced by the precession of the spin of one worldline on that of the other: .
Note that due to the global supersymmetry of the QED worldline Lagrangian (whose explicit form can be found in Eq. (227) of Appendix B), the action between any two charges S ab in QED can be generated through the application of the 0+1-dimensional N = 1 SUSY algebra to the scalar QED term S BB ab [136,137]. For d = 4, after fixing Feynman gauge ζ = 1 and the einbein, these individual contributions can be put together in a compact expression as previously shown in Eq. (34) for the vacuum-vacuum amplitude. This can be seen by substituting Eqs. (64)-(67) into Eq. (62), and thence into Eq. (61).
Having obtained explicit expressions for all the elements necessary to compute W (r, ) in Eq. (61), we will now relate this to to the Dyson S-matrix element in Eq. (37). We start by evaluating the r = 1 case, describing the creation of a free particle at past infinity and its subsequent evolution via the QED Hamiltonian up to plus infinity. For a positive energy charge propagating from past infinity, one can represent the fermion dressed self-energy in terms of a Dyson S-matrix element, given by 15 where α, β and γ are spin indices, and we introduced the following shorthands for positive energy plane-wave functions: The QED generating functional Z[η, η] in physical time is given in Eq. (193) of Appendix B. The S-matrix can be reexpressed as which involves replacing the dressed Euclidean Green function in Eq. (60) with its Minkowski (228) of Appendix B multiplied by a factor of i. The integration of the dressed Green function over the gauge field configurations produces analogous to the one in Eq. (22) for the vacuum-vacuum amplitude, with the -loop contribution given by . 15 The external lines are created by imposing the limits x 0 f → ∞ and x 0 i → −∞ in the dressed N-point Green functions, which is equivalent to the standard LSZ reduction [83]. The consequences of this truncation will be very relevant for the discussion in Section 5 on the IR behavior of the S-matrix.
The first term = 0 contains the tree-level self-energy graphs attached to a single charge to all orders in perturbation theory. The = 0 terms encode self-energy contributions, also to all orders in perturbation theory, with a fixed number of virtual fermions. As discussed before, the disconnected vacuum polarization graphs are removed by the vacuum-vacuum terms in Z/Z MW .
Similarly, for an anti-particle created at past infinity, introducing the shorthands for negative energy plane-wave solutions the Dyson S-matrix element reads which can be thought of as the dressed Green function of a positive energy particle propagating backwards in time fromx f tox i , and multiplied by a factor of −1. Using Eq. (71), it can be rewritten as the loop expansion Turning now to the many-body scattering problem with r real particles, the calculation of the S-matrix elements involves just the reduction of the 2r-point Green functions to products of the (r = 1) 2-point functions above, as emphasized previously. We begin by computing explicitly the worldline Dyson S-matrix elements for the r = 2 case, depicted schematically in Fig. 2. In particular, Möller scattering (e − e − → e − e − ) between two fermions is given by the expression Proceeding along the same lines and using Eq. (71), this gives two contributions, related to the one possible permutation of the single particle states in a final state with two fermions,

S
(2) × lim × lim Along the same lines, Bhabha scattering case (e − e + → e − e + ) in the worldline formalism is given by the Dyson S-matrix element, , and contains two contributions The first contribution yields the amplitude for the scattering of the particle and anti-particle pair in the in state to the out state, to all orders in perturbation theory: The second contribution, the annihilation of the real pair in the in state to produce a different real pair in the out state, is given by These expressions can be generalized to scattering processes involving an arbitrary number of real charges r. For example, an in state with initial n i positive energy particles and an out state with n f = n i = r positive energy particles, has the Dyson S-matrix element, where the -loop contribution is given by Eq. (84) is the central result of this section. It provides a closed form expression for the matrix element for any α → β transition between r real fermions to all loop orders. It can be reinterpreted as the sum over first-quantized many-body worldline configurations of real and virtual pointlike spinning charges whose numbers grow with each order in the loop expansion. We observe that it exponentiates the hard and soft interactions of a many-body system, including the spin degrees of freedom, and defines the expectation value corresponding to averages over worldline trajectories that are seen explicitly to obey energy-momentum and helicity-momentum constraints.
In the next section, we will flesh out the perturbative framework obtained that follows from the formal expressions for the amplitudes derived in this section and in Section 2. We will specifically point to its advantages relative to the conventional Feynman diagram framework for computing amplitudes.

QED perturbation theory without Feynman diagrams
The central results of Sections 2 and 3, Eqs. (22) and (84), respectively, are valid to all orders in perturbation theory, providing a first-quantized many-body description, summarized in the latter equation, of an arbitrary number of photon exchanges between virtual and r real worldlines. To provide a familiar context in which to interpret these findings, we will connect our results to the conventional language of Feynman diagrams, and thereby examine in detail the perturbative expansion in the coupling of the worldline amplitudes found in Sections 2 and 3. We note that by considering general multi-loop and multi-open worldlines we are here going beyond the well-known one-loop approach to QFT without Feynman diagrams pioneered by Bern and Kosower, [78] and further elucidated by Strassler [76], albeit discussed here in the limited context of QED.

Multi-loop diagrams in the QED vacuum
We begin by recalling Eq. (23): with where the photon propagator is the expression in Eq. (16) and the Fourier transforms of the virtual fermion currents in Eq. (19) arẽ Because the required normalized expectation value in Eq. (86) is separable, loop-by-loop, to the product of one-loop normalized expectation values, one can alternately rewrite this expression as whose building block is the interaction functional representing single photon exchange between a pair of virtual fermion loops i and j, and Z ( ) (n) represents the vacuum sub-graphs with n photon lines connecting in all possible ways fermion loops; in other words, Z ( ) (n) is built from -loop n-photon polarization tensors. Thus the perturbation expansion in the coupling of the vacuum-vacuum amplitude can be expressed as (1) + · · · multi-photon, 1-fermion loop (1) + · · · multi-photon, 2-fermion loops (1) + · · · multi-photon, 3-fermion loops Note that here we have factored out the zero-point energy of the vacuum from all diagrams 16 , To compute in full generality terms in this expansion at arbitrary order, we first observe that 16 We omit in what follows the phase common to all diagrams.
where n 11 is the number of photon lines attached at both ends to the i = 1 loop fermion.
Repeating the procedure recursively for the remaining term in brackets in the r.h.s., this gives where n ij (n ji ) are the number of virtual photon lines flowing from the fermion loop j to i (i to j). Inserting this back into Eq. (88), the n! out in front cancels, and one gets Using now Eq. (89), collecting the currents of each virtual fermion i, and noticing the normalized expectation values are independent, the terms in the coupling expansion of the vacuum-vacuum amplitude of QED are given by Through these formal manipulations, we now see that the general -loop, n-photon, amplitudes are reduced to computing a product of one-loop, n th i = j=1 (n ij + n ji ) rank, vacuum polarization tensors.
In the worldline framework, these reduce to evaluating the normalized expectation value of a product of n i currents, specifically, the terms in brackets in the second line of Eq. (95). Here, each k of the n ij photon lines of momentum q ij k are attached at both ends to the i and j virtual fermions. The signs in the arguments of the currents −q ij k and +q ji k indicate, respectively, the absorption and emission of the momentum of the photon by the fermion line. Because there are n ij ! ways of attaching the two ends of n ij virtual photon lines and 2 momentum flow directions for each of the photons, the sub-graph is normalized by the factors n ij ! and 2 n ij . We absorbed the −1 parity factors in each vacuum polarization tensor, and the 1/ ! accounts for the ! identical permutations. The detailed computation of the worldline path integrals in these general n th rank tensors is spelled out in Appendix C.
To understand how the perturbative expansion works, let us first consider the fermion one-loop graphs ( = 1) depicted in Fig. 3. Beginning with the simplest graph, setting n = 0 in Eq. (95), the noninteracting single free fermion loop diagram can be expressed as where the infinite volume phase space factor V = (2π) 4 δ (4) (0) is a consequence of translational invariance, and where we used in the r.h.s. the result Eq. (237) derived in Appendix C. The remaining integral in ε 0 is UV divergent. However, one must pair this amplitude of finding a free one-loop in the vacuum with the corresponding counterterm in the zero-point energy of the vacuum, whose series will produce The UV divergence 1/ε 0 is then removed, leaving the mass-shell branch point In the absence of real photons, there are no diagrams involving tadpole subgraphs in the one-loop contributions = 1. They will appear however in multi-loop diagrams with > 1.
Since Π µ (k) = 0, as seen from the derivation following Eq. (246) of Appendix C, those containing tadpoles vanish as required.
The first correction, in Z (1) , is the amplitude for finding the vacuum contribution from a single fermion loop with the internal exchange of a photon. Setting n = 1 in the expression in Eq. (95), we get where we used Eq. (258) in Appendix C.
Higher order contributions are systematically obtained by substituting Eq. (259) from Appendix C in Eq. (95). For instance, the amplitude for finding the vacuum in a configuration with a single fermion loop exchanging two photons, where Π µνρσ (q 1 , q 2 , q 3 , q 4 ) is the light-by-light scattering graph that can be obtained from precisely the same master equation (Eq. (259) in Appendix C) as the vacuum polarization tensor. The one-loop diagrams Eqs. (96), (99) and (100) It corresponds to graph (a) in Fig. 4. More generally, all the free loop amplitudes satisfy Z ( ) (0) / !; using Eq. (97), they can be collected, resummed and UV regularized as From Eq. (95), for the two-loop graphs (graphs (b), (c) and (d)) with one photon exchange, we get Z The first term corresponds to graphs (b) and (c) in Fig. 4, and the second term to graph (d). Since iΠ µ (q) = 0, this last diagram does not contribute, and we can express the result as Z (1) = Z where Z (1) is defined in Eq. (99). The vacuum-vacuum amplitude for two fermion loops with two photon internal lines, illustrated in Fig. 4 graphs (e)-(h), can be read from Eq. (95) to be, (dropping the diagrams with tadpole subgraphs (e) and (f ) in Fig. 4) The first contribution corresponds to graph (g) in Fig. 4, the second term corresponds to the square of the second diagram in Fig. 3, and the final term is depicted in graph (h) of Fig. 4. The above exercise is very suggestive of the efficacy of the worldline framework for high order perturbative computations. Their advantage can be summarized as follows. Firstly, from the examples we considered, we can deduce that all multi-loop diagrams, to any level of complexity, can be universally generated from the master equation -Eq. (95). The explicit, yet compact, elements for their computation given by Eq. (259) in Appendix C. In general d dimensions, the results of applying this procedure to the calculation of amplitudes to fixed order are worldline diagrams that can be equivalently represented in terms of Schwinger/Feynman parameters, with their genuine UV loop divergences then regularized as in conventional perturbation theory.
Secondly, while each Feynman graph in the latter corresponds to one particular permutation of the photon insertions in Z (2) , shown in Fig. 3, is given by the integral in Eq. (259) of Appendix C. In the usual perturbative calculation, it requires the evaluation of two Feynman graphs, corresponding to the dressing of the fermion loop with the two photons crossed or uncrossed, respectively. At large n in perturbation theory, the situation becomes more complex, with the number of diagrams exploding factorially.
A fundamental question then arises as to whether Eq. (259) grows factorially with the increasing number of diagrams in the large n photon limit, or if cancellations do in fact occur between diagrams and reduce the scaling of Eq. (259) at large n from the naive estimate of the number of diagrams involved. For instance, Cvitanovic and Kinoshita [138] showed that in calculations of the higher order corrections to the electron anomalous magnetic moment [138][139][140] large cancellations occur at higher orders within certain gauge invariant sets of diagrams. Given that the number of diagrams inside these sets is very large (∼ hundreds to sixth-order in perturbation theory for g − 2) this led Cvitanovic to conjecture [138] that g − 2 may not satisfy an asymptotic expansion (at least in the quenched contributions with no virtual fermion loops ( = 0)) as anticipated [141], and that instead the series is rendered convergent by extensive cancellations between graphs in perturbation theory.
A gauge invariant set of diagrams contributing to g − 2 consists of all one-particle irreducible graphs with a fixed number n of virtual photons dressing the electron vertex and no virtual fermion loops. For the further case of the dressings of a closed fermion loop, all the diagrams in a gauge set are then encoded in the single integral in Eq. (259) of Appendix C. This "quenched convergence" in the large n photon limit has been addressed previously for loop diagrams within the worldline formulation of sQED and QED using Borel analysis [142][143][144][145].
Finally, we emphasize that in the conventional Wilson loop approach to gauge theories, genuine fermion loop contributions are not included within the framework itself, but are introduced "by hand". As a specific example, the computation of the expectation value of the Wilson loop that appears in cusp anomalous dimensions (which we shall discuss further in Sec. 6) in standard perturbative Wilson loop computations includes only the boson-boson contributions in Eq. (63) but not the boson-fermion and fermion-fermion contributions contained in the worldline computation. These have to be added separately, adding to the complexity of the computation. We will illustrate this point in detail in Paper II by performing an explicit computation of the two-loop cusp anomalous dimension in the worldline approach and comparing it to the prior computation in standard perturbation theory [36].

Resumming the loop expansion
The multiplicity of virtual worldlines in Eq. (22) suggest that the sum in the number of closed loops can be performed formally for any amplitude. Equivalently, Eq. (84) suggests that these loop-resummed amplitudes can also be summed in the number r of real particles. This would then produce a "field-free" worldline generating functional, depending on a external boson source, which makes it suitable for the evaluation of amplitudes using semi-classical expansions.
To perform the sum in the number of loops , we introduce an external gauge boson source B µ (x) that can be regarded as a dummy gauge background field, and define the functional operator This generating functional connects virtual photon lines to all orders in perturbation theory. We can thus rewrite the interaction term in Eq. (23) as With this, the sum in Eq. (22) over all vacuum-vacuum diagrams becomes, where we used the fact that the currents can be decomposed into the sum of individual currents, as in Eq. (19). Summing in and using Eq. (20), we can reexpress this result as where Γ[B] is the one-loop ( = 1) effective action in the presence of the background gauge field B µ (x), One has thus reduced the multi (infinite)-loop problem of the QED vacuum to that of computing the one-loop worldline effective action in an arbitrary dummy background field. Indeed, since Eq. (109) is valid for any choice of B µ (x), one could look for gauge configurations for which one could obtain a non-perturbative result for Γ [B]. Specifically, one can look configurations that minimize the action, allowing one to perform a semi-classical expansion around this value. For instance, for a constant electric field B µ (x) ∝ x µ , Eq. (110) can be expanded in so-called worldline instanton configurations [85][86][87]. Note that this form of the "dummy" background field gives the rate for e + e − pair production in strong external (not dummy) electric fields.
The generating functional technique for amplitudes with real particles is analogous. We can similarly rewrite Eq. (61) as We consider first the case of a single real particle r = 1. The gauge functional average of the dressed 2-point function in Eq. (71) then becomes Since the currents and actions are additive sums of their single-particle counterparts, the path integrals of the r external and the virtual particles can be separated. Using Eq. (58), we get The first quantity in brackets (in the first two lines) is the amplitude for a real fermion to propagate from The quantity in brackets in the third line above is the amplitude for a single virtual fermion loop in the presence of B µ (x). Polarized virtual fermions from the vacuum and real particles are connected by the action of G. Summing this expression to all orders in , The functional average over gauge fields on the l.h.s. is replaced on the r.h.s by the action of G on the product of a single real particle's dressed Green function D F B (x 1 f , x 1 i ) times a one-loop effective action Γ[B], with B being set finally to zero.
To sum in r, we introduce two dummy anti-commuting sourcesη(x) and η(x) and define the quantity From Eqs. (41) and (42) For no real fermionsη(x) = η(x) = 0, we recover Eq. Note also that we factor out the infinite sea of disconnected photon loops Z MW .
In Section 4.1, we discussed in detail the efficiency of worldline calculations order-byorder in perturbation theory. Some comments are now in order regarding possible strategies to extend these calculations to all orders by using expressions like Eqs. (109) and (115) in this subsection.
The formal expressions we obtained as the principal results of Sections 2 and 3 for QED amplitudes require the computation of many-body path integrals describing the interactions between many point-like currents, with their multiplicity growing with the loop order in the expansion. Monte Carlo Euclidean techniques for the computation of these path integrals for scalar theories have been addressed in [88,89], and more recently in [146]. These have thus far only been in the quenched approximation.
Here we have instead reduced the previous quantum many-body problem to the one of solving the quantum dynamics of a single spin-1/2 particle in a given external gauge field, that is later set to zero. This has some advantages. The number of particles is fixed, one must only compute the path integral of a single fermion in an external field and without backreaction, as it is the case of Eq. (109). Further, the dummy field can be chosen conveniently so as to perform a semi-classical expansion of the worldline action around the corresponding instanton configuration.
This worldline action is given by Eq. (278) of Appendix D, where the semi-classical expansion is also discussed at length. We see that replacing the fields sourced by the other particles in the loop expansion, A

IR structure of QED worldline amplitudes
To summarize where we are in the paper, we provided in Sections 2 and 3 compact expressions to compute amplitudes in QED, to all loop orders in perturbation theory, in a formalism in which matter and gauge fields are integrated out explicitly and replaced by super-pairs of 0+1 dimensional worldlines that represent on equal footing real (external) and virtual (loop) particle propertime trajectories in position and spin. In this first-quantized approach to QED, the quantum degrees of freedom of the theory are exactly accounted for as 0+1dimensional boson and Grassmann path integrals over all possible worldline trajectories.
As noted, the computation of the required path integrals in these amplitudes is a difficult task. One approach is an -loop, n-photon, perturbative expansion; another is via a semiclassical expansion about so-called worldline instantons using Monte Carlo methods. The efficiency of each of these approaches, respectively, was discussed in Sections 4.1 and 4.2.
We will now discuss in this framework the problem of the IR safety of the S-matrix outlined in the introduction to the paper. As we will show, the worldline formalism has significant advantages in tackling this problem. The interactions between real and virtual worldlines, to all orders in perturbation theory, are given by the expressions in Eqs. (64)-(67), representing in the IR the long range classical Coulomb dynamics of pairwise interactions of charged worldline super-pairs. This will allow us to examine clearly the IR divergences of the scattering amplitudes, and construct IR safe amplitudes to all orders in perturbation theory.

Curing the IR divergences of the Dyson S-matrix
To study the IR structure of QED, we first consider the Dyson S-matrix element for r positive energy particles 17 in Eq. (84): Recall that W (r, ) (x r f , x r i , θ r , . . . , x 1 f , x 1 i , θ 1 ) encodes the amplitude for the propagation of the r particles from points placed on the celestial sphere at minus infinity to final points placed on the celestial sphere at plus infinity, It also contains virtual fermions attached to the r real particles, and has the structure of a normalized worldline expectation value of the path integrals in Eq. (58) over the r worldline trajectories connecting the points at past and future infinity in Eq. (117), and over the closed worldline trajectories of the virtual fermions. Here S ab denotes the interaction of worldline a with worldline b, and is given by the Wick rotation to Minkowski spacetime of the Lorentz forces in Eq. (63). The indices a, b = 1, . . . , r + run over real and virtual particles.
To explore the IR asymptotics of the S-matrix, and to translate the previous results to the language of Feynman diagrams, we express the Fourier transform of S ab in Eq. (63) as (d = 4 and ζ = 1) where k is the 4-momentum of the virtual photon, and the currents in momentum space and Minkowski time areJ 17 The discussion here can be extended straightforwardly to scattering amplitudes involving particles and antiparticles.
The problem of the IR structure of the S-matrix is therefore reduced to that of examining the low energy behavior (k → 0) of a theory of worldline currents.
In momentum space, a contribution to S ab will be IR divergent when both currents entering Eq. (119) contain 1/k terms that survive when k → 0. In this k → 0 limit, the fermion component of the current in Eq. (120) is sub-leading relative to the bosonic contribution and need not be considered further in this context. For the scalar component of the current, we can sample N points of the worldline x µ (τ κ ) = x µ κ for κ = 1, . . . , N , and connect pairs of points separated by a 4-distance δx µ κ = x µ κ+1 − x µ κ with a straight line path for the propagation of the particle in a time δτ κ = τ κ+1 − τ κ , with τ 1 = 0 and τ N = 1, as shown in Fig. 5 for the worldline of a particle starting at x µ 1 and ending at x µ N . As emphasized before, the worldline has a closed path (x µ N = x µ 1 ) for a virtual particle. For a real particle, these points are external and placed at plus and minus infinity, following Eq. (117). In this discretized form, the path integral over worldline configurations in Eq. (58) reduces to the sum over all possible positions of the intermediate points in the trajectory x µ κ . Figure 5: Discretized representation of a bosonic worldline x µ (τ ) accounting for one possible configuration of a particle's trajectory in the path integral.
In the soft limit k → 0, the total current gives, using Eq. (121), The form of the current above, within a discrete interval, holds the key to understanding the emergence of IR singularities in the Dyson S-matrix and how to cure them. For fixed x µ κ and x µ κ+1 , the term in the parenthesis is of order O(k) and cancels the 1/k in the denominator when k → 0. It follows that a worldline configuration with all points fixed -a particle confined to a finite 4-volume -does not introduce 1/k poles in S ab leading to IR divergences in the S-matrix. However, if the real or virtual particle explores infinity in some step x µ κ → ±∞, which is a contribution accounted for by the path integral, as the one shown in Fig. 6, the corresponding term e ik·xκ in Eq. (122) should be dropped due to the rapid oscillations in the phase in the asymptotic region. The remaining term in parenthesis in Eq. (122) is then left unpaired, and contains a 1/k contribution in k → 0 that introduces an IR divergence in S ab .
x µ κ → ∞ Figure 6: Discretized configuration of a bosonic worldline x µ (τ ) with an internal point at infinity. A dashed line is drawn to indicate that the fermion line extends to infinity and becomes external as a result.
To examine the k → 0 behavior of these particular configurations, we will distinguish between internal and external points placed at infinity. If the κ-th point at infinity is internal, like the one in Fig. 6, it will appear twice in the sum in Eq. (122), which one can then reorganize as Noticing δx µ κ /k·δx κ = δx µ κ−1 /k·δx κ−1 as x µ κ → ∞, and pairing the κ + 1 and κ − 1 terms left unpaired when x µ κ → ∞, this gives which has no 1/k terms when k → 0. Hence contributions within internal intervals of the currents can never introduce 1/k terms leading to IR divergences in S ab , even if the particle trajectory propagates to infinity at intermediate times δx µ κ → ∞. Since for virtual loop particles it is further verified that all points are internal (x µ N = x µ 1 ), it follows from the previous observation that none of the virtual particles in the -th order of the S-matrix loop expansion can introduce 1/k terms in S ab leading to IR divergences. We thus conclude for any virtual fermion i participating in the process lim k→0J µ i (k) = const. + O(k), i = r + 1, . . . , r + .
The only remaining terms potentially introducing 1/k contributions in the k → 0 limit to S ab are the elementary contributions in Eq. (122) containing the external points x µ N = x µ f and x µ 1 = x µ i of the r real particles. A worldline configuration for a real particle is shown in Fig. 7. The external points, unlike internal ones, appear only once in the sum in Eq. (122) Figure 7: Discretized configuration of a bosonic worldline x µ (τ ) with the external points placed at infinity. A dashed line is drawn to indicate that the fermion lines extend to infinity and become external as a result. and are further placed at infinity by the boundary conditions in Eq. (117). Proceeding as previously, reorganizing the sum in Eq. (122) we get In the x 0 N → ∞ and x 0 1 → ∞ limits we then find Then in the k → 0 limit there are two non-vanishing 1/k contributions, encoding the attachment of an IR virtual photon line to the two external legs of the real particle, and where we neglected the phases in the k → 0 limit. Owing to the fact that the bosonic worldline currents are time reparametrization invariant, dividing and multiplying numerator and denominators above by the particle's proper time s and mass m, respectively, produces and where p µ f and p µ i are the worldline momenta of the external legs. We can therefore rewrite Finally, the QED worldline interaction functional has the IR structure (d = 4 and ζ = 1) r+ a,b=1 where the sum runs only over the real particles, contains only the external lines of these real particles, and Λ denotes the IR region in which the previous soft approximations hold. The above statements are nothing but a worldline reformulation, to all orders in perturbation theory, of Low's soft theorem 18 [15].
We turn now to the cure of the IR problem stated above. In standard perturbation theory, the virtual IR divergences are cancelled in the cross-section by real IR photons attached to the squared S-matrix element. As we noted above, the origin of the virtual IR divergences can be directly traced back to dropping the asymptotic currents in Eq. (126), when the x 0 N,n = t f n → ∞ and x 0 1,n = t i n → −∞ limits are imposed by the standard S-matrix construction in Eq. (117). However, as implicit in the discussion of Faddeev and Kulish [29] (building on previous work by Dollard [23], Chung [24] and Kibble [25][26][27][28]), the order of limits matters.
We can therefore, following Faddeev and Kulish, define as an alternative to the Dyson S-matrix, the equivalent worldline FK S-matrix: with the limits t n f,i → ±∞ taken only after all the IR divergences of the diagrams generated by this new S-matrix are cancelled in the k → 0 limit.S (r) f i has the structure of a dressed 2r-point Green function with no truncated legs since time is kept finite until after k → 0. It is therefore IR finite, by construction, to any order in perturbation theory. This is seen clearly by keeping the asymptotic charged currents that consist of terms entering with phases k·x N and k·x 1 in Eq.
18 Tests of Low's soft theorem by several experimental analyses to date show excellent agreement between Low's prediction and soft photon data in QED processes. On the other hand, a significant excess of soft photons is observed for a wide variety of processes involving hadron final states; these so-called "anomalous soft photons", provide a phenomenological handle to study the interplay of QED and QCD sources of photon production. For a recent discussion, see [147]. instead of Eq. (131). The three terms above are all of order O(1/k). However, after taking the k → 0 limit, they cancel exactly amongst each other and there are no 1/k contributions left in S ab . We conclude therefore that there are no IR divergences in the FK S-matrixS (r) f i , to all orders in perturbation theory.

An example: Möller scattering
We will now discuss as a specific example of the previous general discussion, Möller scattering (e − e − → e − e − ), which is illustrated in Fig. 8. We will first analyze the IR divergences arising from the standard Dyson S-matrix, and subsequently contrast the resulting expressions with those of the FK S-Matrix. We begin by writing the pair-wise worldline interactions as where the soft and hard contributions are defined as where Λ delimits the separation between the IR and UV regions. We first consider the Dyson S-Matrix. S ab UV is IR finite by construction. Following the previous discussion, the 1/k terms responsible for IR divergences in S ab IR are introduced by the t f → +∞ and t i → −∞ legs of the scalar components of the worldline currents of the two real electrons. The virtual fermion contributions, as we discussed previously, are IR finite. Thus in the infrared, the sum over all real and virtual worldlines to any loop order , can be approximated by a sum over the real worldlines 19 : Further, for sufficiently small Λ, virtual photons in S nm IR attach to the asymptotic regions of the worldlines corresponding to t f → +∞ and t i → −∞. In these limits, the weakly interacting worldlines approach asymptotically the trajectories of free particles inJ n µ,IR (or equivalently,J m µ,IR ), with initial and final velocities β n i and β n f , As shown in Fig. 8, any arbitrary configuration of worldlines at short distances, indicated by the cusps at x n c , are accounted for in S ab UV . Using Eq. (138), the scalar currents in the IR region split into two parts, where we introduced an i prescription to kill the incoming and outgoing currents at infinity. Performing the time integrals, and taking t n i → −∞, t n f → +∞, this gives For sufficiently low Λ, the phase in the above expression, of order 1 + O(k), can be replaced by unity. The interactions in the IR region are then given by Because for weak long distance interactions (represented by the above expression in the asymptotics Λ → 0) the n and m worldline paths are well approximated by free classical trajectories, they can be factorized out of the path integral in Eq. (118). This gives where we defined W (2, ) UV (· · · ) as the normalized worldline interaction functional of two real particles (r = 2) with virtual fermions and an arbitrary number of virtual photons with energies above Λ. Introducing this expression in the Dyson S-matrix in Eq. (116), one obtains Note that the S-matrix element for virtual photons with k ≥ Λ is IR finite by construction and given by Eq. (143) corresponds to Weinberg's exponentiation of virtual Abelian IR divergences [5]. It depends only on the directions of the charges at infinity, is independent of their spin dimensions, and encodes their interactions at late times with their self-generated classical gauge fields. This IR factor dresses the hard S-matrix S (2) f i,UV containing photons with k ≥ Λ outside the IR region. Note that this separation of scales is also assumed implicitly in [5]. To evaluate the diagrams with these IR divergent dressings we keep the phases of the currents in Eq. (140), depending on the cusp positions, and finally take the soft photon k → 0 limits. The diagrams in Eq. (141) are shown in Fig. 9 and denoted where S(β n , x n , β m , x m ) = g 2 2 (146) and η n = +1 or η n = −1 for a charged current incoming to or outgoing from the point x n .
The diagonal terms n = m correspond to IR self-energy dressings of the S-matrix occurring at asymptotic times. The non-diagonal n = m, to virtual IR photons exchanged by the two electrons. Consider the initial-initial photon exchange diagram, for which η n = +1 and η m = +1, and given by Defining ω 2 k = k 2 , the integration in the virtual photon energy k 0 , as shown in Fig. 10, has poles at The integration contour is determined by the sign of t n c − t m c . For t n c > t m c , the clockwise contour C + in Fig. 10 should be chosen. The residue of the pole at k 0 = +ω k − i encodes the interaction of the electron n in the radiative mode of the gauge field of the electron m, created by the overall change in the in and out directions. The residue of the pole at k 0 = kβ n i − i corresponds to the Coulomb mode in the interaction of one electron in the Liénard-Wiechert field of the other. This gives where we introduced an IR cut-off λ to capture the logarithmic IR divergence in the resulting integrals.
For t n c < t m c , the anti-clockwise integration contour C − should be chosen instead, with the residues corresponding now to the response of the electron m in the radiative and Coulomb fields of n. In the IR k → 0 limit, the long distance interactions accounted for in S nm IR are however independent of the relative placement of the cusps of Fig. 8, and more generally, of the particular short distance configurations of the actual electron worldlines contributing to the path integral.
Following the same procedure for the rest of diagrams, it can be shown that the IR divergent, long distance interactions between any two particles n and m in S nm IR can be always split, following Weinberg [5] and Faddeev and Kulish [29], into a real and an imaginary contribution where and Performing then the remaining integrals, using Eqs.
with the cusp angles between external legs of the charged-particle worldlines defined by Recall further that η = +1 and η = −1 for an incoming and outgoing external chargedparticle leg. The indices n , m in the r.h.s. indicate an in or an out leg of the real worldlines, and the pairs (n , m ) and (m , n ) are counted separately. The prime in the imaginary part sum indicates that the sum runs only over the out-out and in-in pairings [5] since the Liénard-Wiechert or Coulomb modes in all crossed diagrams, namely S(β n i , x n c , β m f , x m c ) and S(β n f , x n c , β m i , x m c ), vanish in the soft limit.
We turn now to analyze the IR behavior of the FK S-matrixS (r) f i in Eq. (133). The currents at long distances are well approximated again by the cusped worldlines of Fig. (8). Using Eq. (138) and keeping t n f,i finite, the currents givẽ As observed previously in Section 5.1, in the standard approach [5], the first and the last terms in Eq. (155) are dropped due to the rapid oscillations in the phases occurring at asymptotic times. This would give back Eq. (140). As implicit in Faddeev and Kulish's discussion on infrared safety [29] the order of limits matters. For any fixed large x f,i , the phases go to unity as k → 0, and the asymptotic contributions cancel with that from the cusps in this limit.
In the standard construction, these new terms correspond to the (p ± k) 2 = m 2 poles not accounted for by the LSZ truncation of the amplitude, order-by-order in perturbation theory. The asymptotic currents encode virtual photons pinched to the x µ f and x µ i points, and they can be regarded then as IR dressings of the in and out states in the gauge field to all orders in perturbation theory.
With the aid of the asymptotic currents in Eq. (155), the new Lorentz-covariant worldline interaction functional includes photon exchanges at infinity that capture and remove IR singularities to any order in perturbation theory. This can be confirmed explicitly by plugging Eq. (155) into Eq. (136), and finding it to be free of IR singularities in d = 4. Indeed, in the FK S-matrixS where S nm IR are the standard diagrams in Eq. (145) and the new asymptotic diagrams S nm as include a virtual photon attached in at least one of its ends to one of the asymptotic boundaries (t i → −∞ and t f → +∞) where each diagram S(β n , x n , β m , x m ) is also given by Eq. (146). A subset of these asymptotic diagrams are displayed in Figs. 11 and 12. The evaluation of the asymptotic diagrams, follows the same procedure as for the diagrams discussed previously, and also lead to real and imaginary contributions from the asymptotic regions, −iS nm as = R nm as + iΦ nm as (158)  with R nm as = where Λ denotes the IR region of the asymptotic interactions, and Performing the remaining integrals using Eqs.
which is independent of the IR cut-off λ. Our result in Eq. (161) therefore proves that the FK S-matrixS (r) f i in QED is IR finite to all orders in perturbation theory. A few remarks on the cancellation of the IR divergences in Eq. (161) are now in order. Firstly, because the cancellation of these IR divergences in the conventional Dyson S-matrix approach occurs between virtual and real photon IR divergences in the cross-section, one might ask if virtual IR photon diagrams becoming IR finite when asymptotic interactions are included means that real IR divergences they canceled previously now survive in the cross-section. This is fortunately not the case because diagrams with the emission of real photons are also accompanied by corresponding asymptotic diagrams, which similarly lead to IR finite cross-sections. Real photon emissions will be considered in detail in Paper II.
Secondly, in the Dyson S-matrix the imaginary Coulomb phase iΦ exists at the amplitude level, and becomes manifest in the calculation of the Dyson S-matrix element itself starting at next-to-leading order in the coupling g. It should therefore not be ignored in the context of our previous discussion of the IR finite FK S-matrix 20 .
Further, we have included the n = m terms in both in the conventional IR divergent Dyson and IR finite FK S-matrices. These Coulomb IR divergences in external lines (wherein the photon is attached to the same external charged-particle worldline) are rarely discussed because it is presumed that they are included in the counterterms corresponding to wavefunction renormalization of the external currents. However as noted by Weinberg [83], this is not the case because the presumed counterterms are themselves cancelled by counterterms arising from the internal lines and vertices. Therefore these n = m IR Coulomb terms that are manifest in the worldline formalism should be kept. In the conventional approach, their IR divergences will cancel with the IR divergences of the real photons attached to the cross-section; in the FK approach, their IR divergences are removed by the corresponding counterterms, as is manifest in the result obtained in Eq. (161).
Finally, as also previously observed in [36], there is a subtle point regarding the expressions surviving the cancellation of the IR divergences in Eq. (161). The Coulomb phase is ill-defined for the cases in which γ n m = 0. For the case of the long distance interactions between different charged particles (n = m ), this would correspond to the scenario in which the relative velocity of a charged-particle in the rest frame of the another, defined by the invariant relative rapidity v(β n , β m ) = 1 − β 2 n β 2 m /(β n ·β m ) 2 , vanishes in the in and/or out state, reflecting a singularity of the scattering amplitude. For photons attached to the same external line of a charged-particle, it is always true that γ n n = 0. This divergence corresponds to the interaction of a charge in uniform motion in its own Liénard-Wiechert field. In the Lorentz frame of the particle at rest, this is just the divergence of a Coulomb field at short distances and will not impact scattering matrices at any finite relative velocity.

Relating asymptotic worldlines to Faddeev-Kulish soft factors
We will now relate our construction of IR finite amplitudes in the worldline formalism to the Faddeev-Kulish S-matrix in their original Hamiltonian approach [29], and to the recent discussion in the modern language of Wilson lines by Schwartz and Hannesdottir [32]. The fundamental issue noted in [29] is that to construct an IR finite S-matrix operator, the evolution of the states in the far past and future have to account for contributions from long-range interactions that are neglected in the Dyson S-matrix.
We begin with the discussion in [32], where where i∂ t U as (t, t ) = H as (t)U as (t, t ) and H as (t) is the QED Hamiltonian in the asymptotic soft photon limit. In the traditional construction of the Dyson S-matrix, it is assumed that this asymptotic dynamics is well approximated by U 0 (t, t ), with i∂ t U 0 (t, t ) = H 0 (t)U 0 (t, t ), with H 0 (t) being the free QED Hamiltonian. The previous equation can be rewritten as The operator relating the out and in states defines the FK S-matrix: where the relations U † (t i , t) = U(t, t i ) and U(t f , t)U(t, t i ) = U(t f , t i ) were employed. Equivalently, for the Dyson S-matrix, one must replace above U as (t, t ) with U 0 (t, t ): Thus one finds the former can be expressed in terms of the latter as In the FK S-matrix (referred to as the "Hard" S-matrix in [32]), an in state is dressed by asymptotic interactions from t to t i = −∞, subsequently evolves and scatters in the central region from t i = −∞ to t and from t to t f = +∞ as per the Dyson prescription, and returns to t dressed analogously by final asymptotic interactions, with t = 0 chosen in [32]. Feynman rules for the leading IR divergent diagrams including contributions from the asymptotic regions, are shown explicitly to cancel with the IR divergences generated by the Feynman graphs corresponding to the Dyson S-matrix in the central region.
In the original FK S-matrix discussion in [29], the asymptotic amplitudes likewise enter as conjugated to the ones in the central region, but with phases at t f → +∞ and t i → −∞ instead of t = 0 in [32]. Since the role of the asymptotic interactions is to transform a free state at t = 0 to a dressed in and out states at t f → +∞ and t i → −∞, respectively, the matrix elements of the Hard S-matrix [32] are equivalent to the FK form in [29].
The expression in Eq. (166) was further reexpressed by Faddeev and Kulish as where the real and imaginary parts of the IR divergences of the diagrams from the asymptotic regions are encoded in the action of the operators and on the in and out states. Here is the fermion number operator of the charges in the in and out states, b s p (d s p ) and b s, † p (d s, † p ), annihilation and creation operators respectively, of a fermion (antifermion) with momentum p and spin s, a λ k and a λ, † k are annihilation and creation operators of a photon of momentum k and polarization λ. The properties of these operators were considered in some detail in the original Faddeev and Kulish paper, where R is a coherent state operator spanning the Hilbert space of asymptotic charged particle states. However as noted in [29], and further emphasized in [32], these states, dressed with clouds of soft photons, can never be thought of as direct products of similarly dressed single particle states. As also noted previously, the phase Φ(t) is the relativistic generalization of the Coulomb (or Eikonal) phase 21 , which is crucial in canceling the infinite phase factors that arise in the Dyson S-matrix. In our worldline formulation of scattering amplitudes, we integrated out the gauge and matter fields explicitly, giving rise to non-local interaction functionals quadratic in the fermion currents, as in Eq. (119). This prevents us from defining a FK S-matrix in the same way as in Eq. (167), where the asymptotic dressings are linear in the gauge field, interactions are local, and currents can be factorized out from interactions in the central region to dress either the traditional S-matrix, or the in and out states, as desired.
In the operator formulation however, it is clear that the action of R(t) -more generally, the action of R(t) in conjunction with terms coming from the central region in the Dyson S-matrix S, order-by-order in perturbation theory -will introduce the radiative modes of the asymptotic interactions that, once exponentiated, will lead precisely to the R nm as terms we found in Eq. (159). Likewise, the action of Φ(t) on the asymptotic Coulomb modes, once exponentiated, will lead to the Φ nm as contributions we found in Eq. (160). It is worth noting that each charged worldline in this framework is not separately dressed into an individual cloud, but this cloud depends on the rest of charges: the dressings involve photon exchanges between different charges in the asymptotic region, and also photon exchanges between the central and the asymptotic regions. As discussed previously, these features are crucial to completely cancelling the IR divergences inS. Thus the identification of charged worldlines with charged particles, while a useful mnemonic, is not fully robust.
Our form of the IR finite FK S-matrixS presents potential advantages. In the usual FK formulation [29,32], the S-matrix elements are only defined perturbatively, and the IR finiteness of amplitudes can only be proven order-by-order in perturbation theory. In the worldline approach, the amplitudes are defined to all orders in perturbation theory, with an arbitrary number of virtual photon and fermion loops attached. In particular, we showed that the exponentiation of IR divergences obtained from perturbation theory [5] appears naturally, with the added intuition of a theory of classical currents of nearly free first-quantized particles as they approach the far past and future. More importantly, the emergence of the IR dressings of the central and asymptotic regions as non-local worldline interaction functionals between classical electromagnetic currents is transparent. This allowed us to derive an all order proof of the IR finiteness of the FK S-matrix of QED.
Another potential advantage presented by the worldline framework, as discussed in Section 4, is its efficiency in high order perturbative calculations. Each order in perturbation theory can be generated in the worldline from a single one-loop ( = 1) master expression, which contains n-factorial combinations of Feynman graphs corresponding to different permutations of photon insertions into the fermion lines. We will discuss this further in the next section in the context of perturbative computations of the anomalous dimension of a cusped Wilson loop.
An interesting recent development is the interpretation of the asymptotic interactions encoded in the FK S-matrix in the language of symmetries. It was pointed out very early by Weinberg [152] that the universal features of soft theorems in Abelian theories [5] might arise directly from conservation laws of some fundamental symmetry. In his interpretation, the form of soft theorems is dictated from Lorentz invariance to simply enforce charge conservation and the principle of equivalence. More recently, Low's theorem has been shown to be equivalent to the Ward identity [58][59][60] of an asymptotic symmetry of the amplitudes generated by the group of large U (1) gauge transformations acting over the in and out states with angle-dependent constants at null infinity [54]. In gravity, Weinberg's soft graviton theorem can be understood as the Ward identity [54,55] of the BMS group of asymptotically flat spacetime symmetries. [51,52] The invariance of amplitudes under large gauge transformations implies that any scattering process is accompanied by a transition between degenerate vacua containing a differing number of IR photons. These shifts between different vacua have been found [60] to be implemented precisely by the same asymptotic interactions encoded in the dressings of the Dyson S-matrix proposed by Faddev and Kulish in [29]. In the conventional construction of amplitudes through the Dyson S-matrix, the in and out states are assumed free chargedparticle states, and then the vacuum of the theory is assumed to be unique. Because this would violate the conservation law, the Dyson S-matrix has to be identically zero. The emergence of IR divergences in the form of Low's soft photon theorem, order by order in perturbation theory, can be therefore interpreted as a necessary condition to avoid spoiling the BMS-charge conservation laws associated to the BMS-like symmetry of the theory.
As emphasized before, since we have integrated out the gauge field completely, the FK S-matrix in the worldline form of Eq. (133) prevents us from disentangling the dressings in the gauge field of the charged-particle states in the in and out states at plus and minus infinity, from the interactions of these charges with the gauge field occurring at intermediate times in the central region. This indeed has to be the case, since the interaction of these charges to all orders in perturbation theory is encoded in a non-local interaction functional of all charged-particle currents. In this worldline framework then, the vacuum transitions mentioned previously are automatically implemented within the definition of the scattering amplitude itself, by means of the asymptotic currents corresponding to the terms with phases k·x N and k·x 1 in Eq. (126). We plan to return to the discussion of asymptotic symmetries in the worldline context in future work.

Cusp anomalous dimension from worldlines
The fundamental objects in the worldline formulation of QED are normalized expectation values of non-local interaction functionals of charged currents created by super-pairs of scalar x µ (τ ) and Grassmann ψ µ (τ ) worldlines. These expectation values are expressed as generalized Wilson loops, containing the non-local dynamics of point-like super-pairs of particles. Worldlines therefore provide a natural framework to examine the renormalization properties of Wilson loops in QED 22 . In particular, we will show how an efficient calculation of the cusp anomalous dimension can be carried out in the worldline formulation, with all the diagrams generated, to any order in perturbation theory, from the master formula discussed in Section 4.
To dress a cusped worldline to all orders in perturbation theory, we consider the interactions of a single real charged-particle (r = 1) and an arbitrary number of virtual fermions . This reduces to evaluating normalized worldline expectation values of interaction functionals of their respective currents, which are given by Eq. (55) to be The asymptotic behaviour of a real charged-particle is dominated by its scalar degrees of freedom, and the evolution of its worldline x µ (t) is well approximated by a classical trajectory x cl µ (t) possessing a cusp. Following Eq. (50), this is where d = 4 − 2 and we absorbed the mass dimension of the coupling in a parameter µ with mass dimension unity, and where, as previously, β i and β f are the initial and final velocities. The cusp is placed at x cl µ (0) = 0 without loss of generality. The dressings of this scalar charged-particle's Wilson line in the Abelian gauge field are encoded, to all orders in perturbation theory, in the dressed 2-point function of Eq. (71), with the spin degrees of freedom and the path integration over the real worldline configurations omitted. We define the amplitude W accordingly as where x f = x cl (t f ) and x i = x cl (t i ) as t f,i → ±∞, and the th contribution in the loop expansion is given by the normalized worldline expectation value taken only over the worldline configurations of the virtual fermions Further expanding Eq. (173) in powers of g, precisely in the same way as discussed in Sect. 4.1 for the vacuum-vacuum amplitude, gives the series Here the disconnected fermion loop graphs are removed by the corresponding terms Z ( ) (n) in the loop expansion of the vacuum-vacuum amplitude. The resulting n th and th term in the power expansion is of order α n and contains n virtual photons and virtual fermions attached to the cusped charged-particle worldline.
To one-loop order (n = 1 and = 0), Eq. (175) gives where we used the photon propagator in Eq. (17) in d = 4 − 2 dimensions and Minkowski signature, This one-loop contribution contains UV divergences from virtual photons attached to equal points in the worldline x cl (t 1 ) = x cl (t 2 ), and IR divergences coming from the virtual photons attached to the far past and future δ → ∞. Using Eq. (172) the integral can be split into four contributions, their corresponding diagrams shown in Fig. 13, giving where the required integrals I(β f , β i ) are defined as The Schwinger parameters t = t 1 +t 2 and s = t 1 /t are introduced in the second equality. The scale-less t integral receives both the t → 0 (UV) divergence and the t → ∞ (IR) divergence from photons attached to equal and very distant points in the worldline, respectively, in Eq. (176). To capture the IR divergence, we regulate the integral over t as The s integral remains finite. Setting = 0 then yields the geometrical factor where one identifies the 4-dimensional cusp angle between the final f and initial i directions as γ f i , as shown previously Fig. 13, allowing us to write Eq. (181) as The diagrams with initial-initial and final-final photon loops are given by the expressions −I(−β i , β i ) and −I(β f , −β f ), corresponding to the substitutions γ f f = 0 and γ ii = 0 in the equation above 23 . Collecting all the terms in in Eq. (176), and expanding around = 0, gives where γ E is the Euler-Mascheroni constant. Plugging all these results into Eq. (176) gives, to one-loop, Eq. (185) contains both the cusp IR ∼ log(1/λ) and UV ∼ 1/ divergences. The former corresponds precisely to the self-energy contributions (n = m) in the IR dressings of a scattering amplitude found in Eq. (153). It can therefore be cancelled following the Faddeev and Kulish procedure discussed in Section 5.2, leading to Eq. (161).
To renormalize the UV divergences in W, we follow the procedure in [154]. In Eq. (185), the prefactor of the UV divergence 1/ has a dependence on the cusp angle related to the short-time integration in the vicinity of the cusp. The required renormalization factor Z(γ f i ; g, ) is then also angle dependent. This suggests, following [154], where g R is the renormalized coupling and the renormalized amplitude W R (encoding the dressings of the classical cusped scalar worldline) satisfies the renormalization group equation, (RGE) with the anomalous dimension Γ(γ f i , g R ) introduced by the cusp. The QED β-function is defined by (d = 4 − 2 ) where α = g 2 R /4π. Expanding Z in powers of g R and matching, to subtract to one-loop order, the 1/2 UV divergence of W gives Hence the one-loop cusp anomalous dimension is then The previous computation can be extended to higher orders. We first note that all treelevel dressings ( = 0) in Eq. (175) exponentiate, as anticipated, by setting = 0 in Eq. (173) and avoiding therefore the power expansion in Eq. (175). One thus obtains, where the ellipses indicate that the = 0 contributions coming from the attachment of virtual fermion loops remain to be included. The corresponding graphs for these in W are generated in exactly the same way, as the = 0 term by the compact expression in Eq. (175). The building blocks required of this computation, for arbitrary n and , are given by the multi-loop worldline polarization tensors of Eq. (259) in Appendix C.
For instance, at two-loop order (n = 2 and = 1), the first finite correction to the tree level graphs in Eq. (191) using Eq. (175) contains as its building block the expression 1 2 where the polarization tensor Π µν (k 1 , k 2 ) is given by Eq. (259), as shown in Appendix C. As we discussed in Sec. 4.1, automatically encoded in a single and compact expression are multiloop diagrams with many virtual fermions contributing to the cusp anomalous dimension at higher orders. In contrast, in the standard computations of the two-loop cusp anomalous dimension [155], the exact form of a contribution like the one in Eq. (192) has to be obtained separately from the corresponding Feynman graphs that in the UV limit include the spin structure of fermion loops. The efficient computation of higher order loop contributions to the cusp anomalous dimension in the worldline formulation of QED will be addressed in detail in Paper II.

Conclusions and outlook
We developed in this paper a reformulation of QED as a fully covariant, first-quantized, many-body theory of 0+1 dimensional super-pairs of bosonic and Grassmann worldlines. In particular, the QED path integral can be expressed as a power series expansion in generalized Wilson loops, where the Wilson loop at a given -th order (corresponding to virtual worldlines describing loops in propertime) contains a path integral over all possible paths of pointlike boson and fermion super-pairs. The dynamics of the virtual fermions in this path integral is encoded in the exponentiated photon exchanges between the loop currents that are functionals of the boson and fermion super-pairs. We further extended our formalism to construct scattering amplitudes for charged wordlines interacting via the exchange of arbitrary numbers of virtual photons and charged fermion loops. One can similarly express the amplitude for the scattering of r worldlines to -th loop order in terms of a generalized Wilson line/loop structure corresponding to the path integral over first-quantized, many-body worldline configurations of the r real worldlines and virtual worldlines. As previously, the dynamics of their interactions is contained in the exponentiated photon exchanges between the real and virtual currents, with additional energy-momentum and helicity-momentum constraints imposed on the real worldlines.
We outlined the usefulness of this formalism for high order perturbative computations by showing that terms to arbitrary order in a loop expansion of the vacuum-vacuum amplitude can be generated from a compact expression for multi-loop polarization tensors. This master formula can further be applied to efficiently perform higher order calculations of any given QED process, automatically encoding the factorial growth of diagrams in conventional perturbation theory. We also demonstrated that all-order worldline results can be extracted from the generating functional for the QED one-loop worldline effective action in an background field. This suggests that the non-perturbative structure of worldline amplitudes (an example being Schwinger pair production in strong electromagnetic fields) can alternatively be explored by the use of semi-classical expansions in appropriate background fields.
We also showed that the worldline formalism provides a powerful framework to explore the infrared structure of QED. By examining the classical long-range interactions of the charged-particle worldline currents, we showed that virtual IR divergences and Abelian exponentiation theorems naturally emerge in the worldline expressions as soft dressings of the Dyson S-matrix. We showed that these soft dressings that set the Dyson S-matrix to zero can be subtracted by interaction terms with the exchange of soft photons between charged currents, with at least of one them in the asymptotic region of very early or late times governed by the Coulomb dynamics of classical currents. This observation allowed us to construct a Faddeev-Kulish worldline S-matrix including these asymptotic contributions and prove it to be free of IR singularities to all orders in perturbation theory. We illustrated our general argument for the more specialized case of Möller scattering and discussed the mapping of our worldline results to the original Faddeev-Kulish results and more recent discussions in the literature.
We finally addressed the renormalization of cusped Wilson lines which naturally correspond to worldline cusps in our framework. We showed how the one-loop cusp anomalous dimension is recovered in this framework and outlined how this formalism enables efficient computation of cusp anomalous dimensions to higher orders.
Indeed in Paper II, we will discuss an explicit computation of the QED two loop cusp anomalous dimension and compare and contrast this with previous computations in the literature. We will also discuss in Paper II the extension of the computation of scattering amplitudes with all order virtual exchanges discussed here to that of amplitudes involving the emission of real photons. This will complete our proof of the infrared safety of the worldline Faddeev-Kulish S-matrix.
It would also be interesting to understand if the infrared behavior of the worldline Faddeev-Kulish S-matrix can be recast in the universal language of asymptotic symmetries of gauge theories. In particular, since QED can be exactly rewritten in this first-quantized albeit non-local formalism, how does one interpret the non-local FK dressings as Goldstone bosons of spontaneously broken asymptotic symmetries? We plan to pursue such questions in follow-up work.
A further application of our work is to explore whether worldline Monte Carlo methods can be employed specifically to compute the -th loop, n photon polarization tensors that we identified as building blocks of perturbative computations in this framework. If so, this would have immediately applications in precision QED tests for a wide range of observables.
Further important extensions of our work are to QCD and gravity, where semi-classical worldline formulations have been known for decades. In the case of QCD, we noted that the powerful string amplitude and worldline formalisms have been shown to be equivalent. However it is not known whether the procedure followed in this paper can be extended to QCD with the necessary systematic approximations. This problem has a correspondence in the well known language of non-Abelian exponentiation. A specific goal would be to extend our methods to compute high order cusp anomalous dimensions in QED to QCD. Another would be to apply these ideas to the Regge regime, where worldline computations in the presence of semi-classical background fields is promising both in precision computations and in achieving deeper insight into infrared phenomena such as color memory [156,157]. With regard to the former, as we noted briefly, plans are underway to study puzzling apparent violations of Low's theorem at the LHC.
Finally, looking towards the longer term, the fact that QED can be formulated in full generality as a many-body theory of first-quantized worldlines may have powerful implications for quantum computing. Since Grassmann worldlines are natural qubits, it would be interesting to explore whether the supersymmetry between Grassmann and boson superpairs can be exploited to simplify digitization of the large Hilbert space occupied by bosonic worldlines. also like to thank Niklas Mueller for previous collaborative work on worldlines, as well as Fiorenzo Bastianelli, James Edwards, Christian Schubert and the "worldliners" seminar group for useful discussions on this topic. with η µν = diag(+1, +1, +1, +1), µ = 1, 2, 3, 4. It follows that, from the relations above, and and F E µν an antisymmetric tensor. While not affected by the Wick rotation, it will be convenient for us to define Hermitian Euclidean Dirac matrices (γ E µ ) † = γ E µ and (γ E 5 ) † = γ E 5 -so that (γ E µ ) 2 = 1, (γ E 5 ) 2 = 1 and the squares of Dirac operators are positive operators, where / p E = p E µ γ E µ . This requires For the purposes of the discussion in the main text, it will be also convenient for us to introduce the Hermitian matrices One can verify that {γ E 5 ,γ E µ } = 0 and {γ E µ ,γ E ν } = 2η µν . Since (γ E 5 ) 2 = 1, the algebra can be completed by The Euclidean representation of spinors, and their Wick rotataion is subtle; for a nice discussion, see [158]. Since we are primarily interested in Minkowski scattering amplitudes, we will employ for spin the definitions It follows that In the rest of the paper, we will drop the superscripts in fields and coordinates, unless explicitly required. All amplitudes in physical time or imaginary time can be recovered by the replacements in Eq. (196), Eq. (197) and Eq. (199). In particular, the starting point for our discussion, the Euclidean vacuum-vacuum amplitude in Eq. (1) can be obtained in this way from Eq. (193) and setting the external sources to zero.

B Worldline path integral for the dressed fermion propagator
We will present here the derivation of the dressed spin-1/2 propagator in the worldline formalism. This approach can be straightforwardly extended to higher-point correlators.
We will follow the logic of the previous computation by Fradkin and Gitman [125] but shall further clarify some subtle features as well as adapt the derivation to our conventions.
We begin with the equations of motion of the Euclidean dressed fermion propagator in the gauge background field, and seek to represent D A F (x, y) as a sum over worldline trajectories connecting x and y. Following [125], to explicitly account for the helicity-momentum constraint that the worldlines must satisfy, one considers alternativelyD A F (x, y) = D A F (x, y)γ 5 , with γ 5 = γ 0 γ 1 γ 2 γ 3 . Multiplying both sides of Eq. (204) with iγ 5 , this modified propagator satisfies It will be also convenient for us to define the Hermitian matricesγ µ = +iγ 5 γ µ and complete the algebra {γ λ ,γ κ } = 2 η λκ with λ = µ, 5 andγ 5 = γ 5 . (Eqs. (200) and (201) of Appendix A.) Inverting Eq. (205), the problem reduces to computingD A F (x, y) = x|D A F |y with the operator defined as The operator has been decomposed into a bosonic piece (γ µ D µ + iγ 5 m) 2 and a fermionic pieceγ µ D µ + iγ 5 m. The former can be reexpressed as the integral on the r.h.s in terms of the commuting propertime variable ε. The latter can be represented, as shown, by the Grassmann integral in the anti-commuting propertime variable χ. Identifyingp µ = i∂ µ , the bosonic contribution to the propagator yields Further, since the quantity above is bosonic, we can interchange the order of integrands and write the dressed propagator as where the Hamiltonian operator including the propertime super-pair (ε, χ) is given bŷ The first bracket represents the mass-shell constraint satisfied by the worldlines, which becomes apparent after anti-symmetrizing iγ µγν into the spin tensorσ µν , while the second bracket corresponds to the helicity-momentum constraint 24 .
To obtain the path integral representation of the dressed propagator, one needs to discretize this matrix element in n steps. We will however first generalize the parameters ε and χ to commuting ε(τ ) and anti-commuting χ(τ ) functions of the worldline parameter τ , on equal footing with x µ (τ ) and p µ (τ ). This implies thatε andχ can be represented by operators (albeit with trivial dynamics) acting on states, just as forx µ andp µ . We definê π as the operator conjugate toε, andν as the corresponding conjugate ofχ. With these modifications, we obtain This then allows us to write the expectation value as with the boundary conditions ε n = ε and χ n = χ. Upon the trivial action of the Hamiltonian on the states, and using Eq. (210), we get with δa k = a k+1 − a k , and we divided and multiplied by 1/n = δτ k to take the continuum limit. By integrating in π k and ν k we obtain n commuting and n anti-commuting Dirac deltas that set all the ε k = ε and χ k = χ to the constant conventional Schwinger parameters, which can be understood to be the zero-modes of these einbeins. Taking δτ k → 0, the gauge-unfixed path integral in Eq. (208) is given bȳ This exponentiation is still incomplete since the path integral involves the time ordered exponential of an operatorĤ containing Dirac matrices. To get rid of the time ordering and exponentiate spin dynamics, Fradkin and Gitman suggested introducing five auxiliary anticommuting sources {ρ λ (τ ), ρ κ (τ )} = 0, which also anti-commute withγ λ . Then Eq. (213) can be rewritten as Now using the Baker-Campbell-Hausdorff relation recursively, we get the following Magnus expansion: where the ellipses denote a series of higher order nested commutators. However since forms a bosonic subalgebra, it commutes with ρ λ (τ )γ λ . The higher order nested commutators in Eq. (215) therefore vanish and we get T exp Symmetrizing the double integral to make τ 2 run from 0 to 1, and using the fact that {ρ λ (τ 2 ), ρ λ (τ 1 )} = 0, we get for the second term on the r.h.s. of Eq. (217): Noting that τ 2 |(d/dτ ) −1 |τ 1 = 1 2 sign(τ 2 − τ 1 ), the above scalar product of ρ λ (t) can be cast as into a path integral over an anti-commuting worldline ζ λ (t) with the normalization 25 given by the path integral in absence of sources (ρ λ (τ ) = 0): Anti-periodic boundary conditions in the Grassmanian path integral are assumed (ζ λ (1) = −ζ λ (0)) to make it invariant under shifts of the worldlines. The other term on the r.h.s. of Eq. (217) can be written as where θ λ , analogously to the ρ λ , are five anti-commuting auxiliary constant sources anticommuting with theγ λ matrices too. Hence after these manipulations, we obtain for the time ordered exponential Eq. (217) (which appears in the dressed propagator in Eq. (214)) the result .
The form of the interaction term suggests that we can shift the path of integration by ψ λ (τ ) = ζ λ (τ ) + θ λ . Since ζ λ (1) = −ζ λ (0), we then get ψ λ (1) + ψ λ (0) = 2θ λ . Using this, and the commutation relations, we get Hence in these new variables, Eq. (222) becomes Finally, inserting Eq. (224) into Eq. (214), we obtain our result for the dressed propagator, This expression provides the exact worldline Hamiltonian representation of the dressed fermion Green function in the background of an arbitrary Abelian gauge field A µ (x). Shifting p µ (τ ) + gA µ (x(τ )) = p µ (τ ), and performing the path integral in the momentum worldline, we get the equivalent Lagrangian form, where the reparametrization and super-gauge invariant Euclidean action S = 1 0 dτ L(τ ) of the relativistic spinning particle is given by the Lagrangian [69,72] To obtain the classical equations of motion of a single spinning charge in real time, we now need the Green function in Minkowski spacetime. Implementing in the expression above the Wick rotations in Eq. (196) and Eq. (199), and replacing the propertime einbein ε(τ ) → iε(τ ) and its conjugate π(τ ) → −iπ(τ ), we get with µ = 0, 1, 2, 3 now. The trace over the spin degrees of freedom can be expressed as since θ 2 = 0 and N 4 = 2 2 . Setting now θ µ = 0 in what remains produces in the boundary conditions ψ µ (1) + ψ µ (0) = 2θ µ = 0; hence the boundary term in the action ψ µ (1)ψ µ (0) = −ψ µ (1)ψ µ (1) = 0 cancels. Finally taking the remaining trace over spatial positions, we get We emphasize that in the normalized open worldline expectation value open boundary conditions hold and the Grassmanian worldline has 5 components corresponding to the action of the 5 Dirac matrices, ψ λ ↔γ λ . In contrast, in the normalized closed worldline expectation value of the fermion loop determinant above, the Grassmann worldline has 4 components ψ µ ↔ γ µ . Integrating over the einbein variables leads to the well-known reparametrization gauge-fixed closed worldline form for the fermion loop determinant [76] derived employing a wide variety of very different methods in the literature.
The structure of these tensors admits an easy generalization to higher orders since in this formalism the currents for each new photon vertex are additive. Thus to any n-th order in perturbation theory, the corresponding polarization tensor of rank n is given by Π µ 1 ...µn (k 1 , . . . , k n ) = −2 g n n! (2π) 4 δ (4) which (up to differences in choice of conventions) recovers the results in [78] (see also [76]) for the Abelian one-particle-irreducible fermion loop diagrams with n photon vertices attached. In this construction of Abelian loop-diagrams, the integration in proper time of the terms ∼ k 1 ·k 2 G B (τ 1 , τ 2 ) in Eq. (259) containing no Grassmanian structures will lead to the standard denominator in a scalar loop integral in Feynman parameters. This piece carries no information on the spin and color degrees of freedom of the particles involved. The rest of the terms (involving derivatives of G B (τ 1 , τ 2 ) and the fermion Green functions G F (τ 1 , τ 2 ), accompanied always by Grassmanian variables) will contain the numerators of the loop diagrams with all the relevant information about the spin dependence of the particles in the loop. It is worth noting that in a conventional construction of these diagrams in Feynman perturbation theory, one has to go through lengthy algebraic manipulations (of matrix structures and the loop-momentum integrals) to extract spin-dependent pieces of the diagrams.
The prior discussion can be extended to non-Abelian gauge theories. In this case, the gluon vertex involves the extra quadratic term in the gauge field in F µν , which is introduced and evaluated in Eq. (259) through a pinching procedure [76], whereby different gluons are pinched to the same point τ i = τ j in the fermion loop; the extra kinematic factor in the spin-dependent numerator in Eq. (259) can be extracted in a systematic way leading to the well-known Bern-Kosower rules [77,78] to construct one-particle-irreducible fermion loop diagrams in QCD.

D Generalized Bargman-Michel-Telegdi equations and QED worldline instantons
The integration of the matter and gauge fields out of the action formulates QED in terms of many-body first-quantized spinning charges described by even x(t) and odd ψ(t) 0+1dimensional worldlines which interact through the Lorentz forces, plus a set of einbein functions containing the (super) gauge symmetry constraints of the worldline action under propertime reparametrizations. Quantum dynamics is implemented via the sum of over all possible worldline paths and einbein configurations. The classical equations of motion of the spinning particle [69,71,72,135,159] are recovered as a special case of the full quantum scenario [125,126]. We briefly review here the strategies to obtain so-called "worldline instanton" configurations from variational principles [85,86,92,113,133,160] and the role of einbeins in setting the energy-momentum and helicity-momentum constraints. The steepest descent method applied to the amplitudes we found yield a slight more complicated picture than the case of a spinning particle moving in an arbitrary background. Gauge fields are fully dynamical now; hence, the instanton configuration of each particle back reacts on to the instanton configuration of the others, entangled by the particular nature of the spin-spin couplings.
We consider first the instantons of the fermion determinant containing the loop electrons with χ = 0. The einbein saddle point is trivial; from ∂ π H = 0, we findε = 0. Thus ε = ε 0 is a constant of motion. It is constant as well for all worldlines contributing to the amplitude, since the gauge-fixing term, trivially integrated in π, producesε = 0 for all configurations.
we get the energy-momentum constraint. Similarly, allows one to interpret Eq. (263) as fixing the einbein: In the last step above, we changed to arbitrary parametrizations τ (τ ). If τ is chosen as the propertime s thenẋ 2 (s) = 1; substituting Eq. (265) into Eq. (262), we recover the gauge-fixed action in [72], The gauge-fixed Lagrangian as a function ofẋ µ can be obtained now by substituting Eq. (264) in Eq. (262), and substituting the expression for the fixed ε in Eq. (265). Since H cl vanishes when this is imposed, we get which is the classical Lagrangian for the spinning charge in a closed loop [72]. Helicitymomentum constraints are not present since the fermion determinant has been squared without loss of generality. The gauge-fixed equations of motion are where we defined σ µν = i[ψ µ , ψ ν ]/2 and used again an arbitrary τ (τ ). For the propertime, this simplifies to and the spin dynamics of the worldlines is described by For homogeneous backgrounds, Eqs. (269) and (270) reduce to the Bargman-Michel-Telegdi (BMT) equations [103] in covariant form, with wide applications, an example being their fundamental importance in describing the properties of highly relativistic beams of leptons and hadrons in accelerator physics [161][162][163][164]. For arbitrary backgrounds, Eqs. (269) and (270) contain corrections to the BMT equations. The generalization to colored fields are the Wong equations [104,105], and in gravity, they are the Papapetrou-Mathisson-Dixon equations [106][107][108].
This worldline instanton configuration can be used to evaluate the action. Since H vanishes for these configurations, the integration in ε 0 introduces a trivial (infinite) normalization to the fermion determinant times a pure phase factor, with the latter given by the coupling of the scalar charged current to the backgroundẋ µ A µ , as expected.
We consider now the instanton configurations of real electrons (χ = 0), Hence in this case, Inserting back Eqs. (272) and (271) into Eq. (262), using ∂ ν H =π = 0 and noting that the energy-momentum relation cancels the even ε contribution of the Hamiltonian, we obtain which is the Lagrangian of a classical spinning charge with helicity constraints [72]. Further expanding the odd contributions in the Lagrangian in χ, and noting χ 2 = 0, we get The saddle-point configuration for the odd einbein produces the constraint with s(τ ) propertime, so that the instanton configurations are those of the Lagrangian in Eq. (267) supplemented now with the helicity constraint. The helicity-momentum constraint, is as expected, modified in presence of background fields with magnetic components for which m R = m.
We address now the most general case in which the background field is sourced by dynamical worldlines. The worldline action in the -th loop order contribution to the functionally averaged n-point function is, The first term is the free action of the n external particles, the second the free action of the loop particles and the last term all their Lorentz interactions. The Lagrangian is non-local hence the saddle points will be solutions of integral equations. Expanding the action S = S cl + δS cl + δ 2 S cl + · · · (277) we get to first order where the field strength tensor is given recursively in terms of the scalar and spinorial charged currents and the einbeins F µν all (x) := ∂ µ A ν all (x) − ∂ ν A µ all (x) with The condition δS = 0 yields the classical equations of motion satisfied by n + spinning charges in interaction (the ones present in the -th order contribution to the functionally averaged n-point function) plus the energy-momentum and helicity-momentum constraints. However, unlike the case of an external field, in which the interactions are given for all particles at once, the motion described here addresses a much more involved scenario. Since the worldlines themselves -including the classical self-interaction term β = α -do participate in creating the field at all retarded and advanced times, the classical equations for the α-th particle worldlines x α (τ ) and ψ α (τ ) are given by integral equations.
Using that 2πδ(q) = lim then J(β p , β q ) = i 8π 2 (β p ·β q ) The first of these two Dirac deltas above does not contribute, since the argument is always positive. The second can be rewritten as where y 0 is the solution of |β p − y 0 β q | = (1 − y 0 ). Then finally where the relative cusp angle is the one in Eq. (154) of Section 5.2 again.