Collisional strong-field QED kinetic equations from first principles

Starting from nonequilibrium quantum field theory on a closed time path, we derive kinetic equations for the strong-field regime of quantum electrodynamics (QED) using a systematic expansion in the gauge coupling $e$. The strong field regime is characterized by a large photon field of order $\mathcal{O}(1/e)$, which is relevant for the description of, e.g., intense laser fields, the initial stages of off-central heavy ion collisions, and condensed matter systems with net fermion number. The strong field enters the dynamical equations via both quantum Vlasov and collision terms, which we derive to order $\mathcal{O}(e^2)$. The kinetic equations feature generalized scattering amplitudes that have their own equation of motion in terms of the fermion spectral function. The description includes single photon emission, electron-positron pair photoproduction, vacuum (Schwinger) pair production, their inverse processes, medium effects and contributions from the field, which are not restricted to the so-called locally-constant crossed field approximation. This extends known kinetic equations commonly used in strong-field QED of intense laser fields. In particular, we derive an expression for the asymptotic fermion pair number that includes leading-order collisions and remains valid for strongly inhomogeneous fields. For the purpose of analytically highlighting limiting cases, we also consider plane-wave fields for which it is shown how to recover Furry-picture scattering amplitudes by further assuming negligible occupations. Known on-shell descriptions are recovered in the case of simply peaked ultrarelativistic fermion occupations. Collisional strong-field equations are necessary to describe the dynamics to thermal equilibrium starting from strong-field initial conditions.

Present and upcoming laser facilities [1][2][3][4] promise unprecedented insights into the strong-field regime of quantum electrodynamics (QED). Strong dynamical electromagnetic fields are also generated during the early stages in off-central collisions of heavy nuclei at the Large Hadron Collider (LHC) at CERN or the Relativistic Heavy Ion Collider (RHIC) at BNL. The presence of strong electromagnetic fields and their dynamical decay can lead to a wealth of intriguing quantum phenomena, such as related to quantum anomalies which can be probed also in condensed matter systems [5]. Strong fields are also essential for the description of highly charged systems, where the net fermion charge induces strong field configurations also in equilibrium [6]. While experiments pioneered by the Stanford Linear Accelerator Center (SLAC) [7][8][9] have since been developed further [10,11], they are not yet able to enter the full strong-field QED regime by means of lasers. Meanwhile, experiments employing crystals have been found to be a competitor to laser experiments [12][13][14][15].
For the weak QED coupling α = e 2 /4π ≈ 1/137 (we use natural units with = c = k B = ε 0 = 1), the strong-field regime may be characterized by a photon field that is parametrically as large as For a laser field [16] that is described by an electric field amplitude E and frequency ω, the counting rule (1) corresponds to a large (Lorentz-invariant) non-linearity parameter [15][16][17], |e|E/(mω) 1 .
For a macroscopic photon field that varies on the time scale of the Compton length 1/m, the counting rule (1) corresponds to electric fields of the order of the critical field, which induces electron-positron pair creation from the vacuum [18][19][20][21][22]. Despite the smallness of the QED coupling, the theoretical description of strong field phenomena provides important challenges. Standard simulation techniques, such as based on Monte Carlo importance sampling, cannot be applied to general nonequilibrium problems. Rigorous simulations are difficult even in equilibrium in the presence of a net fermion charge leading to non-vanishing fields. As a consequence, the development of suitable approximate treatments is indispensable.
For instance, the decay times of strong electromagnetic fields in the medium created by a heavy ion collision and the role of the fields for the subsequent nonequilibrium dynamics is still poorly understood. Even the idealized problem of how an initially supercritical homogeneous electromagnetic field approaches thermal equilibrium in QED has not been answered yet. The strong field regime at early times may be accurately described by classical statistical field theory techniques [23,24], while the late time behavior at high temperature in the absence of a field is successfully described using standard kinetic theory [25]. In particular, the dynamics of avalanches in which large amounts of fermions are produced can be captured by a kinetic approximation of QED [22,[26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43]. However, to describe in a single approach the evolution all the way from strong fields to equilibrium, or in the presence of a net fermion density, involves the interplay of strong fields and collisions beyond state-of-the-art approximations [44].
As an important step in this direction, we derive in this work dynamical equations for strong fields in a kinetic description including collisional processes to order O(e 2 ). Our ab initio derivation starts from nonequilibrium quantum field theory on a closed time path [45,46]. We derive coupled equations for the spatio-temporal evolution of the field expectation value and correlation functions for commutators and anti-commutators of fields using two-particle irreducible (2PI) generating functional techniques [47,48]. The expectation values of field commutators (anti-commutators) for bosons (fermions) describe the spectral functions of excitations, whereas their anticommutators (commutators) characterize their transport behavior.
Applying a gradient expansion for two-point functions, we derive a kinetic description where the strong-field scattering kernel couples the transport equations for photons and fermions to an equation for the fermion spectral function. The latter includes strong-field off-shell corrections in a self-consistent way. Our description incorporates the processes of single photon emission, electron-positron pair photoproduction, vacuum pair production, their inverse processes, medium effects and contributions from the field going beyond the so-called locally-constant crossed field approximation (LCFA) [16]. In fact, we show that our approximation order captures already the complete explicit field-dependence of the problem. To make further contact with the literature, we also consider plane-wave fields. Plane-wave degrees of freedom are identified and it is shown how to recover Furry-picture scattering amplitudes.
The structure of this paper is the following. We introduce the nonequilibrium equations of motion for one-and two-point correlation functions in Sec. II. The ingredients for a kinetic limit of these equations are discussed in Sec. III. We establish the systematics of counting couplings and gradients in the presence of a strong field and present general strong-field transport equations in Sec. IV. In Sec. V, we point out which additional physical assumptions are necessary to reduce the collision kernels of our transport equations to various known expressions and kinetic equations in the literature and how to describe strong-field pair production in our formalism. We conclude and give an outlook in Sec. VI.

II. NONEQUILIBRIUM QED
All possible information about the dynamics of quantum fields is contained in their correlation functions. The latter can be efficiently encoded in terms of a quantum effective action, which is the generating functional for time ordered field correlation functions. Here we employ the two-particle irreducible effective action Γ[A, D, ∆], which is a functional of the macroscopic field expectation value A µ (x) = Tr {ρ(t 0 )A µ (x)} =: A µ (x) (4) with Heisenberg gauge field operator A µ (x) for given density operator ρ(t 0 ) at initial time t 0 , as well as of the time-ordered connected two-point correlation functions ∆(x, y) = T C Ψ(x)Ψ(y) , for gauge fields and Dirac fermions with fermion field operators Ψ andΨ := Ψ † γ 0 , where we suppress spinor indices. The expectation value of the fermion field Ψ vanishes identically for the dynamics considered and plays no role in the following. The symbol T C denotes contour time ordering on the closed time path C [46], which starts at initial time t 0 and runs along the time axis and back as indicated in Fig. 1. Together with a non-thermal, ρ(t 0 ) = e −βH , and not timetranslation-invariant, [ρ(t 0 ), H] = 0, density matrix the contour can be used to facilitate a compact formulation of quantum field theory as an initial value problem that describes non-equilibrium physics. It is convenient to write the 2PI effective action as [72][73][74][75] where Tr C G := x,C G(x, x). This identifies the pure gauge field part of the gauge-fixed classical QED action with the gauge-invariant field strength tensor and gauge-fixing parameter ξ. We use Lorenz gauge, and keep in mind the possibility for residual gauge-fixing.
If computed within the 2PI loop expansion introduced below without a further kinetic limit, correlation functions such as (4) depend on the gauge-fixing parameter ξ (see also Sec. IV E). This gauge-fixing dependence occurs at a higher perturbative order in the coupling than the actual approximation order [76][77][78] and can be absent in the limit of on-shell photons relevant for kinetic descriptions [79] (see also Eq. (101)). In the present paper, we discuss this in the context of Ward identities in the presence of strong fields in Sec. V E 2, where we show that the gauge-fixing parameter drops out in limiting cases.
The semi-classical or 'one-loop' terms in (7) contain the classical photon and fermion propagators in the presence of the macroscopic gauge field with / A := γ µ A µ etc. Our metric convention is η µν = diag(+1, −1, −1, −1). The benefit of the decomposition identity (7) for the full quantum effective action Γ[A, D, ∆] is that the remaining functional Γ 2 [D, ∆] exhibits specific properties that are very useful for the following. For QED, Γ 2 is the sum of all 2PI contributions built from the full two-point functions D and ∆ and there is no explicit dependence on the macroscopic field A, which is further discussed below. A diagram is 2PI if it cannot be disconnected by cutting two propagator lines (see Fig. 2). The 2PI functional integral approach provides a prescription on how to close equations in terms of one-and twopoint correlation functions only. Such a correlation function based description may be used to initialize the system for instance with vanishing photon and fermion particle FIG. 2. Examples of 2PI and two-particle reducible diagrams. number, described by connected two-point functions, but large electromagnetic field or vice versa.
Furthermore, the 2PI formulation is known to facilitate a derivation of kinetic equations [80,81] and may be transformed into other common formulations: Wigner transformations of 2PI two-point functions allow one to make contact with the Wigner operator formalism [50,51,[82][83][84][85]. In particular, equal-time Wigner functions emerge from integration over frequencies [51]. In this way one is also able to make contact with the equal-time Dirac-Heisenberg-Wigner (DHW) formalism [20,86] which has been applied to the discussion of pair production from collisionless equations. Such quantum Vlasov equations [20,53,68,[87][88][89][90] emerge under the so-called 'meanfield' (or 'Hartree-Fock') approximation, Γ 2 ≈ 0. In an operator formulation, this approximation allows one to close operator equations by treating photon operators classically, at the cost of losing access to collisions. In the 2PI formulation, one can easily go beyond this meanfield order e.g. by means of the 2PI loop expansion as is discussed below. This way of arriving at a kinetic description starting from an effective action formulation has the additional advantage that observables derived from that effective action also become accessible under the kinetic approximation.

A. Equations of motion
The equations of motion for the full one-and two-point functions A µ (x), D µν (x, y), ∆(x, y) appearing in the 2PI effective action (7) are obtained from the stationarity conditions 1 These are coupled partial integro-differential equations for the one-and two-point functions on the closed time contour. From them emerge a Maxwell equation, and photon and electron-positron transport equations respectively. 1 These equations are valid in the absence of external source terms. Sources encoding initial conditions are stated accordingly together with the differential equations for the fields and propagators.
In order to discuss the equations of motion, it is convenient to make the time ordering explicit by writing After splitting the equations of motion into equations for the 'statistical functions' (F ) and 'spectral functions' (ρ), the contour C no longer appears and a clear separation into transport and spectral dynamics is achieved. These functions have distinct hermiticity properties, These properties are related to the underlying (anti-)commutator representations in terms of Heisenberg field operators: In particular, the equal-time (anti-)commutation rules are encoded in the spectral functions according to These equal-time conditions imply that spectral functions are normalized and that their initial conditions are fixed by the underlying quantum theory. An important simplification in Abelian theories such as QED occurs because of the absence of 2PI one-point function diagrams, such that Γ 2 [D, ∆] does not explicitly depend on A: The electromagnetic field expectation value enters the 2PI effective action for QED via the 'classical vertex' term which can be depicted graphically as in Fig. 3. Such a contribution cannot be found in the 2PI diagrams contributing to Γ 2 since the two fermion lines emanating from the vertex could always be cut, thus making any such diagram two-particle reducible (see also Fig. 2). For QED, the macroscopic field therefore enters the 2PI effective action (7) only via the classical fermion propagator ∆ 0 [A] and the classical action S[A], while Γ 2 is field-independent. Since 2PI diagrams are at least twoloop, this implies that the complete explicit macroscopic field dependence enters at one-loop order of Γ, Consequently, the field evolution equation always has a Maxwell-type form, i.e.
with the fermion current (see appendix B) irrespective of the approximation order for Γ 2 . This would not be the case, e.g., in QCD or self-interactring Φ 4 scalar field theory, where the two-particle irreducible part of the effective action depends explicitly on the field expectation value, such that the form of the field evolution equation depends strongly on the order of approximation. Because Γ 2 is field-independent in QED, there are no further terms coming from higher order corrections. Approximations to Γ 2 affect the field evolution only implicitly via F Ψ in the fermion current (30). Furthermore, that each 2PI diagram in Γ 2 is separately gauge-invariant in QED [91] remains true in the presence of a macroscopic field due to the field-independence of Γ 2 . Notably, a vanishing field is not in general a selfconsistent solution: if the system is initialized with a finite net charge density, it will develop a field from fermion fluctuations in the Maxwell equation. This field is then necessarily inhomogeneous as dictated by Gauss's law, i.e. the 0-component of the Maxwell equation. Therefore, if the system equilibrates, it has to do so under this constraint for inhomogeneity.
In the equations of motion for the two-point functions, explicitly field-independent self-energies are given by and can be decomposed similarly to two-point functions: With these definitions, assuming Gaussian initial conditions, the stationarity conditions for the propagators in Eq. (13) can be written as 2 [93] While the structure of these equations is determined by causality, details of the underlying theory enter through the differential operators and self-energies, which couple all spectral and statistical functions to each other.
The fact that initial conditions for spectral functions are fixed by the equal-time (anti-)commutation relations (24) - (26), is reflected by the absence of the initial time t 0 in the memory integrals of their equations. In contrast, the evolution equations for the statistical functions have to be supplied by initial conditions. Non-Gaussian quantum fluctuations are built up dynamically but vanish at initial time, x 0 = y 0 = t 0 , by vanishing of the memory integrals.
All equations are considered to be suitably regularized and the renormalization of the 2PI effective action for QED is discussed in detail in Ref. [94]. Since we will finally arrive at a set of finite equations at the level of the kinetic approximation, renormalization will not be further discussed and we refer, e.g., to Refs. [52,53] for details concerning dynamics.
The self-energies, encoding collisions, have leading contributions at Σ, Σ Ψ ∼ O(e 2 ). While self-energies have no explicit dependence on the macroscopic field by their definition in terms of the field-independent Γ 2 , fermion two-point functions introduce an implicit field-dependence when evaluated from their equations of motion. As we will demonstrate, strong-field collision kernels are generated both in photon and fermion kinetic equations in this way. The macroscopic field enters via the terms e / A(x), encoding in particular the Vlasov terms of fermion transport equations, which can be any order depending on the strength of A µ (x). By the smallness of the coupling e, these terms are suppressed in a naive power counting. However, in the presence of a strong field, A µ (x) ∼ O(1/e), these terms are effectively of order eA µ (x) ∼ O(e 0 ) such that the field-vertex (27) has to be resummed. As the macroscopic field decays [31,53,95] from its strong-field initial conditions, the system passes through different power counting scenarios that are all captured by our strong-field counting.

B. 2PI loop expansion
In order to close the equations (35)-(38) one requires explicit expressions for the self-energies (31) and (32). This is achieved by employing a 2PI coupling or 'loop' expansion, which expresses Γ 2 in terms of resummed propagators D µν and ∆ and of free vertices. This selfconsistent treatment of propagators selectively resums perturbative contributions, which helps achieving a nonsecular time evolution with a valid expansion scheme at all times [93,96]. In such an expansion, Γ 2 can be written as where we have suppressed all indices and arguments that are contracted or integrated over. This is diagramatically depicted in Fig. 4. The explicit expressions obey Feynman rules including symmetry factors. Only the free QED vertex appears. Correspondingly, the 2PI loop expansion of the selfenergies (31) and (32) is a series of 1PI diagrams with two amputated external legs (see Fig. 5). The 1PI property of the self-energies can also directly be understood from the definition of Γ 2 as the sum of all closed 2PI diagrams, from which Σ, Σ Ψ are obtained by opening one propagator line, i.e. by Eq. (31), (32).
As long as photon occupations are not too large, i.e. if the statistical photon two-point function obeys the power counting of e from vertices in a 2PI loop expansion can be expected to be a valid approximation scheme and we can truncate by virtue of the smallness of e. Similar conditions for the spectral functions always hold since they are normalized by the equal-time commutation relations. Since fermion occupancies are limited by Fermi-Dirac statistics there are no further corresponding constraints for the expansion scheme. The condition (41) is dynamical such that even if the system is initialized with small occupations, a kinetic description breaks down if too many photons with the same position and momentum are produced. Physically, the assumption (41) may be understood as the requirement for a sufficiently long mean free path in kinetic descriptions: The loop expansion of self-energies in the kinetic limit is an expansion in the number of particles involved in a scattering [97][98][99]. The denser the medium, the smaller the mean free path, and the more likely a collision involving many particles. If the medium is too dense, collisions between arbitrarily many particles become equally likely, invalidating a truncation in an expansion of the number of particles. 3 We emphasize that these considerations do not directly limit the size of the macroscopic field: Because of the field-independence of Γ 2 , higher order contributions to self-energies are negligible also in the presence of strong fields and processes such as eeγγ or eeee scattering do not contribute to a leading-order (LO) description (see also Ref. [101]). As long as (41) is fulfilled, the coupling remains a valid expansion parameter, no matter how large the field is at that time. Thus we may employ the leading order of self-energies to obtain a closed description that is complete at order O(e 2 ).
The corresponding self-energy expressions are where the relative sign originates from the fermion loop in Σ µν . The kinetic equations derived in this paper neglect all higher orders of the 2PI loop expansion. 4 In agreement with the coupling counting in perturbation theory, all possible crossings of eeγ scattering terms emerge from these O(e 2 ) self-energies. The following sections are dedicated to understanding how effective transport and kinetic descriptions emerge from this approach.

III. THE KINETIC LIMIT OF NONEQUILIBRIUM QED
To express the equations of motion in kinetic degrees of freedom, we change to center and relative space-time variables The four-momentum p associated to −i∂ s is the momentum that appears in kinetic equations, while X is the kinetic four-position variable. The momentum p is introduced by a Wigner transform with respect to the relative coordinate s. For an evolution starting at time t 0 at which the initial conditions are given, the Wigner transform of a generic two-point function G may be written as Here t 0 appears in the time integral as a lower boundary for all time variables. Since initially we have X 0 = t 0 , there are no relative times to integrate in this case, which preempts a Wigner transformation starting at initial time.
To nevertheless be able to talk about kinetic variables from the initial time of our kinetic description, we employ a late-time limit described in the following.

A. Late-time limit
For finite t 0 and X 0 the integration range for s 0 is always limited. Only if t 0 → −∞ the relative time variable s 0 can be infinite, which is required for a proper introduction of Fourier frequency modes p 0 . Of course, sending formally t 0 → −∞ while still initializing the evolution at some finite time implies that a general system is initially not accurately described by these late-time equations. However, for sufficiently large X 0 compared to the finite initialization time, the description is expected to become accurate [102]. Therefore, instead of Eq. (46) we consider the late-time Wigner transform which has contributions from all s 0 for arbitrary X 0 . Equal-point objects such as the fermion current (30) can be expressed in terms of such late-time Wigner transforms, The notation p = ∞ −∞ d 4 p/(2π) 4 for momentum integrals is used throughout. The canonical equal-time anticommutator (26) such that the late-time vector-zero component 1 4 tr{γ 0 ρ Ψ (X, p)} may be interpreted as a density of states [103].
In the microscopic description, finite-time Wigner transforms (46) produce factors with finite-width energy-peaks on correspondingly small timescales [104] that reduce to delta peaks at late times via In this late-time regime, the interactions of QED may be described by those of kinetic theory in terms of degrees of freedom that carry a definite amount of energy. Applying the late-time limit, t 0 → −∞, one can write the equations of motion (35) - (38) as with z = d 4 z, where we have introduced the retarded and advanced functions for photons and fermions (A5) -(A7) defined in appendix A. Given the multitude of different nonequilibrium twopoint functions, it is important to remember that there are only two independent two-point functions per field species: the statistical and spectral functions. However, this can be invalidated by approximations, in particular, by the procedure of sending t 0 → −∞ while initializing the equations at a finite time. Wigner functions that include small frequencies via (47) may appear independent of each other because of spurious small frequency contributions that, in an exact description employing finite-time Wigner transforms (46), do not yet exist at early times [105,106].

B. Gradient expansion
As a next step in the derivation of kinetic equations, one considers an expansion in the Lorentz-invariant and dimensionless parameter (s · ∂ X ). An expansion in propagatorgradients is achieved by the late-time identity [102] , which applies to photon and fermion convolutions Expansion of the exponential in Eq. (55) corresponds to an expansion in (∂ p · ∂ X ), i.e. a gradient expansion in Wigner space. While the LO simply replaces the Wigner transform of convolutions by products of Wigner transforms, an expansion to next-to-leading order (NLO) in propagator-gradients would involve Poisson brackets, The truncated gradient expansion leads to equations that are irreversible and local in central time X 0 , as in the case of kinetic equations. Still, gradient expanded 2PI equations contain parts of the memory integrals of the fundamental equations and are non-local in relative time s 0 . This allows for access to unconstrained frequency variables, which are not present in traditional kinetic descriptions as further discussed in the following sections.
The smallness of the expansion parameter (s·∂ X ) can be met in several circumstances. 5 Quantum field dynamics often becomes insensitive to its past, such that correlations are dominated by small s [107][108][109][110]. From the perspective of the spectral function, this damping of correlations in time corresponds to the emergence of a particle picture in momentum space [111,112]. Furthermore, assuming that (s · ∂ X ) is small depends on what the derivative acts on: In the following, we neglect only gradients of two-point functions G, by dropping Poisson brackets while formally keeping gradients of the gauge-invariant field strength tensor, (s · ∂ X ) j eF µν (X) ∼ O(e 0 (s · ∂ X ) j ), to all orders. That is, we count field-gradients as (s · ∂ X )F µν ∼ F µν and propagator-gradients as (s · ∂ X )G G. This allows us to treat a large class of far-from-equilibrium initial conditions of the macroscopic field. Approximations to field-gradients are then discussed in Secs. V A and V F 2, where we make contact with the locally-constant field approximation.
However, field-gradients may be implicit in propagator solutions (see also Sec. V F 2 for the example of planewave fields): Given an explicit field-dependent solution for a two-point function, for example of the form different gradients may be related via In fact, the separation of field and propagator-gradients that is possible at the level of the equations of motion does not ensure that the ratio (61) is small. Nevertheless, we can observe from Eq. (61) that large fermion momenta can facilitate such a separation. When solving the kinetic equations for inhomogeneous fields derived below, the smallness of the ratio (61) should be be checked.
In the presence of chiral symmetry (facilitated by massless fermions or ultrarelativistic momenta), scalar, pseudoscalar and tensor components vanish identically [108]. If a description in terms of free particles is valid, the axial component of the free fermion spectral function would also vanish.
Similar comments apply to the photon distribution function and a decomposition of the Lorentz tensor structures of the photon equations of motion in the presence of a macroscopic field can be achieved with the basis discussed in Refs. [113,114].
With this in mind, one could write without loss of generality for each component of F Ψ (X, p): The change from a description in terms of F Ψ,S...T (X, p) to a formulation in terms of f Ψ,S...T (X, p) is convenient because in characteristic limits f Ψ,S...T (X, p) can be interpreted as distribution functions.
In particular, in thermal equilibrium all distribution functions are time-independent and equal the Fermi-Dirac distribution, i.e. f Ψ,S...T (p 0 ) = 1/(e βp 0 + 1) (and correspondingly a Bose-Einstein distribution for the photon case). For a thermal theory this is valid no matter how strong the interactions are and holds even in the absence of a dispersion relation between frequency and spatial momenta, p 0 = ω( p ).
Phenomena such as the chiral magnetic effect [115][116][117][118], chiral kinetic theory [119][120][121][122][123] or spin transport [124] should become accessible from first principles by using (68) - (72) in the equations of motion (51) - (54). However, for our current purposes of strong-field kinetic equations and to make contact with existing limiting cases in the literature, we consider a single distribution function f Ψ (X, p) for fermions and f (X, k) for photons by writing [93,125] For the fermion distribution function one has Pauli's principle [108], In order to distinguish fermion and anti-fermion distribution functions, it is convenient to define [126] f Ψ (X, p) =: In a charge conjugation invariant system, the fermion distribution function obeys [73,74] such that the system is charge neutral, While the vacuum is CP-invariant, the general initial conditions which we want to discuss in this paper break CP-invariance by introducing a net total charge, such The photon identity analogous to (77) reads [73,74] −[f (X, −k) and does not rely on CP-invariance.

On-shell particle picture
In general, the distribution functions introduced in Eqs. (73) and (74) depend on the off-shell frequency variable p 0 that is not restricted to any dispersion relation, p 0 = ω( p ). However, they only appear in combination with the respective spectral function. As a consequence, if the physics can be approximately described by free spectral functions, then the distribution functions can be restricted to their on-shell values. Whether an on-shell description is possible is determined self-consistently by solving the equations of motion (52) and (54) for the spectral functions. At initial time, the photon (fermion) spectral functions are determined by the equal-time (anti)commutation rules and each subsequent time step is determined by the equations of motion. If and when on-shell spectral functions emerge depends on timescales and initial conditions for statistical propagators and the macroscopic field. As we argue in Sec. IV B, the free fermion spectral function (81) is in fact not complete at order O(e 2 ) in the presence of general strong fields, A µ (x) ∼ 1/e, such that a standard on-shell kinetic description breaks down. Instead, we propose in this paper a less restrictive 'transport' description that includes off-shell frequencies of fermions (but not of photons) in terms of the off-shell distribution function f Ψ (X, p). The frequency dependence of this function is then determined dynamically by the equations of motion and independently of its momentum dependence p. An electron and positron particle picture is assumed only in Sec. V B to compute particle production at asymptotic times when the field has decayed. With this application in mind, it is instructive to compute the fermion current (48) for the free fermion spectral function (81), i.e.
with on-shell electron and positron distribution functions, The zero-component (82) can be interpreted in terms of the conserved electric charge which then reads on-shell Similarly, on-shell, j i gives rise to the fermion pair number density which is related to the total pair number via This expression will serve us to define an asymptotic particle number of strong-field systems in Sec. V B.
In contrast to the fermion case, the photon spectral function may be set to its free form also in the presence of a strong field (see Sec. IV A). We can then identify the on-shell photon distribution functions of kinetic theory by integrating over frequency k 0 , i.e.
as we discuss below. The total number of photons is then

IV. STRONG-FIELD QED TRANSPORT EQUATIONS
We now apply the procedure of Sec. III to the equations of motion (51) -(54) for the statistical and spectral functions and the equation of motion (29) for the macroscopic field. To ease the notation, we refer to the left sides of the two-point function equations as (F LHS) µν (x, y), (ρLHS) µν (x, y) and (F LHS) Ψ (x, y), (ρLHS) Ψ (x, y) respectively, and similarly to the right hand sides ('RHS') or to entire equations ('EOM').
To reveal the 'gain-minus-loss' structure of collision terms, we identify the '+ − / − +' or 'Wightman functions' (defined in appendix A) by making use of the identity and an analogous identity for fermions. Then Eqs. (73) and (74) can be expressed in terms of the Wightman functions as From the '+−' functions, one readily observes the appearance of Bose-enhancement terms (1+f (X, k)) for photons and Pauli-blocking terms (1 − f Ψ (X, p)) for fermions. In collision terms, these emerge attached to outgoing particles, while ingoing photons and fermions, associated with '−+' functions, are not distinguished in terms of their statistics.

A. Photon spectral function & gauge-fixing independent photon drift term
The photon transport equation is related to the evolution equation of the statistical photon propagator via i.e. by a Wigner transformation of Lorentz-traced differences. Combined with the change of variables to X and p, the Boltzmann derivative operator is recovered from the d'Alembertian in a Lorentz-invariant way by the identity By use of the convolution identity (55) at LO gradient expansion as well as of symmetry properties of the Wigner transforms given in appendix A, one finds that Eq. (97) reads The tracing over Lorentz indices reduces the ten equations for the components of F µν to a single scalar equation.
In combination with the introduction of the distribution functions (73) and (74), which reduces the amount of independent tensor structures, (99) is then sufficient to close the dynamics. This transport equation (99) is valid to all orders in the coupling of the 2PI loop expansion. To obtain a leading order collision term, we neglect terms of order O(e 4 ) to this equation. There are two types of such higher order terms: a) terms of order O(e 4 ) in Γ 2 discussed in Sec. II B; b) terms of order O(e 2 ) in equations of motion for spectral functions contributing to the transport equations only at order O(e 4 ). Terms of the latter type appear in the analogous expression (97) for the photon spectral function, i.e.
The O(e 2 ) terms of this equation contribute only at order O(e 4 ) to the transport system, because the self-energies in Eq. (99) are already of order O(e 2 ) themselves, before being multiplied with the photon two-point function containing ρ µν . It is therefore sufficient at order O(e 2 ) to employ the free O(e 0 ) solution of Eq. (100) in transport equations. In this way, transport equations self-consistently resum statistical functions, but not spectral functions. The additional 'collisional broadening' of spectral peaks, that does not enter the LO strong-field transport description explicitly, can then be estimated from its solutions, e.g. by evaluating spectral self-energies in terms of distribution functions or computing the decay rate. The gradient expansion further supports this special treatment of spectral functions, as we discuss in appendix C.
Employing the free photon spectral function (81) by this reasoning, the gauge-fixing dependence of the LHS of the photon transport equation (99) drops out due to a cancellation between (D µν 0,ξ ) −1 and ρ µν 0,ξ : To obtain Boltzmann-type equations, one finally integrates over frequencies k 0 , leading to the appearance of the on-shell distribution functions f (X, k) defined in (90), where we have made use of (101). This integration explicitly reduces the information that is redundant because an on-shell dispersion relation is valid for photons.

B. Fermion spectral function
Similarly to the photon case discussed around Eq. (100), terms that are of order O(e 2 ) in contribute only at order O(e 4 ) to the transport RHS that is already of order O(e 2 ) itself. Crucially, the fielddependent term in Eq. (103) is of order O(e 0 ) for strong fields and may thereby not be neglected in the O(e 2 ) transport description. In particular, this implies that a simple fermion particle picture may not exist in general strong-field systems. From a kinetic perspective, this is the essential way in which strong-field systems differ from weak-field systems that may still be described by free fermion spectral functions.
The O(e 0 ) solutionρ Ψ [A] of Eq. (103) has a functional dependence only on the macroscopic field A µ . This is in contrast to the exact spectral solution which would be a functional also of F Ψ , F µν and ρ µν . Nevertheless, because of the field-independence of self-energies, the approximate spectral equation (103) contains the complete explicit field dependence. This includes in particular infinite orders of field-gradients: For instance, the traced LHS reads In spacetime regions where the field vanishes one recovers from Eq. (103) the free particle description (81), i.e. ρ Ψ [A = 0] = ρ Ψ,0 . Eq. (103) may therefore be understood as the strong-field generalization of a fermion particle picture. In particular, since the difference between ρ Ψ,0 and ρ Ψ [A] is of order O(eA), one would be allowed to exchange the two in a leading order description for weak fields in accordance with a near-equilibrium quasiparticle picture [127,128].
Rephrasing the equation of motion for the fermion spectral function into an equation for the retarded propagator, In the zero-field case, this equation implies that the spectral function has a peaked shape with a 'width' given by the square of the spectral O(e 2 ) self-energy [129] −iΣ which is indeed O(e 4 ), as anticipated by our counting of couplings in the equations of motion. In the strong-field case, the off-diagonal momentum structure of the field term in Eq. (105) highlights the absence of a simple peak structure of general solutions of Eq. (103). Eq. (105) further shows that the physical reason for this more complex structure is four-momentum exchange between the retarded fermion propagator and the macroscopic field.
We give an analytical solution of Eq. (103) under the assumption of strong external plane-wave fields in Sec. V D, which allows us to showcase the appearance of exponentials exp(O(eA)), that resum the field-vertex (27) as desired. By employing the solutionρ Ψ [A] of Eq. (103) in transport equations, one recovers O(e 2 ) strong-field scattering amplitudes in limiting cases (see Sec. V E 2). A particle picture emerges only in special cases and can change with time (see Sec. V H).
C. Strong-field photon transport equation

Collision term
To obtain the O(e 2 ) strong-field photon collision term from the expression (99) we need the leading order self-consistent photon self-energies, i.e.
The structure of the strong-field photon transport equation is that of Eq. (97) integrated over positive frequencies, ∞ 0 dk 0 /(2π). Spectral functions are evaluated from their equations of motion with the reasoning discussed in the previous paragraphs, i.e.
where the O(e 2 ) strong-field photon collision term is with the trace and the longitudinal projection of the eeγ-collision kernel This general expression derived from quantum field theory plays the role of a generalized scattering amplitude squared that has its own equation of motion [Eq. (103) or equivalently Eq. (138) below] and is adapted to the properties of the macroscopic field at each instance of time. This goes beyond previous approaches that have so far been restricted by additional assumptions on the macroscopic field. In particular, it provides a prescription of how to implement an inhomogeneous macroscopic field in local transport equations. We achieved this by describing collisions in terms of a dynamical strong-field fermion spectral function, which includes all leading order effects. This approach allows for many links to existing literature as we demonstrate in Sec. V. In particular, the collision kernel (115) may be reduced to scattering amplitudes computable from Feynman rules in strong-field QED (see Sec. V E 2).
The collision term (112) features the factorization of interaction terms into a collision kernel and a gainminus-loss term familiar from traditional kinetic equations. While the photon distribution functions can be reduced to on-shell distributions (90) by virtue of the delta function δ(k 2 ) in P µν , this is not in general possible for the fermion distribution function. The off-shell frequency dependence of the latter is computed dynamically by solving the transport system coupled to the fermion spectral equation (103). This allows the collision kernel, to adjust in time to a self-consistent macroscopic field as the system evolves, while still being local in the kinetic position variable X without relying on locally constant fields.

Strong-field photon decay rate
By linearizing and integrating the photon transport equation over position and external momentum, one may find the field-dependent decay rate γ of a photon with momentum k i and position X i at time t := X 0 , with the photon number (91). Such a linearization may be achieved e.g. under the assumption that the system is close to vacuum (i.e. for small distribution functions, see Sec. V G) or in linear response theory around equilibrium [104,109,130] In equilibrium, gainminus-loss terms vanish by energy conservation, resulting in the photon equilibrium decay rate D. Strong-field fermion transport equation Here, we derive the fermion equations that close the transport system in terms of off-shell fermion and on-shell photon distribution functions.

Gauge-invariant fermion correlation functions
The presence of a macroscopic field complicates the gauge-invariance of approximations such as the gradient expansion. This was not an issue in the case of the photon equations where the field is only implicit viaρ Ψ [A] and the photon self energies are gauge-invariant. In the following, before repeating the analogous steps for the fermion transport equation, we express all fermion equations in terms of the gauge-invariant field strength tensor or equivalently in terms of electric and magnetic fields, This is necessary, in particular, in order to identify a gauge-invariant fermion drift term that contains the gaugeinvariant Lorentz force. One can achieve gauge-invariance (as opposed to covariance) by introducing Wilson lines 6 with Γ indicating the path of integration from y to x. The gauge transformation of a Wilson line exactly compensates the gauge transformation of fermion two-point functions, such that the quantitieŝ are gauge-invariant (but path-dependent). It is well known that straight Wilson lines, W := W Γ=[x,y] , facilitate a derivation of gauge-invariant transport equations [50,56,82,83]. Following this approach, we employ and express everything in terms of gauge-invariant latetime Wigner functionŝ Invariant and covariant Wigner functions are related by with the real differential operator By virtue of this relation is simple for small field-gradients (which we discuss in Secs. V A and V F) in which case it becomes the translation One now has to decide whether to identify fermion distribution functions in terms of F Ψ and ρ Ψ as in (74) or in terms ofF Ψ andρ Ψ , i.e.
In principle, f Ψ andf Ψ are arbitrary definitions which can be translated into each other. In particular for small field-gradients one would havẽ In photon equations, the distinction between co-and invariant fermion functions is redundant. This is because, by virtue of one may replace co-and invariant Wigner functions in the gauge-invariant photon self-energy that features a fermion loop, i.e.
In Wigner space this involves two fermion momentum integrals and a delta function. In particular, the fact that implies that one may replace f Ψ withf Ψ if ρ Ψ is replaced withρ Ψ in the photon collision kernel (115). Similarly, because such that this may also be done for the current (160) in the Maxwell equation (29). In this way, one obtains a closed set of equations in terms of fermion distributions of thef Ψtype to any order of field-gradients. We stress that these replacements do not work in reverse (going fromf Ψ to f Ψ ) for the fermion equations to be discussed below, such that a practicable description in terms of f Ψ -type distributions would have to rely on small field-gradients by relying on Eq. (131).

Gauge-invariant equations of motion: 2PI vs. Wigner operator formalism
Having introduced gauge-invariant correlation functions, we can express the gauge-covariant 2PI fermion equations of motion in a gauge-invariant way. We start with the equation for the fermion spectral function, (137) explicitly at our order of interest, Here we have employed the commuting, real and gaugeinvariant differential operators introduced in Ref. [50], Using anti-hermiticity, Eq. (19), one may verify in particular that solutions of Eq. (138) satisfy The second condition, which is satisfied by any strongfield solution, is much weaker than the on-shell condition in the absence of a field, ( / p − m)ρ Ψ,0 (X, p) = 0. Eq. (138) is proven as in the Wigner operator formalism of Refs. [49][50][51]. While the Wigner operator formalism has not been able to provide closed collision terms, the 2PI formalism is able to achieve this: Instead of discussing equations for the normal-ordered product, :Ψ(x)Ψ(y): , resulting in real and imaginary parts with different differential operators [50,51], we distinguish real and imaginary parts of the time-ordered product T C Ψ(x)Ψ(y) (6), i.e. statistical and spectral functions. Their 2PI equations of motion (37) - (38) do not differ by their differential operators, but by the integral structure of their RHS, which automatically ensures the correct hermiticity properties of their solutions, (17) and (19). Because of the absence of these RHS integrals in the approximated spectral equation (103), the anti-hermiticity (19) of the approximate solution has to be prescribed. In fact at 1-loop, i.e. by neglecting collisions, the equations for F Ψ and ρ Ψ without (anti)hermiticity constraints are equivalent and the equation for the fermion statistical function alone is sufficient to discuss transport phenomena as has been done e.g. in Ref. [20]. Going to order O(e 2 ), the self-energy terms of the 2-loop equations for the spectral functions still do not contribute to the kinetic equations as discussed in Sec. IV B, but the self-energy terms of the statistical equations provide collision terms.

Quantum Vlasov term
In order to obtain a gauge-invariant fermion transport equation, we consider , where (anti-)commutators are taken in Dirac space. By building differences, the fermion mass drops out of this expression, but enters again via the spectral equation (138). By taking the trace of (143) we obtain the all-order in field-gradients quantum Vlasov term to which the commutator term with Π µ does not contribute. In (144) we have indicated the fermion collision term, which we compute to leading order below. Employing Eq. (141), the fermion transport equation (144) in terms off Ψ then reads The off-shell all-gradient drift term of this equation goes beyond a Lorentz force description, which it contains as its on-shell contribution (see Secs. V C and V H). The emergence of this fermion drift term is distinctly different from the photon case, because fermion derivatives involve the macroscopic field and are first order already in the fundamental equations of motion. In particular, the momentum factor of (p · ∂ X ), that emerges automatically for photons via the identity (98), has to be provided by the vector component of the free fermion spectral function.
Without an on-shell approximation, momentum derivatives of the spectral function in Eq. (145) are physically regulated by the macroscopic field.

Collision term & charge conservation
Having discussed the LHS, we now derive the gaugeinvariant collision term already indicated in Eq. (145).
In general, gauge-invariance of the convolutions on the fermion spectral and statistical RHS is achieved by writing where we have identified the (triangle) Wilson loop By virtue of Eqs. (55) and (126), the LO of the gradient expansion of this gauge invariant convolution is [56] For weak fields near equilibrium the additional term as compared to the covariant convolution, , is effectively of order O(e 4 ) and compatible with a kinetic description [56]. To focus on the part of the fermion RHS that contains the collision term indicated in Eq. (144), we drop terms of the type (149) also in the presence of strong fields. We stress that the validity of dropping these terms in a far-from-equilibrium system requires further investigation. 7 7 As discussed in Ref. [56] terms of the form (149) have the effect of accounting for further off-shell corrections and replace the spatial derivative ∂ X → ∂ X − eF µ ν ∂ ν p in Poisson brackets. Alternatively, one may think of dropping these terms as setting the Wilson loop to one, L ≈ 1. Because of the group properties (132), (135) and At leading order, the gauge-invariant self-energies in Eq. (150) may be written aŝ The strong-field O(e 2 ) fermion collision term then reads whereP is obtained from the collision kernel (115) by exchange of ρ Ψ →ρ Ψ with the solutionρ Ψ of Eq. (138), or at LO in field-gradients viã As anticipated in Sec. IV D 1, while the photon collision term is gauge-invariant also without this replacement, the fermion collision term is not. This is because gaugeinvariance requires integration over both fermion momenta according to (134). Indeed, if we integrate the fermion transport equation over its external momentum, subtleties of gauge-invariance are absent and, with and using (136), we can recover the Maxwell current (48) in the fermion transport equation via As a consequence of the U(1) symmetry of QED, this current is conserved by the fundamental equations, as well as by our approximate transport equations, such that the total electric charge (86) is constant, with t := X 0 . To see this, one may verify that the relabeling q ↔ p and k → −k leaves both the delta function and the gain-minus-loss term invariant [by virtue of (79)], but changes the sign of the collision kernel (also without tilde), i.e.

E. Transport Maxwell equation & gauge-fixing dependence
The free photon propagator D µν 0,ξ (11) and spectral function ρ µν 0,ξ (81) introduce a gauge-fixing dependence. This ξ-dependence is distributed over several equations of motion by virtue of P ξ (114) and the solution of the Maxwell equation (29) with the late-time current There are two ways in which ξ-dependence is controlled. Firstly, starting from the 2PI effective action, a perturbative coupling expansion shows that the total ξ-dependence of P ξ [A ξ ] is always of higher perturbative order in e [76][77][78]. Indeed, for a free fermion spectral function ρ Ψ,0 , leading order collisions are trivially gauge-fixing independent, Secondly, the ξ-dependence can drop out for on-shell photons [79] (see also Eq. (101)). We demonstrate this also in the strong-field case by virtue of Ward identities for scattering amplitudes that emerge in the kinetic approximation and play the role of redressed 1PI vertices. We make contact with such strong-field Ward identities [131][132][133] in the case of plane-wave fields in Sec. V E 2, where the ξ-dependence then drops out in a corresponding limit. A general proof for cancellations between (D µν 0,ξ ) −1 in A µ ξ and ρ µν 0,ξ in P ξ in the self-consistent strong-field case P ξ [A ξ ] seems highly non-trivial.
A summary of the interconnections among the extended transport system which we have now arrived at is graphically presented in Fig. 6. The transport equations for photons [Eq. (111) for f (X, k)] and fermions [Eq. (145) forf Ψ (X, p)] couple to eachother via the collision terms (112) and (153) (160). The macroscopic field enters the fermion spectral and transport equation explicitly via the strong-field derivatives (139) and (140), and the photon and fermion transport equations implicitly via the strongfield fermion spectral function in the scattering kernel (115).

V. STRONG-FIELD QED KINETIC EQUATIONS
In this section, we investigate ways to further approximate the transport system of Sec. IV and how to reduce it to Boltzmann-type equations with scattering amplitudes by considering limiting cases of the collision kernel. To this end, we discuss various common additional approximations in strong-field QED, namely small fieldgradients (locally constant fields), classical fermion propagation (Lorentz force), external plane-wave fields (Volkov states), near-vacuum physics (small occupations), as well as fermion distributions that are peaked at large momenta (ultrarelativistic limit). In particular, the ultrarelativistic limit finally allows us to make contact with fermion on-shell descriptions [e.g. Ref. [22]], which are valid if a long-lived separation of scales exists (see Sec. V H).
A. The case of small field-gradients So far, our transport equations have been infinite order in gradients of the macroscopic field. In a physical situation with small field-gradients, one can simplify the collision kernels and the fermion drift term. We demonstrate how to do this at the level of the equations for the fermion spectral and statistical functions in the following. 8 For this purpose we assume in this section that This means we only keep LO terms O e 0 (∂ p · ∂ X ) 0 and truncate the NLO O(e 0 ∂ p · ∂ X ) of gauge-invariant fieldgradients (see appendix D for a comparison of approximations to invariant and covariant field-gradients). We can simplify the fermion spectral equation of motion (138) and in turn the collision kernel (115) by using (162). The exponential derivatives of the differential operators (139) and (140) allow for an expansion in terms 8 A collisionless discussion of field-gradients can be found in Ref. [20], where it is shown that field-gradients can enhance pair production rates in particular for low momenta.
of gradients of the field-strength tensor. Thereby one can explicitly compute the first orders of the λ integrals, i.e. [50,51] ∇ µ (X, p) = ∂ ∂X µ − eF µν (X)∂ ν p (163) Note in particular, that the leading order of ∇ µ , is the classical Vlasov derivative which contains the Lorentz force as its on-shell contribution (see Sec. V C).
Neglecting gradients of the field-strength tensor, the gauge-invariant spectral equation (138) becomes Solutions of equation (167) neglect field-gradients, but are exact in the field strength. This implies in particular that, even for a constant strong field strength tensor, the fermion spectral function is not a delta peak and does not allow for a simple particle picture. 9 The fermion transport equation (145) for small fieldgradients then reduces to 9 This is an essential difference to Yukawa theory [134] or scalar λφ 3 theory (which are diagramatically very similar to QED): the LO equation of motion, e.g. for the scalar spectral function with a strong constant scalar macroscopic field φ 0 ∼ O(1/λ), does have a delta-peaked particle solution Moreover, the equation for the scalar statistical propagator [102] 2(p · ∂ X )F (X, p) has a force term, ∂ X M 2 (X) with M 2 (X) = m 2 + λφ(X), which is NLO of the gradient expansion and vanishes for constant fields.
where we have used the fact that in contrast to ∇ µ , which contains higher order derivatives, D µ satisfies the Leibniz product rule and that a solution of (167) satisfies Plugging the solution of the approximated equation for the spectral function (167) into the collision kernel (115) one obtains photon and fermion collision terms for fields with small gradients. In Sec. V F 2, we demonstrate how the locally-constant field approximation arises from such spectral functions in the special case of plane-wave fields. There, instead of solving the approximated equation (167), we will first solve the infinite order gradient equation (103) (or equivalently (138)) and approximate gradients in the solution in the end.
B. Asymptotic (Schwinger) pair production from unequal-time correlations In this section, as an application of the above small fieldgradient approximation, we discuss how pair production is implemented in the present formalism. We start in the regime of the collisionless Schwinger pair production yield per volume V and time-interval T [18], i.e. the regime of constant fields at 1-loop, and end this section with a general collisional expression for inhomogeneous fields. Under the asymptotic assumption In order to extract the fermion pair number (89), we integrate Eq. (175) over negative and positive energies separately and subtract the resulting integrals (instead of summing them, which would instead give the trivial total charge (157) For the momentum derivatives ∂ i p of D µ we exploit (155), and for its frequency derivative ∂ 0 p we note that 0 −∞ dp 0 − ∞ 0 dp 0 ∂ 0 pFΨ (X, p) = 2 dp 0 δ(p 0 )F Ψ (X, p).
This term eventually acts as a source term to the asymptotic number of fermion pairs. For the position-space derivative we use Finally, the time-derivative in D µ allows us to identify the pair number (89) in the asymptotic past and future, where we have employed the asymptotic assumption (174) and identified the on-shell electron and positron distribution functions (76), (84) and (85) in the asymptotic past and future. Applying the above identities to the 1-loop transport equation (176) gives the result Importantly, this expression relates pair production to self-consistent spectral and field dynamics. The asymptotic assumption (174) only fixes a boundary condition at X 0 → ∞ and interacting spectral dynamics [Eq. (138)] contribute to (180) at all finite times X 0 . In particular, the above expression shows that pair production from the vacuum occurs off-shell at the time of the creation event: the fermion yield (180) vanishes for a free (on-shell) fermion spectral function, because massive fermions can not have zero energy, i.e.
It is only the subsequent evolution, that brings these off-shell contributions from vacuum pairs to the on-shell regime in the asymptotic future, where a particle number N Ψ (∞) is well-defined. Furthermore, the expression (180) vanishes for E = 0, even if B = 0, in accordance with the general statement that magnetic fields can not produce fermion pairs. In our derivation, this is a consequence of the vanishing of momentum derivatives at infinity, i.e. Eq. (155). The structure of the expression (180) is reminiscent of the time-integrated source term of the quantum Vlasov equation from which particle production at zero energy is well known [53]. Such a source term is not manifest in Eqs. (175), (171), or (145), 10 but we have demonstrated here that vacuum pair production is nevertheless contained in these transport equations by coupling to the dynamics of the fermion spectral function.
Since practicable procedures at 1-loop already exist in literature, we want to stress that the significance of Eq. (180) does not stem from its 1-loop practicability but from the fact that it may be systematically generalized and thereby put in the context of thermalization, while other procedures have struggled to do so: At 1-loop, where the equations for spectral and statistical functions are decoupled, one may compute the asymptotic fermion particle number by ignoring spectral dynamics and solving the complete tensorial system for the statistical function. In literature, this is often done in terms of the equal-time 'DHW' functionF Ψ (t, t, x, y )γ 0 , or dp 0F Ψ (X, p)γ 0 in Wigner space. In fact, existing transport derivations of the Schwinger result typically employ such equal-time formulations [20,51,86,137], in which spectral informationn such as a distinction between on and off shell contributions is not explicitly accessible due to spectral functions being constant at equal times [see Eq. (26)]. Equal-time equations can be closed, e.g. by neglecting collisions, but how to close an equal-time description for general strong fields in a controlled approximation is an open problem. From an unequal-time perspective, the equation for the fermion statistical function is not self-sufficient at O(e 2 ), but couples to the fermion spectral function (103), which is not on-shell for strong fields. The unequal-time approach closes by including this equation for the spectral function and is thereby systematically generalizable to higher loop orders that are essential for the approach to equilibrium.
Simply by keeping field-gradients and the collision term, i.e. starting from Eq. (144) instead of Eq. (175), one obtains Due to the presence of higher-order frequency derivatives, the identity (177) is not sufficient to treat inhomogeneous fields, which are able to transfer momentum and produce occupations with finite energy, p 0 = 0. A complete selfconsistent solution of the set of equations in Fig. 6 is necessary to obtain a numerical result for the asymptotic pair number in this way. In particular, the collisional part of (182) contains contributions from 0 → 3 (2-loop vacuum pair production) and 1 → 2 processes ('seeded cascades'), the latter of which dominate over vacuum pair production in subcritical fields [22,27,31]. In contrast to Eq. (182), the 1-loop result (173) describes the effect of a constant external electric field with no feedback from the dynamics of the photon sector.

C. Lorentz force & classical propagation in isolated systems
The Lorentz force, emerges from the quantum Vlasov term of Eq. (145) in the case of a free fermion spectral function and small field-gradients via where the factor of p µ is provided by the vector component Therefore, on-shell particles may be described by the Lorentz force. However, the validity of employing a free spectral function in Eq. (184), i.e. whether on-shell particles indeed dominate the dynamics, depends on the details of the strong-field system: Typical experiments where on-shell particles dominate the dynamics are, for example, those where an electron beam or material target is initially in a zero-field region and then collides with a strong field such as a laser beam [10,11]. In such a setting fermion distribution functions are initialized with occupations only in the on-shell region and the subsequent deviations from on-shell occupations induced by the strong field often remain small even when fermion pairs are produced: This is because these systems feature a separation of time scales due to the typically very large values of the parameter ξ = |e|E/(mω) [16], implying that particles (target or produced) are transported in momentum space to relativistic energies in much less than a laser period. Thereby, the fermion distribution function is typically peaked at an ultrarelaticistic scale and far away from its equilibrium (Fermi-Dirac) shape. At such high energies, off-shell effects can be suppressed [15] and can remain suppressed, if the ultrarelativistic peak in the fermion distribution function is long-lived (see Sec. V H).
In the presence of such long-lived peaked distribution functions, one may then distinguish two kinds of quantum effects [15]: One class is related to the recoil that a fermion experiences during collisions (i.e., emissions of photons). This is controlled by the (spacetime and momentum dependent) parameter [16] which may be small even for large ξ or vice versa. Systems that have small χ may be described completely (both drifting and collisional interactions) in terms of the classical radiation reaction force [16,[138][139][140] that includes collisional corrections to the Lorentz force [22]. The other class of quantum effects is related to how accurate a classical description is between collisions. This is commonly discussed in terms of the de Broglie wavelength /p * , which is then required to be small enough such that the quasiclassical approximation applies [141], and smaller than the mean-free-path such that a separation between propagation and interaction is possible [142]. In our context, p * is the characteristic momentum of the fermion distribution function. At higher and higher energies, the de Broglie wavelength decreases whereas the parameter χ increases, such that quantum effects remain important during collisions for ultrarelativistic fermions and no radiation reaction force description exists [15]. These parameters are not manifest at the level of the equations of motion, but become accessible by analysis of its solution (see e.g. Secs. V F 2 and V H). In the absence of peaked distribution functions, the medium may not be completely described by a single de Broglie wavelength and no such separation of scales may be identified. In fact, a peaked fermion distribution describes a farfrom-equilibrium situation that does not survive indefinitely in an isolated system. Thereby, systems for which an on-shell Lorentz force description is typically insufficient are those which are initialized with a supercritical field, E E c , and which are then isolated and left on their own. In such systems, fermions are produced from the vacuum -off-shell and at low energies according to Eq. (180) -and then transported in momentum space by the gain-minus-loss structure of the collision terms towards a distribution that is not sharply peaked at any single scale. To describe the evolution towards such a distribution, one requires a description that is valid over a wide range of energies. Thus, the separation of scales from the case of an external field may not be exploited to argue for a Lorentz force description of the equilibration of isolated strong-field systems.
Near equilibrium, a weak field again favors on-shell descriptions, because the field term e / A in the equation of motion of the fermion spectral function (103) then contributes to the transport description only at higher orders and collisions may be added to the on-shell Vlasov equation [see Eq. (191) below] perturbatively in the field vertex (27). However our analysis suggests that for intermediate times, at which off-shell contributions from vacuum pair production equilibrate in the presence of a depleting field, one requires a description of off-shell drifting beyond the Lorentz force. The description derived in Sec. IV can capture this evolution of off-shell contributions inf Ψ (X, p) as they move in phase space towards p 0 = | p | 2 + m 2 to become on-shell particles in the asymptotic future.
It is then instructive to follow how the Lorentz force emerges from the off-shell drift term of Eq. (171), which contains the frequency derivative term ρ µ Ψ (X, p)eF µ0 (X) ∂ 0 pfΨ (X, p) .
As we have shown in Sec. V B, in the asymptotic future the effect of this off-shell frequency derivative is fermion pair production. In the on-shell regime, where pair production is forbidden via Eq. (181), this off-shell frequency deriva-tive is controlled by the dispersion relation, p 0 = ε( p ): The term then contains the group velocity such that, by chain rule, one may replace and recover the classical Vlasov equation Making use of the fact that and applying definitions for on-shell electron and positron distribution functionsf − Ψ andf + Ψ analogously to Eqs. (76), (84) and (85) If we interpret X and p as functions X(λ) and p(λ), then the curves along whichf Ψ is constant, i.e. the characteristic curves d dλf Ψ (X(λ), p(λ)) = 0 , solve the Lorentz equation [50] dp µ dλ = L µ (X, p) , with the Lorentz force (183). Adding collisions that are non-linear inf Ψ makes this method of characteristics inapplicable and the concept of trajectories breaks down. We reiterate that, for general strong fields and fermion distribution functions, the limit of classical propagation (184) is not controlled by an expansion in a small parameter and a combination of the Lorentz force term with the O(e 2 ) collision term (153) is not in general complete to leading order O(e 2 ). To be complete in a general situation, the Lorentz force term should be replaced by the quantum Vlasov term of Eq. (145) (or that of Eq. (171) for small field-gradients).

D. The case of strong external plane-wave fields
We assume in the following that the macroscopic field is of the one-dimensional 'plane-wave' form as originally employed by Volkov [143]. 11 We drop the label 'v' where the context is clear. Assuming (197) means we suppress the parts of the dynamics of the macroscopic field that deviate from a plane-wave field form. The plane-wave approximation is widely used in studying the interaction of laser fields with matter and is valid if the laser beam is not tightly focused in space such that the wave front is approximately flat. Even under such a relatively controlled setup, but especially in isolated systems, one has to take into account that the validity of the planewave approximation can be limited in time. The validity time-scale then depends on the back reaction [53,68] of the matter on the field via Maxwell's equation (29). A simple parametric estimate suggests a large range of validity up to times of t v ∼ O(1/e 2 ). However it is well known [15] that strong macroscopic fields can further decrease this timescale. Below, we assume that the plane-wave approximation is valid for the times under consideration. Although this assumption significantly simplifies the equations, we stress that it does not restrict the discussion of a multitude of common experimental field configurations, such as (linearly or elliptically) polarized fields, (long or short) pulses, monochromatic or polychromatic fields, and (constant or strongly varying) crossed fields.
The field strength tensor of plane-wave fields can be written as where a dot stands for a derivative with respect to the argument. From this it follows that plane-wave fields necessarily satisfy Therefore, the magnetic field B is always perpendicular to and of equal absolute value of the electric field E, such that it is sufficient to only talk about electric fields in the context of plane-waves. In particular, the topological term (200) associated with CP violation [144,145] vanishes identically. This has the implication that the pseudoscalar component of the spectral function (which we introduce in Sec. V D 1, see also appendix F 3) vanishes. 11 Other integrable cases include external fields such as the Coulomb potential (leading to hydrogen levels), homogeneous magnetic fields (leading to Landau levels), constant crossed fields (leading to Airy functions) and constant non-crossed electric fields (leading to Weber parabolic cylinder functions).
Plane-wave systems are most conveniently described using lightcone coordinates that use the special direction n µ of the field, x ⊥ := (x 1 , x 2 , 0) .
Lightcone coordinates have metric tensor η We work in Lorenz gauge [∂ · A(x) = 0] and use the residual gauge freedom to also fix temporal axial gauge [A 0 (x) = 0]. In lightcone coordinates that use the physical direction n µ of the field, this is equivalent (for vanishing asymptotic boundary conditions) to so-called lightfront gauge [146], i.e.
In this frame and gauge, the electric field is simply In particular, this allows for a simple form of the (symmetric) energy momentum tensor from which the energy density of the plane-wave field can be read off. A peculiarity of the plane-wave field is that the classical quantity (208) coincides with the exact vacuum expectation value of the energy momentum tensor up to fermionic contributions [18]. For any function K(X, s − ) of n · s =: s − , one has This is can be written compactly as with the one-dimensional Wigner transform The field dependence enters via the Ritus matrices R q , R q [17,149,150] which are defined as 13 The essential property of the Ritus matrices is that they translate the strong field Dirac operator in position-space into the free Dirac operator in momentum space, i.e.
The plane-wave spectral function contains the strong-field dressed mass [148,153] (see Sec. V F 1) and recovers the Airy-type scattering amplitudes for small field-gradients (see Sec. V F 2). For the proof that (213) solves (103), and satisfies the symmetry constraint (19), as well as for the computation of its Dirac components, we refer to the appendices E and F. The non-perturbative nature of the plane-wave spectral function can be observed from the exponential e iSp : The field-dependent part of the exponent is small for not too strong fields and an expansion in powers of e could be truncated in that case [corresponding to perturbation theory with the vertex (27)]. However for strong fields, A ∼ O(1/e), the exponent is O(e 0 ) and all orders in e, have to be taken into account as depicted in Fig. 7. 12 This plane-wave spectral function ρ Ψ,v is the antisymmetric part of the time ordered 'Volkov propagator' [17,147,148] (see appendix E). By plugging ρ Ψ,v into our transport equations we resum the symmetric part of the fermion propagator to selfconsistent 2-loop order. 13 Sp[X(λ)] is the classical action for the trajectory X µ (λ) of a test particle in a plane-wave field [151]. This fact gives rise to an interpretation of plane-wave scattering probabilities in terms of a stationary phase principle [132,152].
For the application to our transport equations we need the late-time Wigner transform From this expression we can observe that the plane-wave spectral function captures off-shell effects: the external momentum p is not restricted to on-shell values but becomes on-shell in the limit A v → 0, which recovers the free spectral function via With the identity (211) we can discuss the emergence of plane-wave fermion degrees of freedom in strong fields by writing with the field-dependent Dirac matrix and the field-dependent phase factor N q (n · X, n · s) (223) While the phase in terms of S q fully depends on x µ and y µ , the phase N q only depends on n · X and n · s via S q (x) − S q (y) = −q · s − N q (n · X, n · s) .
From (221) we observe that, by the integration over l, the on-shell condition for free fermions (l = 0) is modified to the condition (with l unconstrained) where we have defined the plane-wave relation This expression depends explicitly on the z-component p z := p 3 in which the plane-wave field varies, is positive and satisfies From the context of plane-wave collision terms one further observes in Sec. V E that the parameter l corresponds to the energy exchanged between fermions and the macroscopic field during quantum processes. Since l is integrated over, the relation ε l ( p ) does not on its own restrict the external momentum of the fermion spectral function. Its interpretation as a dispersion relation is thereby not straightforward. Depending on the details of the macroscopic field, the integration over l may have different effects such as broadening the peak structure or adding more peaks. The plane-wave spectral function thereby describes interactions with different l-modes of the macroscopic field, where the lowest mode, l = 0, describes freely propagating particles via In particular, if the macroscopic field is periodic in s − with frequency ω, K(X, l; q) has support only for l = jω with j ∈ Z and a countable peak structure emerges via (see also Ref. [70] for a similar discussion at the level of amplitudes). If l is continuous, the infinitely many delta peaks may merge to form a function with finite-width peaks, such as the function computed in Ref. [49]. As we will see in Sec. V E, the l-modes can be kept track of as individual degrees of freedom by defining appropriate distribution functions that are summed or integrated over in the collision terms. A traditional on-shell description in terms of only the l = 0-mode is then favored as long as a separation of scales in terms of ultrarelativistic fermions exists, as we discuss in Sec. V H.

Collision kernel
Plugging the plane-wave spectral function (218) into (115) we obtain the strong-field plane-wave collision kernel where we have defined the pre-exponential such that, together with the phase (223), the trace over the Ritus matrices becomes We discuss the familiar case of s 1 + s 2 = 0 that emerges in the absence of a medium in Sec. V E 2. For zero field, the phase N p vanishes and the pre-exponential becomes the on-shell amplitude squared s1,s2 such that p → p and q → q as A v → 0. In the presence of a field, the plus-components of p and p, and q and q do not coincide and p + and q + are not on-shell.

Off-shell vs. on-shell kinematics
Classically, particle motion in a plane-wave field [described by the classical Vlasov equation (191)] is characterized by the conservation of the two transverse and the minus-component of the canonical momentum. The plus-component, that is conserved for free particles, is no longer conserved in the presence of a plane-wave field that exchanges energy with particles in this longitudinal direction.
We can derive this interpretation of the field as an energy reservoir from our plane-wave collision kernel also in the off-shell quantum case. By applying identity (211), we may write with the remaining kernel The collision terms therefore contain delta functions enforcing the kinematic conditions where l 1 is the Fourier conjugate to (n · s 1 ) and l 2 to (n · s 2 ). An equivalent set of equations is Equations (245) and (246) make explicit that the physical momenta p, q (carried by the fermion distribution functions) contribute with arbitrary off-shell values, where the 'off-shellness' 2l 1 (n · p ) and 2l 2 (n · q ) is integrated over in the collision terms. In this way, the macroscopic field provides the momenta l 1 n µ and l 2 n µ , preventing the collision terms from vanishing kinematically. Furthermore, the auxiliary momenta p , q are not conserved and k − p + q is not always zero, but corresponds to the energy exchanged with the field, (l 1 − l 2 )n.
In comparison, the zero-field kinematic conditions are which are 'forbidden' because = 0 (on-shell) (253) for massive fermions can only be fulfilled for the trivial case of k = 0, while otherwise Thereby, for vanishing macroscopic field, the delta functions have vanishing overlap and zero-field collision terms vanish at leading order O(e 2 ).
E. Plane-wave photon kinetic equation

Collision term
Employing the plane-wave collision kernel (234), the photon transport equation (111) obtains the following collision term: We can identify the crossings of eeγ scattering depicted in Fig. 8 by taking the frequency integrals over Identifying plane-wave degrees of freedom in terms of the plane-wave fermion and anti-fermion distribution functions and making use of Eq. (227), the plane-wave photon collision term may equivalently be written as with the collision kernels for the different crossings of eeγ scattering 14 14 The sign in the kernels involving one positron recovers the positron term ( / p − m) from the electron term ( / p + m) = −(− / p − m) after the change of sign p → −p (see also Ref. [154]).
As anticipated in our discussion of the plane-wave spectral function (221), we can observe from the energy conserving delta functions in the collision term (259) that the Wigner variables l i correspond to the energy that is exchanged with the macroscopic field by degrees of freedom with energy ε li ( p i ). By means of the changes of variables p z → p z − l 1 and q z → q z − l 2 , this energy exchange can be written in the Lorentz covariant form which clearly relates four-momentum conservation to the structure of the plane-wave field. From the delta functions of the 0 → 3 and 3 → 0 processes in Eq. (259), we observe that they are forbidden for plane-wave fields, since the combination of energy and momentum conditions, can not be fulfilled. In the absence of a macroscopic field, such processes are already forbidden by energy conservation alone. For an arbitrary macroscopic field such a 0 → 3 term would act as a source term for vacuum pair production since it does not come with any distribution function and therefore would not vanish for f + Ψ = f − Ψ = 0. The fact that this contribution vanishes for plane-wave fields is in agreement with the general statement that plane waves are not able to produce pairs from the vacuum [18,155]. In general, such 0 → 3 terms contribute to vacuum pair production at 2-loop order via Eq. (182).

Emergence of a gauge-invariant vertex and gauge-fixing independent amplitude in plane-wave vacuum
Electromagnetic interactions are often described in terms of probabilities for scattering events built from S-matrix amplitudes, which are computed in terms of Feynman rules with free on-shell asymptotic states in vacuum, i.e. vanishing or single mode distribution functions. Such an S-matrix based formulation is not able to resolve real-time dynamics between in-medium states, for which the interaction is not adiabatically switched off. In this section, we follow the emergence of such amplitudes and thereby highlight limitations to their ability to capture collective dynamics of strong-field systems.
General considerations about the trace T , Eq. (231), can be found e.g. in the reviews [16,17,156,157] (see also Ref. [71]) for the special case of s 1 + s 2 = 0 and with X integrated over, which is needed for the computation of probabilities. In this section we identify a scattering amplitude that is local in X and demonstrate that the reduction in terms of relative variables s 1 and s 2 is related to vanishing or single mode plane-wave fermion distribution functions, which we refer to as the 'plane-wave vacuum'. Importantly, such vacuum approximations to distribution functions may only be applied once the relevant degrees of freedom are separated from quantum vacuum fluctuations, because general off-shell distribution functions contain the 'quantum half' which can never physically vanish [see e.g. the constant terms in Eqs. (77) and (79)]. We start in-medium, i.e. without the assumption (267), where the collision kernel (234), may be factorized in terms of Volkov spinors, and written as a 'spin sum' by introducing spin labels σ and σ via By amputating the free Dirac spinorsū pσ , u pσ of U qσ (X − s2 2 )γ µ U p σ (X + s1 2 ) (272) =:ū qσ V µ qp (X, s 2 , s 1 ) u pσ , in Eq. (270), we may identify the vertex This expression differs from the well-known local and gauge-invariant plane-wave vertex [17,158], by its spacetime structure. This difference arises because the local vertex Γ µ pq is constructed from the time-ordered Volkov propagator (see appendix E), which is a vacuum object, i.e. assumes the absence of a medium by vanishing distribution functions, while our vertex V µ pq is constructed in the presence of distribution functions from the antisymmetric part of the Volkov propagator alone. The additional s-dependence of V µ pq , which is integrated over in the collision kernel thus implements the fact that the effective interaction in a strong-field medium is non-local.
While Γ µ pq is gauge-invariant, V µ pq is not. 15 We stress that our collision term is nevertheless gauge-invariant, such that this is not a flaw of our description, but simply exhibits the physical limitations of the concept of scattering probabilities. Electromagnetic interactions in the presence of a medium, i.e. arbitrary distribution functions, can not in general be described by assigning probabilities to individual events. While the photon collision term (259) is gauge-invariant by virtue of additional momentum integrals, the collision kernel and the vertex (273) are not gauge-invariant on their own. Without further assumptions, we can not identify scattering probabilities from them. As we now demonstrate, gauge-invariant amplitudes can be defined in plane-wave vacuum.
First, we investigate how V µ pq reduces to Γ µ pq . For this we need to put the vertex back into the context of the collision term: In general, the photon collision term is of the form p,q δ(k − p + q) g(X, p, q, k) P µν (X, p, q, k) , with the gain-minus-loss term g. If we now assume the absence of a medium, i.e. Eq. (267), there are no other objects carrying fermion momentum dependence other than the kernel itself. We may then write Here, the underlying structure that is simplified by the vacuum assumption is the product of Wigner space fermion spectral functions ρ Ψ (X, p)ρ Ψ (X, q), that can in general not be factorized in real-space via ρ Ψ (x, y)ρ Ψ (y, x) =ρ Ψ (x, y)ρ Ψ (y, x) in the presence of fermion distribution functions, e.g. as in expression (275). However, in the vacuum case, Eq. (276) contains such a factorization, where the δ(k − p + q) has been expressed in real space to invert the Wigner transforms of the spectral functions as in Eq. (134). Eq. (276) then allows us to identify the auxiliary momentum labels p , q of the collision kernel with the physical fermion momenta p, q in the vacuum case. The emerging vacuum collision kernel is gaugeinvariant on its own, and has contributions only from values of s 1 and s 2 satisfying the condition s 1 + s 2 = 0. The momentum labels p and q of the scattering kernel are now on-shell, but there are no fermion distribution functions left. Correspondingly, k − p + q = 0 because momentum is exchanged with the macroscopic field as the amplitude would otherwise vanish kinematically as in the zero-field case. In case of fermion distribution functions that vanish almost everywhere, except e.g. a few ultrarelativistic modes, the dominant contributions from the collision kernel still come from the region of s 1 + s 2 = 0. By taking the collision kernel out of the context of the transport equation in this way, medium effects from more complex fermion distribution functions such as the nonlocal interaction via Eq. (273), and the difference between the on-shell labels p , q and the off-shell labels p, q are missed.
We can now make contact with the language of amplitudes by writing the vacuum collision kernel (277) in terms of the local vertex (274), From this we may read off the local amplitudẽ It is tempting to go one step further and identify the square of the well-known global amplitude [17], by integrating over all X and returning to microscopic position variables via However, it is important to note that such an integration over all times X 0 generally includes late times outside of the range of validity of both the plane-wave field and the plane-wave vacuum approximation. Even if one makes assumptions such as (197) and (267) at initial time, the macroscopic field does not remain a plane wave and distribution functions do not in general remain negligible, but backreact on the field, such that a self-consistent description of both becomes essential to describe equilibration.
What type of interactions take place in a strong-field system depends also on the details of distribution functions and is a time (and space) dependent questions. To determine this time dependence, one has to solve the transport system including the dynamics of distribution functions away from the plane-wave vacuum. Instead, a common approach in literature is to rely on the S-matrix in the Furry picture [147], which takes amplitudes out of the context of in-medium evolution equations. To then extract local probabilities from this S-matrix, the LCFA is a necessary approximation, as otherwise what is supposed to be the local probability density may turn out to be negative (see e.g. Ref. [71]). This problem occurs because the probability that a scattering will take place in an external field and in the absence of a medium at any time is not a self-consistent concept in general. The gaugeinvariant amplitudes such as (280) are not observable and probabilities for individual scattering processes need not exist to compute statistical observables such as electrical conductivity [79] or pressure [78]. No further approximations such as the plane-wave vacuum or locally constant fields are necessary to compute also e.g. the photon decay rate (116) from the equations discussed in Sec. IV. While the physical interpretation of the strong-field amplitude (280) is problematic, that object is very useful to understand the ξ-dependence of our equations. The amplitude (280) is known to obey the modified Ward identity [131][132][133] k · M σσ =ū pσ / nu qσ with the phase relating gauge-fixing to boundary terms at n · x = ±∞.
Vanishing boundary terms then lead to gauge-fixing independence, P ξ ≡ 0.

F. Plane-wave fermion kinetic equation
Because the fermion collision term (153) relies on the gauge-invariant fermion spectral functionρ Ψ (as opposed to the covariant function ρ Ψ ), we start this section by investigating this function for plane-wave fields. The wellknown plane-wave momentum and dressed mass emerge automatically in this function. These gauge-invariant expressions then serve us to approximate field-gradients in a gauge-invariant manner, equivalently to Sec. V A, but at the level of the solution rather than the equation of motion.

Gauge-invariant spectral function: plane-wave momentum & dressed mass
The covariant plane-wave spectral function (221) transforms as any other fermion two-point function. The ambiguity [50,56,82] for the choice of the path of integration in the Wilson line is not present in the plane-wave case because there is only one path in one dimension from n · x to n · y. Thereby, the Wilson line automatically emerges with a straight path of integration, dλ a(n · X + λ) (285) (not to be confused with the ensemble average (4)) for any plane-wave function a(n · x), we can make the gaugeinvariance ofρ Ψ,v manifest. By employing Eq. (285), the plane-wave Wilson line (284) can be written as The Lorentz equation for plane-wave fields is solved by the gauge-invariant momentum of an electron in a plane-wave [151] π µ q (n · X) (287) which is related to the Lorentz action (215) via The plane-wave momentum obeys π 2 q = q 2 and n · π q = n · q and is related to the free mass m and the gauge-invariant dressed mass [148,153] m 2 (X, s − ) := m 2 − e 2 (n · s) n·s 2 − n·s 2 dλ A 2 (n · X + λ) via [148] π 2 q = m 2 and π q 2 =m 2 (291) for any q with q 2 = m 2 (which in our context is ensured by the delta function under the integral e.g. in Eq. (293)). We can identify this plane-wave momentum in the exponent of the gauge-invariant spectral function via such that an exact solution of Eq. (138) in the plane-wave case may be written aŝ with the gauge-invariant Dirac matrix .
While the covariant spectral function (221) makes manifest the energy exchange with the field and facilitates a formulation in the plane-wave degrees of freedom (257) and (258), the invariant function (293) makes manifest the solution of the Lorentz equation (287). The scalar and pseudoscalar components ofK q are The vanishing of the pseudoscalar component is a direct consequence of the crossed nature of plane-wave fields, i.e. Eq. (200). The vector component, which plays a crucial role in the quantum Vlasov term, contains the plane-wave momentum also in the pre-exponential and is given in Sec. V F 3. The axial and tensor components can be found in the appendix F. The tensor and scalar components vanish for massless fermions in agreement with chiral symmetry. Similarly to the identity (211) one has Computation of the scalar component (see appendix The corresponding symmetric component has been computed in Ref. [49] (see also Ref. [21]) for various choices of plane-wave fields, in the context of scalar QED.

Plane-wave fields with small gradients
In this section, we investigate the gauge-invariant approximation of field-gradients using the example of the scalar spectral component (298).
For plane-wave fields, the gradient expansion becomes an expansion in longitudinal gradients via where A (j) µ is the jth derivative with respect to n · X. In the scalar component (298), field-gradients are carried only by the gauge-invariant mass whose expansion is gauge-invariant order by order [148].
Similar to the fact that the equation of motion (167) has contributions from constant gauge-invariant fields, the second term of the dressed mass is also non-trivial for constant electric fields and generally not small compared to unity. In fact, introducing the dimensionless and Lorentz-invariant quantities with the Schwinger critical field E c = m 2 /|e| and a characteristic field amplitude F 0 and frequency ω, we can write the constant-field contribution from this exponent as −e 2 1 12Ȧ 2 (X) 2(n · p) (n · s) 3 = 1 24 with the notation E(X) := | E(X)|. Equation (304) reveals the significance of ξ 3 0 /χ 0 for locally constant fields, which is well known in laser physics [159][160][161]. All higher order gradient contributions to the exponent of the gaugeinvariant spectral function from the dressed mass, e.g. the next order terms e 2 1 720 are suppressed by gauge-invariant gradients. Similarly, one may explicitly verify that under this locally-constant field approximation, the Wilson line can be approximated in such a way that the relation between covariant and invariant fermion two-point functions, Eq. (129), and the relation between f Ψ -andf Ψ -type fermion distribution functions, Eq. (131), indeed holds.
Keeping the LO of the dressed mass, we find the gaugeinvariant LO scalar component The ϕ integral leads to the Airy function 16 where the local parameter χ defined in (186) amounts to χ 0 with F 0 replaced by E(X) for plane-wave fields (for the computation of Eq. (307) see appendix F 1). The LCFA strong-field scattering probabilities [149,150] that are used as input in the kinetic equations e.g. of Ref. [22] also feature such Airy functions. As anticipated in Sec. V C, these functions may be further reduced to on-shell delta-peaks by virtue of the identity lim χ→0 1 χ Ai(x/χ) = δ(x), consistent with a classical radiation reaction regime.

Quantum Vlasov term
To discuss the quantum Vlasov term for small fieldgradients of Eq. (171) for plane-wave fields it is useful to switch to lightcone coordinates, For plane-wave fields, the lightcone components of the Vlasov derivative simplify to with ∂/∂X − = (∂/∂X 0 − ∂/∂X 3 )/2 and ∂/∂X + = ∂/∂X 0 + ∂/∂X 3 , and analogous definitions for momentum derivatives. A p − derivative is absent as it comes with F µν n ν , which vanishes for plane-wave fields.
The all-order in field-gradients plane-wave spectral vector component iŝ whereπ q is the plane-wave momentum in the field 1 2 [A µ (x) + A µ (y)] explicitly stated by Eq. (F17) in the appendix. The computation of the pre-exponential makes use of the fact that and can also be found in the appendix F 2, alongside the leading order in field-gradients. These are all the ingredients one needs for the quantum Vlasov term for locally constant plane-wave fields. In principle, with Eq. (312) available, the drift term of the all-order fieldgradient equation (145) is also accessible. The lightcone components of the vector spectral function (312) obeŷ These identities are particularly useful in the classical radiation reaction regime (χ → 0), to recover the on-shell Lorentz force drift term of the classical Vlasov equation (191) by expressing the vector component in terms of the scalar component and then using lim χ→0 1 χ Ai(x/χ) = δ(x). In systems with a long-lived separation of scales in terms of ultrarelativistic fermions, it is possible to reduce Eq. (308) to a Lorentz force term also without sending χ → 0 (see Secs. V C and V H).

Electron and positron collision terms
Inserting the plane-wave collision kernel (234) into the fermion collision term (153) and making use of identities (131) and (154) we may write the fermion collision term for small field-gradients as Here, we have relied on small field-gradients to write the invariant spectral function of the fermion collision term (153) in terms of the plane-wave delta function of the covariant specrtal function (221), Defining electron and positron collision terms, the frequency delta functions in Eq. (316) then allow for explicit computation of the frequency integrals and to recover the structure in terms of the strong-field scattering processes depicted in Fig. 8.
The sign in the definition (319) accounts for a sign that arises when substituting p → − p. The factors of 1 2 account for the absence of a factor 2 in the identity for the firstorder derivatives of fermions i(/ ∂ x +/ ∂ y ) = i/ ∂ X as compared to the identity for the second-order d'Alembertians for photons (98).
The appearance of p+eA in Q is resolved in the vacuum limit, where scattering kernels become gauge-invariant on their own as discussed in Sec. V E 2. Since the fermion self-energy is not gauge-invariant, the emergence of gaugeinvariant scattering amplitudes with no Wilson lines is far from obvious. However, in the vacuum case, a gaugeinvariant fermion loop emerges from the product of the fermion self-energy and the fermion propagator under an additional momentum integral. The ultrarelativistic limit discussed below in Sec. V H then resolves any re-maining obstructions to a description in terms of on-shell distribution functions.

G. The case of small occupations
The complexity of collisional kinetic equations is largely due to the nonlinearity in distribution functions of collision terms. However, many physical situations allow for an assumption of small distribution functions, i.e.
implicit for example in the kinetic equations of Refs. [22,43]. For such settings close to vacuum, one may drop 2 → 1 and 3 → 0 processes entirely, since they contain no linear terms and are therefore suppressed 17 by Thereby, strong-field systems with small occupations single out a direction in time -the direction of energy transport from the macroscopic field to the particle sector by 1 → 2 and 0 → 3 processes -even if the corresponding scattering matrix elements and the fundamental equations of motion are symmetric under time reversal. Similarly, one may simplify all Bose-enhancement or Pauli-blocking terms in 1 → 2 and 0 → 3 processes via 1 + f (X, k) ≈ 1 and 1 − f ± Ψ (X, l, p ) ≈ 1. In this way, small distribution functions lead to a linearization of collision terms. In contrast to a linearization around equilibrium [104] which keeps thermal distributions as in Eq. (118), collision terms linearized by small occupations violate detailed balance and are thereby no longer able to describe the approach to thermal distribution functions. Charge conservation [Eq. (157)] is still exact.
The linearized near-vacuum plane-wave photon collision term (259) then reads [after a substition to recover covariant energy conservation as described around Eq. (264)], × δ(k − p + q − (l 1 − l 2 )n) Q e + →e + γ (X, l 2 , l 1 , q + l 2 n, p + l 1 n ) f + Ψ (X, l 1 , p + l 1 n ) The electron collision term (318) under the same approximation reduces to the three 1 → 2 scattering processes, e − → e − γ with ingoing momentum p , e − → e − γ with outgoing momentum p , and γ → e − e + . Analogously, the linearized positron collision term (319) contains the processes e + → e + γ and γ → e + e − . In all near-vacuum collision terms, each process is weighted linearly by the distribution function of the ingoing particle as in the equations of Ref. [22].
We emphasize that for general macroscopic fields, these near-vacuum collision terms would all additionally contain 0 → e + e − γ source terms with no distribution function, contributing to vacuum pair production at 2-loop O(e 2 ) precision.

H. The case of ultrarelativistic fermions & on-shell strong-field descriptions
Many of the approximations discussed in previous sections are tied together in an ultrarelativistic setting: strong macroscopic fields accelerate fermions to ultrarelativistic energies within small regions of space. Once accelerated, any macroscopic field appears like a planewave field in the Lorentz rest frame of an ultrarelativistic fermion [162]. Therefore, plane-wave fields represent generic qualities of strong fields in an ultrarelativistic setting. Furthermore, ultrarelativistic fermions facilitate chiral symmetry, which in turn leads to a reduction of tensor structures, which is assumed by our definition of the fermion distribution function as discussed in Sec. III C. Additionally, large fermion momenta can facilitate that field-gradients are numerically separated from propagatorgradients as we have seen in Eq. (61). Moreover, ultrarelativistic fermions have a small de Broglie wavelength facilitating classical propagation in-between quantum processes like the emission of photons. From an analysis of the classical propagation of fermions one then finds that ultrarelativistic fermions emit radiation along their instantaneous velocity, within a cone of angular aperture ∼ m/ε( p ) [15,138]. If the particle is ultrarelativistic and its energy is the largest scale in the system, its motion has a pronounced directionality. In strong-field vacuum, i.e. for vanishing occupations, and if the transverse momenta are much larger than m, one can then show that only small patches of their trajectory contribute to scatter-ing amplitudes [15,152] (which is the assumption of the LCFA).
There are several notions of ultrarelativistic limits for fermions in literature. They range from assumptions on kinematic restrictions [22] to expansions in terms of p ⊥ /(n · p) [162] or 1/γ = m/ε( p ) [15]. In the language of the present paper, an ultrarelativistic system is defined by a fermion distribution function that is peaked at an ultrarelativistic scale p * . Such a distribution function then gives meaning to single particle concepts such as the de-Broglie wavelength /p * also in many-body systems. In particular, the structure of the fermion spectral function of such a system only matters for characteristic momenta as it always appears in a product with the fermion distribution function which approximately vanishes away from the characteristic scale.
We now assume that the characteristic scale p * of fermion distribution functions is well separated from the characteristic scale l * of the strong-field spectral kernel (222), i.e.
In such a situation, the strong-field spectral function only contributes with on-shell values, because becomes the on-shell dispersion, independent of l * . As anticipated by our discussion in Sec. V C, this implies that ultrarelativistic fermions may indeed be described by on-shell particles whose energy ε( p * ) then satisfies The ultrarelativistic limit (324) leaves the strong-field properties of the spectral function intact, simplifies kinematic restrictions, and favors a description in terms of free distribution functions. However, in general, there is no mechanism that dynamically controls this approximation, i.e. it may become invalid during the evolution of the system even if it is valid at initial time. An important effect that explicitly breaks the validity of an ultrarelativistic approximation is vacuum pair production, which generates off-shell contributions tof Ψ (X, p) at zero frequency according to Eq. (180).
Indeed, a kinetic description in terms of only on-shell distribution functions is suggested in Ref. [22] for ultrarelativistic fermions in strong (but subcritical E E c ) fields with small gradients. Our off-shell transport description of Sec. IV reduces to that description under the following approximations: a) an approximation of fieldgradients (see Sec. V A); b) an approximation of collision terms for small occupations to neglect medium effects (see Secs. V E 2 and V G); c) an assumption of ultrarelativistic simply peaked fermion distribution functions and subcritical fields to replace the quantum Vlasov term with the Lorentz force term of the classical Vlasov equation (193) (see also Sec. V C) and to justify the on-shell limit of collision terms, dl 1 dl 2 g(l 1 , l 2 )Q(l 1 , l 2 ) (326) where g indicates the gain-minus-loss terms and Q includes the delta-functions such that energy conservation is treated exactly. Together with the Lorentz force term, this closes the strong-field description in terms of the traditional on-shell particle distribution functions (84) and (85) which emerge from plane-wave distribution functions via With all these approximations combined, also subtleties regarding gauge-invariance both of the scattering kernels and the distribution functions are resolved: The scattering kernels become gauge-invariant objects in the vacuum limit (274) and a distinction between thef Ψ -and f Ψ -type fermion distribution functions is not important after the ultrarelativistic limit (327) for distribution functions that are always only occupied in terms of a few on-shell particle modes for which gauge-invariance is assured. Dropping all these assumptions is possible by employing the gauge-invariant off-shell equations discussed in Sec. IV. Starting from this off-shell description, it would be interesting to investigate whether collisional or inhomogeneous contributions to the particle yield (182) can invalidate an on-shell description also for subcritical fields on some significant time scale on the way to equilibrium.

VI. CONCLUSIONS & OUTLOOK
Our work demonstrates how to systematically derive transport and kinetic equations including collisions for general supercritical fields. The equations to order O(e 2 ) include local scattering kernels for strong fields that can also be inhomogeneous. This is achieved by off-shell transport equations that include non-local relative times and all field-gradients in the fermion spectral function, while retaining the gain-minus-loss structure of traditional kinetic equations. To investigate our equations analytically and to make contact to limiting cases in the literature, we have also considered plane-wave fields.
We have shown that the inclusion of fermion spectral dynamics is essential to describe collisions and fermion drifting in the presence of general strong fields. Existing derivations of strong-field Wigner descriptions in the literature have neglected spectral dynamics by limiting themselves to the collisionless regime, in which equations for spectral functions decouple from transport equations. In general, however, the macroscopic field enters the collision kernel (115) via the fermion spectral function (103). This resums infinitely high perturbative orders of the coupling that all become relevant for sufficiently large macroscopic fields. The macroscopic field itself is governed by a Maxwell equation in the presence of a fermion current involving the quantum corrections. The general form of this Maxwell equation turns out to be valid to arbitrary loop and gradient order in our framework. Our approach paves the way for investigations of the thermalization process starting from strong field initial conditions, which requires to go beyond collisionless approximations.
We have pointed out a connection between asymptotic pair production and spectral dynamics. While 1-loop results such as the Schwinger pair production rate (173) assume the macroscopic field to be external and constant in time, our 1-loop result (180) is fully dynamical and generalizable to the expression (182), which in principle includes collisions to 2-loop order and all orders in field-gradients. Our description in terms of distribution functions does not rely on asymptotic expressions, such as total particle numbers or total probabilities in order to compute time-dependent observables such as the strongfield photon decay rate (118).
We solved the LO equation for the fermion spectral function for the special case of an external plane-wave macroscopic field, A µ (x) A µ v (n · x) with a null vector n 2 = 0. This reduces the transport description to only two equations for the off-shell fermion and the on-shell photon distribution function. The plane-wave spectral function is the antisymmetric part of the well-known timeordered Volkov fermion propagator. By employing only its antisymmetric part in the O(e 2 ) transport equations, we self-consistently resum quantum fluctuations to 2-loop order. Thereby, a solution of our equation goes beyond the statistical component of the 1-loop Volkov propagator that implicitly assumes vanishing distribution functions.
Employing the all-order field-gradient plane-wave spectral function in the collision kernel reproduces expressions which are similar to the Furry picture, but have the advantage of being automatically local in the kinetic position variable X while containing contributions from inhomogeneous fields not limited to the vicinity of X. In particular, we have demonstrated that plane-wave scattering kernels emerge with a space-time structure that is more general than the one of local scattering amplitudes that are known from laser applications. The more general scattering kernels reduce to known expressions only if relative times are restricted to certain values. We have recognized this condition as the implicit assumption that the system is in plane-wave vacuum, i.e. that fermion plane-wave distribution functions are negligible or have only single occupied modes. This means that medium effects are missed if a strong-field scattering kernel is obtained from Feynman rules for the Furry picture S-matrix. For negligible distribution functions, known gauge-invariant global scattering amplitudes emerge by integrating over all X. To employ external-field vacuum amplitudes in an isolated dynamical setting is typically inconsistent because it includes times outside the range of validity of external field and vacuum approximations as non-negligible distribution functions develop dynamically and backreact on the field. Nevertheless, these emergent amplitudes allowed us to highlight connections to Ward identities, which remove the gaugefixing dependence of the 2PI formulation of QED in the corresponding limit.
Furthermore, the plane-wave fermion spectral function allowed us to use the energy exchange with the macroscopic field as a parameter l to label strong-field degrees of freedom with energy ε l ( p ), which enable a continuous connection to the free particle degrees of freedom of an on-shell description: When this spectral function is multiplied by a fermion distribution function that is peaked on an ultrarelativistic scale p * that is well separated from the characteristic value of l, its dispersion relation becomes independent of l and reduces to that of free fermions, ε l ( p * ) | p * | 2 + m 2 . This facilitates an on-shell description despite the presence of strong fields. Thereby, strong-field systems in which such a clearly separated scale p * exists for long times may be accurately captured by on-shell descriptions that combine collisions with classical Lorentz force drifting. Since any field appears as a plane wave in the rest frame of a single ultrarelativistic fermion, we expect that most statements that we arrived at under the assumption of an external plane-wave field also hold for more general fields, as long as fermion distribution functions are dominated by a few ultrarelativistic particle-modes.
In contrast, in isolated systems with supercritical fields, initial characteristic scales are dynamically affected by pair production (which occurs off-shell and is largest at zero frequency) and by the transport of fermion occupations towards an equilibrium distribution (which has its maximum at low energies and is not sharply peaked). For such isolated systems, we argued that an initial ultrarelativistic separation of scales is not long-lived, such that an on-shell Lorentz force description introduces an error larger than our desired accuracy of O(e 2 ). In the absence of a long-lived separation of scales, one needs instead a description that remains valid over a wide range of energies to describe the evolution of off-shell contributions induced by vacuum pair production towards the on-shell regime of the asymptotic future. The gauge-invariant fermion transport equation Eq. (145) with its all-gradient off-shell drift and collision term constitutes such a description by coupling to the fermion spectral equation, the photon transport equation and the Maxwell equation summarized in Fig. 6.
Outlook. The leading order equations may give insight into the largely unexplored late-time behaviour of isolated QED systems with finite net charge. If such a system equilibrates, its late-time state can not be the traditional homogeneous thermal equilibrium, because the Gauss constraint for finite net charge prevents the initial field both from decaying completely and from becoming fully homogeneous. The possible approach to such a charged time-translation invariant state may be completely described by our equations, if the equilibrium field induced by the net charge turns out to be sufficiently large.
Such a numerical computation, in particular of the self-consistent strong-field fermion spectral function, will also allow for a more detailed study of the collision kernel and the spectral peak structure. This would, e.g., enable an analysis of spectral widths and to establish under what circumstances they are small. To obtain insight into specific controlled experimental settings, one may employ other external fields in such a computation, as we have done for the plane-wave spectral function with laser fields in mind. Possible other choices of external fields include non-crossed constant electric fields, homogeneous magnetic fields and Coulomb fields.
In the future, dropping our assumption of reduced tensor structures (74) with the help of Ref. [50] could clarify the significance of chiral dynamics [115][116][117][118] and spin transport [124], and extend chiral kinetic theory [119][120][121][122][123] to the collisional regime. To access the transport dynamics of the axial current j µ 5 (X) = −etr{γ 5 γ µ F Ψ (X, X)}, an interacting spectral function that has a non-vanishing axial component such as the strong-field spectral function employed in this paper is required (see the expression for plane-wave fields in appendix F 4). A leading-order collisional description including all tensor structures in the presence of a macroscopic field is now in reach and would open up diverse applications on chiral dynamics reaching from astrophysics [163,164] to semiconductors [165,166].
Acknowledgements. This work is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -Project-ID 273811115 -SFB 1225.
In Wigner space, one may alternatively exploit the Wigner transform of the Heaviside function, to obtain a Källén-Lehmann representation and similarly for the advanced functions with The self-energies obey completely analogous identities. We stress that any singularity associated to the εprescription does not arise in an exact (early-time) description that employs Wigner transforms (46) instead of the late-time Wigner transforms (47) (see also Refs. [105,106]). The photon Wigner functions have the properties F µν (X, k) = F νµ (X, −k) , (A23) ρ µν (X, k) = −ρ νµ (X, −k) , Similarly, fermion Wigner functions obey Given all this, it should be remembered that there are only two independent two-point functions per field species (see also our comment at the end of Sec. III A). The LO O(e 2 ) 2PI loop expansion of the Wightman self-energies is where we used that θ(x − y) + θ(y − x) = 1 and that on the backward branch of the contour, x, y ∈ C − , δ C (x, y) = − 1 2 δ(x − y), which results in in the fact that only the statistical function F Ψ (x, y) = 1 2 (∆ +− (x, y) + ∆ −+ (x, y)) contributes to this term (see also Ref. [65]). The variation of this term then gives the Maxwell current We now change to the gauge-invariant statistical propagator by introducing Wilson lines. For small fields with small gradients, we may expand the straight Wilson exponent via x y dz µ A µ (z) = s µ ∞ n=0 1 (2n + 1)! 1 2 2n (s · ∂ X ) 2n A µ (X) . (D4) The leading order of straight Wilson lines, W(x, y) = e −ies·A(X) + O (e 0 s · ∂ X ) 2 (first order vanishes), is O(e 0 ) for strong fields and produces the missing term that is necessary to identify the gauge-invariant field strength tensor. Changing the prescription (D3) how to derive fermion kinetic equations to include a Wilson line, we recover the Vlasov term via d 4 (x − y) e ip(x−y) W(y, x) 1 4 tr (F LHS) Ψ (x, y) − γ 0 (F LHS) † Ψ (y, x)γ 0 (D6) Appendix E: Symmetric and antisymmetric parts of the Volkov propagator By virtue of Eq. (217), (i/ ∂ x − e / A v (x) − m)E q (x)( / q + m)Ē q (y) = E q (x)(q 2 − m 2 )Ē q (y), such that indeed the plane-wave spectral function solves The Volkov spectral function is antisymmetric because γ 0 ρ Ψ,v (x, y) † γ 0 = −i(2π) q δ(q 2 − m 2 )sgn(q 0 )γ 0Ē † q (y)( / q + m) † E † q (x)γ 0 = −ρ Ψ,v (y, x) .
We may put this spectral function into the context of canonical quantization in the Furry picture [147], which is achieved for plane-wave fields by means of the Volkov states [148] Ψ via the canonical commutation relations for the ladder operators with c s ( p ) |0 v = 0: The plane-wave spectral function can be written as the expectation value of the anticommutator of Volkov states with respect to the strong field vacuum |0 v , such that The symmetric part of the Volkov propagator (that has been discussed e.g. in Ref. [49]) is the commutator F Ψ,v (x, y) = 1 2 Ψ v (x),Ψ v (y) = π d 4 q (2π) 4 δ(q 2 − m 2 ) θ(q 0 ) + θ(−q 0 ) R q (x)( / q + m)R q (y) .
The Volkov propagator [148] is then built via the standard asymptotic state identity where T denotes ordinary time-ordering.
Since here the field appears only under the absolute value, this is simply the sign function sgn(q − ) familiar from lightcone quantization and compensates the absolute value in the identity (F4). The exact scalar plane-wave spectral function in position-space thereby is ρ Ψ,v,S (X + s 2 , X − s 2 ) = i m (2π) q δ(q 2 − m 2 )sgn(p 0 ) e −i[q·s+Nq(X,s)] (F6) With this we can immediately identify the gauge-invariant part via ρ Ψ (x, y) = W(x, y)ρ Ψ (x, y). Without gradient expansion, the Wilson line is exactly canceled and no additional substitution of p → p + eA is necessary. The Wigner transform can easily be computed up to the s − integral viâ where we have used thatm only depends on s − not on s + , s ⊥ . The result (298) mentioned in the main text, follows after taking the trivial integrals over the delta functions. Next we prove that the leading order in gauge-invariant field-gradients of the function (306) is equivalent to the Airy expressions (307). For this purpose we make use of dϕ e iaϕ−bϕ 3 = (3b) −1/3 Ai − a(3b) −1/3 .