Resummed hydrodynamic expansion for a plasma of particles interacting with fields

A novel description of kinetic theory dynamics is proposed in terms of resummed moments that embed information of both hydrodynamic and non-hydrodynamic modes. The resulting expansion can be used to extend hydrodynamics to higher orders in a consistent and numerically efficient way; at lowest order it reduces to an Israel-Stewart-like theory. This formalism is especially suited to investigate the general problem of particles interacting with fields. We tested the accuracy of this approach against the exact solution of the coupled Boltzmann-Vlasov-Maxwell equations for a plasma in an electromagnetic field undergoing Bjorken-like expansion, including extreme cases characterized by large deviations from local equilibrium and large electric fields. We show that this new resummed method maintains the fast convergence of the traditional method of moments. We also find a new condition, unrelated to Knudsen numbers and pressure corrections, that justifies the truncation of the series even in situations far from local thermal equilibrium.


I. INTRODUCTION
Relativistic hydrodynamics plays a fundamental role in the description of a wide range of physical phenomena, from astrophysical plasmas to heavy-ion collisions [1][2][3]. There are two main theoretical frameworks for deriving hydrodynamics from an underlying microscopic theory.
The first one originates from the Chapman-Enskog expansion [4], which involves a systematic power counting of the gradients of the standard hydrodynamic quantities, i.e. temperature, chemical potential, and velocity fields. An appealing feature of this approach is that the gradient expansion can be performed even around a quantum field theoretical local equilibrium background, i.e. it does not fundamentally require a classical approximation. Such an expansion can also be done in the relativistic regime, both in the context of kinetic theory [5] and in the case of strongly coupled relativistic systems, as shown in Refs. [6,7].
The straightforward relativistic generalization of the gradient expansion truncated at first order leads to the relativistic Navier-Stokes equations [8]. These equations violate causality [9] and are linearly unstable around equilibrium [10,11] (see also [12,13]). The series can be extended by including second order gradient corrections [6,7] and also third order terms [14], though precise statements regarding causality in the nonlinear regime and stability are not available in those cases. Mathematically rigorous results about causality and stability in relativistic viscous hydrodynamics were presented in Ref. [15] where it was shown that the gradient expansion can be used to derive a causal and stable theory at first order in gradients if different definitions for the hydro-dynamic fields are used. However, recent work showed that the gradient series has zero radius of convergence [16][17][18][19]. Hence, including higher order terms in the expansion does not constitute a viable path towards a systematic improvement of the description of the system in the far-from-equilibrium regime.
The second method widely used in the derivation of relativistic hydrodynamics employs a moment expansion of the relativistic Boltzmann equation [20]. While Boltzmann kinetic theory is only applicable to sufficiently weakly coupled ("dilute") systems, this approach exploits the idea that hydrodynamics can be understood as an effective field theory describing the long time, long distance macroscopic behavior of the system [7] whose structure is universal. Therefore, details of the microphysics enter into the hydrodynamic equations only through material properties (such as the equation of state and the transport coefficients) whose calculation involves different methods at weak [21] and strong coupling [22].
In this approach the Boltzmann equation is expressed in terms of an infinite set of coupled equations for the momentum moments of the the distribution function. Since (at least in the relaxation time approximation, see below) every moment couples only with a finite number of other moments, one can truncate the set of equations at some order, using some approximation for the leftover moments [20]. In the absence of long range mean fields, this method can be systematically improved by the inclusion of more moments. Under flow conditions of extreme symmetry where the relativistic Boltzmann equation can be solved exactly [23][24][25][26][27], this procedure has been shown to converge rapidly to the exact results of relativistic kinetic theory [18,28,29]. The reasons behind this rapid convergence even in far-from-equilibrium conditions are, however, still poorly understood as, at first sight, all moments seem to contribute to the equations with similar weight.
The purpose of this work is explore these questions in greater depth, by extending the method of moments to a more general microscopic background. By generalizing the Boltzmann equation to Boltzmann-Vlasov form we introduce long-range non-collisional forces that could stem from of an electromagnetic gauge field or from medium-dependent particle masses. In the case of electromagnetic interactions the gauge field is calculated self-consistently by taking into account the contribution of the particle current density to Maxwell's equations. We develop the moment expansion of the coupled set of Boltzmann-Vlasov-Maxwell (BVM) equations for a general collision term in 3+1 dimensions, before solving the theory exactly (i.e. with arbitrary numerical precision) in the Relaxation Time Approximation (RTA) [30] in 0+1 dimensions for a system undergoing boost-invariant longitudinal expansion without transverse flow (Bjorken expansion [31]). This paper is organized as follows: In Section II we briefly review the canonical approach to the method of moments for the Boltzmann equation. In Sec. III we introduce the BVM equations and discuss its expansion in terms of moments. In Sec. IV we propose a new set of moments that resum the contributions of an infinite set of non-hydrodynamic degrees of freedom and show how they can be used to investigate the physics behind the BVM equations efficiently even for massless particle systems. An exactly solvable case of these equations is studied in Sec. V and compared in Sec. V C with results from the expansion in terms of the resummed moments. We find fast convergence of the resummed moment expansion to the exact solution. Conclusions and a final overview are presented in Sec. VI. Various technical discussions are relegated to Appendices A-E.
Unless otherwise stated we use natural units where = c = k B = 1. We adopt the Einstein convention of automatically summing over repeated upper and lower indices, and we represent the contraction (scalar product) between four-vector with a dot: v µ w µ = v · w. The "mostly minus" convention for the Minkowski metric is used, i.e. g µν = diag(1, −1, −1, −1), as well as the convention ε 0123 = 1 for the four-dimensional Levi-Civita symbol. Round parentheses around groups of Lorentz indices indicate symmetrization, e.g. A (µ B ν) = 1 2! (A µ B ν + A ν B µ ).

II. THE METHOD OF MOMENTS
The method of moments, initially introduced in the non-relativistic regime by Grad [32], has been widely used to study a number of properties of the relativistic Boltzmann equation [33] where f (x, p) is the single-particle phase-space distribution function of a gas of particles with mass m, and C[f ] is the collision term which, in general, involves an integral over the distribution function f . In the moments method one replaces the integrodifferential mathematical problem defined by the Boltzmann equation by an infinite set of coupled partial differential equations for the momentum moments of f , which correspond to macroscopic quantities such as, for instance, the fluid's energy-momentum tensor. For the moment we will use Cartesian coordinates such that we do not need to distinguish between the partial derivative ∂ µ and the covariant one d µ [34]. For any definition of the four-velocity u µ , we can split ∂ µ into the time and spatial derivatives in the comoving frame, where D = u · ∂ is the comoving time derivative, also denoted by a dot (Ȧ ≡ DA ≡ u µ ∂ µ A for any quantity A), and ∇ µ ≡ ∆ µν ∂ ν (where ∆ µν = g µν −u µ u ν projects on the spatial coordinates in the comoving frame). Using these definitions, the Boltzmann equation (1) can be rewritten as follows [35] ( from which exact equations of motion for the (reducible) tensor moments with r being an integer and s a non-negative integer, can be derived: The contribution from the collisional kernel is given by To obtain Eq. (5) one only needs uniform convergence of the momentum integrals (4). In (4), (6), and similar integrals below, p indicates the Lorentz invariant momentum integral with g counting the degeneracy of each momentum eigenstate. In the last expression the integrand must be taken on-shell, i.e. p 0 = m 2 +p 2 . The particle number current N µ = p p µ f and energy-momentum tensor T µν = p p µ p ν f are given by the moments F µ 0 and F µν 0 , respectively. For (r, s) = (0, 2), Eq. (5) gives the following exact evolution equation for the stress-energy tensor: Its projection onto the four-velocity u µ yields four equations describing energy-momentum conservation, ∂ µ T µν = 0 (see Appendix A). The remaining six equations provide the exact evolution of the dissipative corrections to the pressure (bulk viscous pressure and shear stress). Unlike the four-momentum conservation equations, the latter couple the components of the stressenergy tensor with other moments of the distribution function.
Since in the comoving frame only the spatial projections of the T µν couple directly to non-hydrodynamic moments, it is convenient to introduce the following notation for the spatial components of the tensor moment in the local rest frame (LRF) where u µ = (1, 0): Angular brackets around a tensor index indicate its spatial components in the LRF obtained by contracting the four-index with the spatial projector ∆ µ ν : p µ ≡ ∆ µ ν p ν . We note that the tensor moments F are not mutually independent. In fact, the projection with u α of a tensor moment of rank (r, s+1) produces a tensor moment of different rank (r+1, s), u α F αµ1···µs r = F µ1···µs r+1 , while projecting all upper indices along the four-velocity yields a scalar moment u µ1 · · · u µs F µ1···µs r = F r+s = f r+s . On the other hand, no analogous relation exists for their spatial components in the LRF, i.e. for f µ1···µs r .
It is useful to rewrite the exact evolution equations (5) in terms of f µ1···µs r . As shown in Appendix B, their exact evolution equations arė where θ ≡ ∇ µ u µ is the scalar expansion rate. These moment equations contain the information we need in this work. For example, by projecting Eq. (8) with u µ u ν we obtain Eq. (10) with r = 2 and s = 0: where we used in the first line that p p µ C = 0 and in the second line that u µ f µ 1 = 0. This describes the conservation of energy. Projecting (8) with u µ ∆ α ν yields the momentum conservation laẇ where ∇ µ f µ α ≡ ∆ α ν ∇ µ f µν . Making use of the general hydrodynamic decomposition of the energy-momentum tensor in the so-called Landau frame [8] Eqs. (11) and (12) can be rewritten in the more familiar formĖ E+P+Π u α = ∇ α P+Π − ∆ α µ ∂ ν π µν + π ανu ν , (15) where E = f 2 is the energy density seen by a comoving observer, E = u µ u ν T µν , and σ µν = ∆ αβ µν ∇ α u β with the traceless and spatial double projector ∆ αβ µν ≡ is the shear flow tensor. The shear stress tensor π µν is the traceless part of f µν 0 while the trace of the latter enters the isotropic pressure through P + Π = − 1 3 ∆ µν f µν 0 where Π is the bulk viscous pressure. The hydrostatic pressure is given by the equilibrium equation of state P = P(E, n) (where n is the conserved charge if applicable). Therefore, f µν 0 contains all nontrivial information about the dissipative pressure corrections Π and π µν . Their exact evolution is, therefore, described by the equatioṅ Differently from the four-momentum conservation equations (11) and (12), the equation above couples the components of T µν with moments of the distribution function that are not part of T µν or the particle current N ν = p p µ f ; this is seen in both the collisional kernel contribution on the left-hand side and the direct couplings on the right-hand side that survive even in the free-streaming limit C[f ] → 0. The most common approximation to obtain a closed set of equations is to assume that the deviations from local equilibrium are small and that the deformation of the distribution function is dominated by the nonequilibrium corrections in T µν and N µ [36]. In this way it is possible to obtain approximate expressions for f αµν −1 and f αβµν −2 that can be expressed entirely in terms of quantities usually associated with relativistic hydrodynamic behavior. This procedure leads to Israel-Stewart-like theories [36] of dissipative relativistic hydrodynamics, as discussed in detail in Ref. [20].
It is important to note, however, that thanks to Eq. (10) one can improve the hydrodynamic description by considering f αµν −1 and f αβµν −2 as dynamical variables in their own right. The dynamical equations for these moments will then couple to moments of different ranks and orders. This provides a way to systematically improve the solution of the moment equations that does not necessarily depend on the hypothesis of small gradients and approximate local equilibrium [20]. In Ref. [18] such an expansion was tested in 0 + 1 dimensions 1 and fast convergence was found to the corresponding exact solution 1 In fact, in [18] the local equilibrium expectation values were sub-of the RTA Boltzmann equation. The reasons behind this rapid convergence are not yet fully understood, especially in cases without a large degree of symmetry.
In the following section we point out that the method of moments becomes much more involved when longrange forces are added to the kinetic system via the coupling to mean fields. Nevertheless, rapid convergence can still be recovered by reorganizing the expansion in terms of the resummed moments presented in Sec. IV.

III. ADDING LONG-RANGE FORCES THROUGH MEAN FIELDS
Extending the results of the previous section to gas mixtures consisting of multiple particle species is conceptually straightforward [33]. New difficulties are encountered, however, when one introduces interactions with long-range mean fields, even for a single particle species. Such interactions are described by adding to the relativistic Boltzmann equation an additional force term, the Vlasov term, describing the momentum drift of the distribution function due to acceleration or deflection of particles [5].
The two most common types of such mean field interactions are medium dependent effective particle masses m(x) and gauge fields. Including them gives rise to the relativistic Boltzmann-Vlasov equation: Here ∂ µ p ≡ ∂/∂p µ is the partial derivative with respect to momentum. In this paper we only consider Abelian U (1) gauge fields with a single conserved charge (electromagnetism) though all mathematical considerations should also apply in the case of non-Abelian gauge fields [37][38][39]. For a recent derivation of dissipative relativistic magnetohydrodynamics from the relativistic Boltzmann equation coupled to Maxwell's equations using the method of moments [20] we refer the reader to Ref. [40].
Equation (17) assumes that all four components of p µ are independent, i.e. the particle momenta are in general off-shell, with the on-shell constraint imposed by the momentum integration measure in Eq. (7) when computing physical quantities. It is sometimes more intuitive to work directly with the on-shell distribution, which depends only on the spatial momenta p, with p 0 being replaced by m 2 +p 2 . As shown in Appendix C, the covariant equation (17) above can be rewritten as the following equation for the on-shell distribution function: However, in this approach one loses manifest relativistic covariance, which makes it difficult to easily switch between the global and fluid rest frames. Also, the derivation of the moment equations becomes more complicated. For these reasons we prefer to use Eq. (17) as our starting point. 2 Equation (17) must be solved together with the Maxwell equations that determine the electromagnetic field strength tensor generated by the moving electric charges in the gas: Clearly, an external electromagnetic field can also be added if desired (for example one generated by the charges of the colliding nuclei in relativistic heavy-ion collisions). Following the same steps as in the preceding section one can rewrite Eq. (17) as the following infinite set of coupled evolution equations for the f-moments: Here we used the following relativistic decomposition of the tensor field strength F µν into the electric and magnetic fields in the comoving frame: 3 with 2 In general, a position dependent mass requires the introduction of an additional mean field B µν for thermodynamic consistency, subject to certain constraints to ensure global four-momentum conservation and with its own equation of motion [41,42]. These extra equations only complicate the algebra but do not change conceptually the following considerations about the method of moments. Therefore we will neglect in the rest of the paper all technical details about the proper implementation of a medium dependent mass and focus mostly on gauge field interactions. 3 This is not the only convention found in the literature. The Levi-Civita symbol εµνρσ is a tensor density of rank 1, not a tensor [34]. The magnetic field defined in this way is not a vector but it transforms with an additional determinant of the Jacobian and gµν B ν = det(g)Bµ. Some authors multiply (divide for the upper case indices) the Levi-Civita by − det(g). In this way the magnetic field behaves like a tensor under orientation-preserving transformations. Although in Cartesian coordinates there is no difference, when considering curvilinear coordinate systems (such as Milne coordinates) this fact must be taken into account.
Compared with Eq. (10), Eq. (20) has five additional terms (the first five terms on the right-hand side) that describe couplings with the mean electromagnetic and mass fields. At first sight it seems straightforward to repeat the procedure described in the previous section to derive evolution equations for the hydrodynamic moments of the distribution function and to systematically improve them by adding the contribution from the nonhydrodynamic moments to which they couple dynamically. However, this simple extension of the moment expansion to the case of the BVM equations (17,19) faces additional difficulties, as we explain now.
In the Boltzmann case (1) the f-moments on the righthand side of Eq. (10) have the same physical dimensions as the ones on the left-hand side: except for the contribution from the collisional kernel, all terms on both sides of the equation have the same sum r+s of the energy and tensor rank indices r and s. In the Boltzmann-Vlasov case, on the other hand, neither the effective mass m nor the electromagnetic field F µν are dimensionless; therefore, they couple to f-moments of different physical dimensions. Specifically, the first five terms on the right-hand side of Eq. (20) involve f-moments whose index sum r + s is lower than that of the moment on the left-hand side. Either the energy index r or the tensorial rank s, or both, are reduced. Therefore the systematic improvement of the solution requires considering the influence of moments with ever-decreasing physical dimensions. We note that this already appears in the case of the Boltzmann equation for a single particle species using irreducible moments as shown in Ref. [20], where particles with constant mass m were considered. In this case, when m is nonzero there are terms that couple irreducible moments with others of reduced physical dimensions. However, in that case this coupling to moments of lower physical dimensions vanishes for m = 0.
In the case of BVM, even for massless particles, this coupling does not vanish because of the electromagnetic interactions. 4 In the massless limit these lowerdimensional moments become ill-defined for r+s < − 2, and even for nonzero but small masses m T their magnitude grows with increasingly negative values of r+s, which should affect the convergence of the moment expansion. This is easily seen from their definition (see Eqs. (4) and (9)) when writing out the momentum integral in LRF components, where the distribution function f is evaluated on-shell, p = |p|, p 0 = m 2 +p 2 , andp µ1 ≡ p µ1 /|p| is a spatial unit vector in the LRF which depends only on the momentum angles (θ p , φ p ), with dΩ = sin θ p dθ p dφ p . Scaling out the particle mass, p = my, this becomes where the angular integral ω µ1···µs ≡ dΩp µ1 · · ·p µs f is some finite dimensionless number. Eq. (24) shows that in the ultra-relativistic limit m/T → 0 the f-moments diverge for r+s < −2; for the equilibrium distribution f eq this is easily verified explicitly. In this limit, a systematic solution of the BVM equations requires taking into account the contributions from an infinite set of moments defined by the tower of moment equations in Eq. (20) as soon as the particle momenta are influenced by nonvanishing mean fields.
In principle, this problem can be addressed by solving Eq. (20) only for the moments with positive energy index, approximating the other (non-hydrodynamic) moments to which they couple by non-dynamic constitutive equations. This approach was pursued in Ref. [20] where irreducible tensor moments of negative r were expanded in terms of the corresponding irreducible tensors with positive r which are assumed to form a complete basis for the non-equilibrium correction to the distribution function (see Appendix D for an illustration).
We here propose a different method to avoid couplings to possibly ill-defined moments with r+s < −2, by introducing a new set of moments of the distribution function that resum an infinite number of f-moments. We show how all the physical information about the system (including both hydrodynamic and non-hydrodynamic moments of the distribution function) can be recovered from these resummed moments, and that their exact dynamical evolution is well-defined even in the ultra-relativistic m/T 1 limit. Furthermore, we show for a simplified physical situation that a solution in terms of an expansion in these resummed moments converges rapidly to the exact result from the BVM equations.

IV. RESUMMED MOMENTS AND THE HYDRODYNAMIC EXPANSION
Let us define and introduce the space-like projections The similarity of these definitions with Eqs. (4) and (9) is obvious: they differ only by the Gaussian weight factor e −ξ 2 (p·u) 2 under the integral. 5 By Taylor expanding this Gaussian 6 one sees that each φ-moment φ µ1···µs r can be written as an infinite sum of f-moments, with the same tensor rank s but different energy indices r. The dimensionful parameter ξ determines the relative weight of f-moments with different dimensions, i.e. with different factors of p · u, in the sum.
By construction, all resummed moments (26) of the same tensor rank s are related. 7 Using the relations one easily finds that The f-moments, in particular the hydrodynamic stressenergy tensor and charge current density, are recovered from the resummed moments of the same order and tensor rank via the relation Structurally, the exact evolution equations for the resummed moments are very similar to those for the fmoments in Eq. (20): Added complications are the extra ξ-dependent coupling terms in the last two lines, and the fact that the resummed moments depend on an additional continuous parameter ξ that effectively adds an extra dimension to the complexity of the numerical solution. This is the price we have to pay in our approach for avoiding the need in the standard approach [20] for expanding moments with negative energy index r in terms of positive energy index moments. However, due to the relations (29) between moments of the same tensor rank with different values of r, the moment equations can be written entirely in terms of moments with a fixed r, thereby removing the need to determine moments with lower r values that could become ill-defined. Hence, in this approach the only truly independent moments are those with different tensor ranks. Crucially, the moment equations (31) can be further simplified by using (29) to express everything through moments with r = 1 such that the terms proportional to r−1 vanish exactly: 5 The underlying idea is that this weight function "generates" inverse powers of p · u via ∞ 0 dξ e −ξ 2 (p·u) 2 = √ π/2 p · u . 6 There is more than one way to perform such an expansion. One can consider, for instance, −ξ 2 (p · u) 2 = x and use the Taylor expansion of e x ; or ξ 2 (p·u) 2 = x 2 and use the Taylor expansion of the Gaussian e −x 2 . All such series correspond to an infinite sum of f µ 1 ···µs r . This is why we call the Φ-and φ-moments "resummed Here the resummed moments under the integrals are understood as functions of the auxiliary variable √ υ instead of ξ. These equations are well-defined for any tensor rank s, even in the massless limit. The expansion in terms of resummed moments whose dynamics follows Eq. (32) therefore provides a well-defined and systematically improvable generalization of the hydrodynamic expansion of the Boltzmann-Vlasov equation in the presence of mean field forces.
It is worth noting that after integrating Eqs. (32) over ξ one recovers Eq. (20) for the f-moments with energy index r = 1 exactly. More generally, with appropriate manipulations using relations (29) one recovers the evolution equations for all the well-defined f-moments, before any approximations (in particular the lowest-order ones that form the equations of hydrodynamics in the presence of mean fields). 8 However, by using the resummed φ-moments instead of the f-moments we avoid the numerical issues associated with f-moments with sufficiently negative energy index.
For the simultaneous solution of the Maxwell equations (19) one expresses the electric current in terms of the resummed moments as follows (see Eq. (30)): To test numerically the convergence properties of the set of equations (32) we study a situation where one can solve the BVM equations exactly. The approximations for physically meaningful macroscopic quantities generated by truncating the moment expansion can then be checked order by order.
We will focus on the massless case where the complications mentioned in the previous section are most relevant. In the absence of mean-field forces an exact solution of the Boltzmann equation in RTA was found in Refs. [23][24][25] for a transversally homogeneous gas undergoing longitudinally boost-invariant expansion, i.e. Bjorken flow [31]. For this case approximations made within the traditional moment expansion [20] for the dynamics of the shear stress tensor are known to provide a very good description of the exact result obtained from the Boltzmann equation [43], for reasonable values of transport moments" of the distribution function. 7 The Gaussian weight e −ξ 2 (p·u) 2 falls off quickly enough to preserve uniform convergence in all the following manipulations. 8 By applying consistent approximations to the higher-order, nondynamical moments (see discussion below) one can ensure that the transport coefficients are the same in both approaches.
coefficients. We now show that this exact solution can be extended to the RTA Boltzmann-Vlasov case, which will be used as our testing ground.

A. Exact solution of the BVM equations for Bjorken symmetry
Because of the symmetries of Bjorken flow it is convenient to formulate the problem in Milne coordinates The quantities τ and η are called, respectively, longitudinal proper time and space-time rapidity or, in short, proper time and rapidity. Correspondingly, we use for the four-momentum For more information about the underlying symmetries of the Bjorken expanding fluid in hydrodynamics and kinetic theory see Refs. [44] and [26,27], respectively.
The symmetry of the expansion imposes that the (onshell) distribution functions of all particle species depend only on the proper time, the transverse momentum p T , and the longitudinal momentum p η . If the electromagnetic field is dynamical and there are no external sources, Bjorken symmetry also constrains the fields and the electric current. The simplest case consistent with all symmetry constraints is an overall charge neutral twocomponent gas of particles f (τ, p T , p η ) and antiparticles f (τ, p T , p η ), expanding in a purely longitudinal electric field E η (τ ) without magnetic components, and evolving according to the following equations for the on-shell distributions and electric field: Here τ r is the relaxation time present in the RTA collision term. To preserve the conformal symmetry of this massless plasma we assume τ r (τ ) = c/T (τ ) where T (τ ) is the effective temperature of the system at proper time τ and the constant c is related to the specific shear viscosityη ≡ η/s (i.e. the ratio of shear viscosity to entropy density) of the system by c = 5η [26,27,35,45]. The equilibrium distribution function reads, for both particles and antiparticles, For a system undergoing Bjorken expansion the flow fourvelocity reduces to u µ = (1, 0) in Milne coordinates. The effective temperature T is defined through the usual Landau matching prescription E = E eq = 48πkT 4 where k ≡ g/(2π) 3 (see Eq. (7)). Appendix E presents an explicit derivation of the evolution equations and the demonstration that the symmetries ensure that in Milne coordinates the flow remains static at all proper times.
The system of equations (35)-(37) is very similar to the one studied in Ref. [46], except that we here consider a U (1) electromagnetic gauge field (which is not related to an Abelian subgroup of color SU (3)) and also neglect the effects from nonperturbative electric field decay, which are induced via the Schwinger mechanism. The latter is exponentially suppressed by the inverse electromagnetic coupling α −1 , and thus in heavy-ion collisions the electromagnetic Schwinger effect acts on a much slower time scale than those characterizing the strong interactions.
The solution of Eqs. (35) and (36) for the particle and antiparticle distributions reads The initial conditions are f = f 0 andf =f 0 at τ = τ 0 . The damping function D and the momentum shift ∆p η are defined as The implicit equations (39) and (40) can be solved, for arbitrary initial distributions f 0 andf 0 , by numerically iterating the temperature and electric field profiles, T (τ ) and E η (τ ), subject to the Landau matching constraint and the solution to the Maxwell equation (37) following the method used in Ref. [46]: One inserts an initial guess for the temperature and electric field profiles into the right-hand sides of Eqs. (39) and (40) and then uses the distribution functions obtained from these equations to obtain new profiles from Eqs. (43) and (44). The process is iterated until convergence of the profiles, at some desired precision, is achieved. The resulting particle and antiparticle distribution functions form (for all practical purposes) an exact solution of the coupled BVM equations for systems with Bjorken symmetry. As a basis for comparisons with the resummed expansion discussed further below, we consider a set of thermal equilibrium initial conditions with parameters that yield energy densities that can be considered reasonable for problems of interest in the field, such as high energy proton-antiproton collisions. We take for the initial temperature T 0 = 0.3 GeV defined at the initial time τ 0 = 1 fm/c. For the initial electric field we assume a temperature-normalized ratio E 0 η /T 0 = (τ 0 E 0 L )/T 0 = 1/5, corresponding to a normalized value for the longitudinal Cartesian component of the initial electric field given by E 0 L /T 2 0 ≈ 2/15. In dimensionful units this initial electric field has the value E 0 L ≈ 0.3 fm −2 . This is a natural order of magnitude in this case because it corresponds to the electric field of a parallel plate capacitor made of conducting sheets representing the central collision between a proton and an antiproton that are both (infinitely) Lorentz contracted, with charge densities σ corresponding to the charge of a proton smeared over a disk of radius 1 fm: Figures 1a-d show the evolution of the temperature T /T 0 , longitudinal/transverse pressure ratio P L /P T , 9 normalized electric current J L /T 3 , and electric field E L /E 0 L , respectively, for four different choices of the specific shear viscosity, 4πη = 1, 3, 10, and 100, covering the range from the KSS result [47] to almost free-streaming. We compare the evolution according to the BVM equations (solid lines) with that of a system of uncharged particles, i.e. without electric field (dash-dotted lines).
According to Fig. 1a the temperature evolution is not significantly affected by the presence of the electric field. This is easily explained by the fact that the energy density is dominated by the particles, with an almost negligible field energy density contribution of E 2 /2 = E 2 L /2 = E 2 η /(2τ ) which, at τ = τ 0 , corresponds to 0.7% of the particles' energy density 48πkT 4 . On the other hand, the pressure anisotropy P L /P T shown in Fig. 1b exhibits a much stronger sensitivity to the interaction of the charged particles with the electric field. This interaction accelerates the isotropization of the pressure, most prominently for intermediate values of the specific shear viscosityη.
Compared to the study performed in Ref. [46] where a similar set of equations was investigated, the effect of the electric field on the pressure anisotropy is significantly reduced in our study. This is due in part to the smaller electric fields considered here, but the lack of a Schwinger term for the spontaneous decay of the fields in our treatment also contributes to this difference. As already stated we neglect the latter (and the emission of photons by the accelerated charges) since their rates are suppressed relative to strong-interaction collisions by the smallness of the electromagnetic coupling constant α = 1/137.
Proceeding to the electric current shown in Fig. 1c one observes at early times that the normalized diffusion current J L /T 3 increases with increasingη (i.e. for larger collisional mean free path), but that for smallerη (i.e. larger microscopic collision rate) it persists longer. For 4πη = 100, i.e. close to the free streaming limit, the electric current is large enough to eventually flip the sign of the electric field, as seen in Fig. 1d. For the initial conditions chosen here this happens about 7 fm/c after starting the simulation. This reversal of the direction of the electric field repeats at larger times, i.e. over long time scales the electric field oscillates in time. This will be further explored in Fig. 4 below.
For later comparison with the resummed moment expansion we also studied the exact solution for more ex-treme cases, corresponding to much larger initial electric fields (×20) and larger initial expansion rates (×4). They are shown in Fig. 2. To increase the initial expansion rate θ 0 = 1/τ 0 we simply decrease τ 0 from 1 to 0.25 fm/c. This entails a larger initial Knudsen number Kn 0 ≡ τ 0 r θ 0 = τ 0 r /τ 0 which increases from Kn 0 < ∼ 1/3 to Kn 0 > ∼ 1 as we decrease τ 0 from 1 to 0.25 fm/c. Typically, non-hydrodynamic behavior is expected to occur when the Knudsen numbers are not sufficiently small. Figure 2a shows that changing the Knudsen number affects the temperature evolution only quantitatively (the system cools faster as the expansion rate and Knudsen number increase) but not qualitatively. When increasing the initial electric field by a factor 20, however, the temperature starts to evolve non-monotonically: the initial cooling stage is soon (after about 1.5-2 fm/c) followed by an extended reheating period where Ohmic heating (i.e. the dissipation of the initially large electric field due to collisions among the charged particles) causes the temperature to increase again, in spite of the longitudinal expansion.
The evolution of the pressure anisotropy shown in Fig. 2b exhibits even stronger sensitivity to the initial expansion rate and, in particular, the initial electric field. The reheating seen in Fig. 2a for the cases with large initial electric fields is accompanied (in fact, preceded) by a reversal of the pressure anisotropy from P L /P T < 1 to P L /P T > 1. Fig. 2c shows that the normalized electric current J L /T 3 increases with both the initial electric field and initial expansion rate. Without the normalization by T 3 , the non-monotonic behavior of T (τ ) seen in panel (a) actually causes J L to oscillate, too, for large initial electric fields (not shown in Fig. 2). Only the electric field, shown in Fig. 2d, decreases continuously, without exhibiting any non-monotonic effects on the same time scales. Reheating, reversal of the pressure anisotropy and oscillations in the magnitude of the electric current do not occur for the smaller of the two initial values for the electric field studied here, irrespective of the initial expansion rate.
A careful comparison between Figs. 1d and 2d leads to an interesting observation: While increasing the specific shear viscosity 4πη at fixed τ 0 (i.e. fixed initial expansion rate) (Fig. 1) and increasing the initial expansion rate (by reducing τ 0 ) at fixed 4πη (Fig. 2) both increase the initial value of the Knudsen number Kn 0 , they have opposite effects on the initial rate of decrease of the electric field. Fig. 1d shows that E L initially decreases more rapidly when 4πη is increased at fixed τ 0 whereas we see in Fig. 2d that E L initially decreases more slowly when the initial expansion rate is increased at fixed 4πη, for both choices of the initial electric field value. Furthermore, we note from Fig. 2d that, at fixed initial expansion rate, the electric field initially decreases relatively more slowly, but at later times relatively more quickly as the initial value E 0 L /T 0 is increased. These observations reflect an interesting interplay between the two transport mechanisms related to the shear stress and charge diffusion currents that control the dissipative effects in this theory.
We point out that for the larger of the two initial electric field values studied here the system is always far away from the Navier-Stokes limit. In order to see the ultimate approach to local thermal equilibrium we show in Fig. 3 the long-time behavior of the macroscopic quantities studied in Figs. 1 and 2, for the most extreme case with both large initial electric field E 0 L /T 0 = 4 fm −1 and large initial expansion rate θ 0 = 4 fm −1 c. Clearly the system remains far from equilibrium up to 50 fm/c when the pressure ratio P L /P T finally returns to values close to 1. Beyond that time, the temperature, electric field, and electric current approach their asymptotic limits as power laws. 10 In this regard, it would be interesting to study how the electric field affects the existence of hydro- dynamic attractor solutions [48][49][50][51][52] for the BVM equations. We already mentioned in the discussion of Fig. 1 that, for very large mean free paths, 4πη → ∞, corresponding to the collisionless Vlasov limit, the electric field flips sign at rather early times, hinting at long-time oscillations. This is further explored in Fig. 4 where we consider 4πη = 100 and show the behavior of the solid red lines in Fig. 1 over a larger time interval τ − τ 0 = 100 fm/c, together with a similar case with 20× larger initial electric field and 4× larger initial expansion rate. One sees that for large collisional relaxation times both the electric current and electric field (panels (c) and (d)) exhibit prominent long-time oscillations, even in the absence of a Schwinger term describing spontaneous electric field decay [46]. For the moderate initial value of the electric field selected, E 0 L /T 0 = 0.2 fm −1 (see red solid lines), these field oscillations have only a minor effect on the evolution of the pressure anisotropy (shown in panel (b)) and leave no visible trace in the evolution of the temperature (panel (a)). This is, as before, easily explained by the small overall contribution of the electric field to the total energy density which is strongly dominated by the particle contribution. For much larger initial electric fields (grey dash-dotted lines), however, the oscillating energy content stored in the electric field and current is note that the current and electric field approach their asymptotic power laws sooner than the temperature. large enough to cause, through Ohmic heating, oscillations even in the temperature (a) and pressure ratio (b).

B. Resummed moment equations for BVM with Bjorken expansion
Since our plasma contains both particles and antiparticles we must consider two generations of resummed moments, φ µ1···µs r for the particles andφ µ1···µs r for the antiparticles. The symmetry of the expansion reduces the number of independent degrees of freedom: Because of homogeneity in the transverse plane, moments with an odd number of x or y indices vanish -only pairs (or, more generally, even numbers) of x and y indices yield nonzero moments. In particular, With this result we can express all of the non-vanishing moments with x and y indices through moments of equal or lower rank with only η indices. At this point it is convenient to use the longitudinal vector z µ with Cartesian components (sinh η, 0, 0, cosh η)  Fig. 1, but evolved over longer times. Grey dash-dotted lines: Similar, but for 20× larger initial electric field and 4× larger initial expansion rate.
or Milne components (0, 0, 0, 1/τ ) as it allows us to express the factor p η /τ covariantly as p η /τ = p · z. We use it to introduce the independent scalar moments For the Bjorken expanding case, the scalars {φ ± 0 , · · · , φ ± l } contain all the information about the resummed moments up to tensor rank l of the system. Their exact evolution can be obtained directly from Eq. (31), taking l projections along z and summing (subtracting) the particle and anti-particle equations. An equivalent method consists in taking directly the τ derivative of the right-hand side of Eq. (46). Either way one finds In the last term of Eq. (47) it is understood that for l = 0 the first term in the square brackets is zero. The equilibrium moments φ − l,eq vanish for all l while the equilibrium +-moments φ + l,eq vanish for odd l = 2n+1. For even l = 2n one has φ + 2n,eq = 8πk T (2n+2)(2n)! (2ξ) 2n+4 U 2+n, where U (a, b, z) is the Tricomi confluent hypergeometric function [53]. Full knowledge of the hydrodynamic moments (i.e. the stress-energy tensor and the electric current) is encoded in the three scalar moments φ + 0 , φ − 1 , and φ + 2 . Therefore, at leading order in the resummed moment expansion only the φ ± l up to l = 2 are considered as dynamical variables, i.e. we truncate the hierarchy of moment equations (47) at l = 2. As higher order corrections we will progressively consider the higher order moments as additional dynamical variables. 11 When truncating the moment hierarchy at some maximum l value l max , there is no unique prescription to approximate the moments of orders l max +1 and l max +2 appearing on the r.h.s. of Eq. (47). The simplest one is to assume φ ± l ≈ φ ± n,eq for n > l max . For numerical purposes this approximation is very helpful because we can use the simple formula (48) to evaluate the non-dynamical higher-order moments. It must be noted, though, that 11 Because of symmetry, the odd φ + 2n+1 and the even φ − 2n never couple directly to any of the other moments. Their evolution couples with the hydrodynamic moments only indirectly through the electric field, which is itself coupled to the electric current. In particular, if the initial values of all these moments are zero (for instance for local equilibrium initial conditions) their evolution is trivial. this approximation is not the most accurate possible.
Recall that the tensor moments φ µ1···µs r are defined on a non-orthonormal polynomial basis of tensors constructed from the locally spatial momentum vectors p µ . Due to the lack of orthogonality, some information on the higherorder moments is already contained in the lower order ones and this could be exploited. 12 We will here ignore this extra information, leaving the most effective approximation scheme for the non-dynamical higher-order moments for future research.
It is numerically convenient to normalize the scalar projections φ ± l in order to confront pure numbers of similar order of magnitude in the evolution equations. For even l = 2n we do so by dividing them by the ξ → 0 limit of their equilibrium values: The right hand side can be extended to the case of odd l, and we can use the same normalizing factors for both plus and minus moments: 12 The situation is analogous the one presented in [18] for the (nonresummed) reducible moments.
Obtaining from Eq. (47) the evolution of the normalized moments (50) is straightforward: It must be noted that for l = 0 the (otherwise undefined) first term ∼ M ∓ −1 inside the square brackets is absent. The effective temperature is defined through the Landau matching procedure (43). Making use of the normalized moments, the effective temperature can be written as 52) or, equivalently, Implementing Landau matching is numerically expensive since the temperature must be matched to the energy density moment by evaluating the integral on the r.h.s. of (53) at each time step. It turns out to be more convenient to instead treat T as an additional dynamical variable. Its exact evolution can be obtained by differentiating the r.h.s. of Eq. (53) with respect to time and inserting Eq. (51) for l = 0. One finds The last equation to solve numerically is the evolution of the electric field (44). Written in terms of the normalized moments it becomes These last two equations, coupled to Eqs. (51) for those moments that we treat dynamically, can be solved with an explicit Runge-Kutta method. In the following subsection we compare the evolution of the macroscopic quantities (i.e. the temperature, pressure anisotropy, electric current, and electric field) between the solution of the resummed moment equations and the exact solution of the Boltzmann-Vlasov-Maxwell equations discussed in the previous subsection.

C. Comparison of the resummed moment expansion with the exact solution of the BVM equations
As an example for studying the convergence of the resummed hydrodynamic moment expansion, we consider the default initial conditions used for the exact solutions shown in Fig. 1 (T 0 = 0.3 GeV and E 0 L /T 0 = 0 or 0.2 fm −1 at τ 0 = 1 fm/c), with minimal shear viscosity 4πη = 1 (black lines). In Fig. 5 we compare these exact solutions with three iterations of the resummed moment expansion, corresponding to truncations at different l max values as described in the preceding section, i.e. with thermal equilibrium values M ± n,eq for all moments of order n > l max . In Fig. 5 the results from the resummed moment expansion are nearly indistinguishable from the exact results already for l max = 6. Fig. 6 shows that a relative precision better than 1/1000 is achieved with l max ≥ 8. We verified the continued convergence of the series up to l max = 34 where we stopped the calculation. Figure 6 shows that for τ −τ 0 < 10 fm/c the resummed moment expansion agrees with the exact solution to better than 10% already at leading order l max = 2, i.e. in the hydrodynamic limit. The largest deviation is seen for the current J L /T 3 where it reaches about 8%, followed by 3.5% for the pressure ratio and 2.5% for the electric field. We found similar levels of (im)precision for all other parameter sets studied in this work. Larger deviations from the exact result at leading order are accompanied by slower convergence towards the exact result as the order of the approximation is increased. But in most cases l max = 12 resulted in relative deviations from the exact result of less than 1%. Large Knudsen numbers (large expansion rates) and large initial pressure anisotropies have surprisingly weak effects on the speed of convergence.
A significant slowing of the rate of convergence was noted only for very large initial electric fields and/or large relaxation times (i.e. large specific shear viscosities, close to the free-streaming limit). In both cases the reason for slower convergence is easily understood by inspecting Eq. (51). There are only two couplings to moments of higher orders, one of which is multiplied by the electric field. No matter the order of the approximation, in the equation for l = l max the moments M ∓ lmax+1 and M ± lmax+2 must be approximated. A large electric field, or better yet, a large (dimensionless) ratio E η (τ )/T 0 , magnifies the coupling to one of the approximated moments. The approximation error then affects the evolution of all the lower moments since they all couple to the higher ones via the electric field. For large relaxation times one must recall that the temperature decreases more slowly than the expansion rate, due to Ohmic heating. Normally the second term on the left hand side of Eq. (51), which is proportional to 1/τ r = T (τ )/(5η), soon dominates the evolution since the scalar expansion rate θ = 1/τ suppresses all terms on the right hand side. This stage is, however, postponed asη → ∞. In the free-streaming limit the right hand side, including the terms coupling to the higher order moments, becomes dominant, and thus the evolution becomes more sensitive to approximations made for the moments of order n > l max .
While illuminating, these considerations are not enough to understand the reasons for the observed generically rapid convergence of the resummed moment expansion, nor do they throw light on possible situations that might spoil this convergence when going from the simple case of Bjorken flow to a more realistic 3-dimensional expansion. It is worth noting that nowhere in the discussion of the convergence of the resummed moment expansion is the influence of large gradients or large pressure corrections mentioned. The present formalism thus appears unrelated to the traditional expansion in powers of Knudsen and inverse Reynolds numbers [20] which have been used to define the effective range of validity of relativistic hydrodynamics. Within the resummed moment expansion, the only approximation lies in the coupling to moments of order n > l max , and only these (non-hydrodynamic) moments are approximated (in our case with the assumption of small deviations from equilibrium 13 ).
So how does approximating these higher-order, nonhydrodynamic moments affect the dynamical evolution of the hydrodynamic moments, i.e. of the charge current and energy-momentum tensor? The following observations may bring us closer to an understanding of the good convergence of the moment expansion even in situations where traditional arguments based on the size of hydrodynamic gradients and/or dissipative flows (i.e. on Knudsen and/or inverse Reynolds numbers) suggest the breakdown of hydrodynamics. Looking at Eq. (51) one sees that the higher order moments always appear multiplied by a factor (ξT 0 ) 2 . A crucial observation is that, in the situations studied in this paper, the resummed moments are all strongly peaked around ξ 0, especially the higher order ones (see Fig. 7). This is because the resummed moments in Eq. (25), and hence their (normalized) (0 + 1)-dimensional versions in Eqs. (46) and (50), contain a Gaussian factor e −ξ 2 (p·u) 2 that effectively suppresses all integrals at large values of ξ. It is especially true for the higher-rank moments which, due to extra momentum factors in the integrand, probe the distribution function at higher momenta where the suppression by the Gaussian becomes effective already at lower values of ξ.
On the r.h.s. of Eq. (51) we already grouped the terms (both those with and without an electric field factor) accordingly. In both lines the term with the lower-order dynamically evolved moment will dominate over the one involving the higher-order moment because, in the ξ region where the moments are substantial, the latter is suppressed by a factor ξ 2 . For l = l max , in particular, this means that the coupling to the moments of order l max +1 and l max +2 that are not evolved dynamically, but instead approximated by their equilibrium values, are suppressed. This suppression of the influence from higher-order moments repeats at each step of increasing the order l, which largely explains the observed rapid convergence of the resummed moment method.  51), plotted as functions of the transformed variable t ≡ (ξT 0 −1)/(ξT 0 +1) whose origin lies at ξ = 1/T 0 . The blue-dashed lines show the subleading second terms in these two lines (the ones that describe the coupling to moments of higher order) as obtained from a calculation where all moments up to order l max = 24 were evolved dynamically. The red-dotted lines show the effect of truncating the resummed moment expansion at leading order (i.e. at l max = 2) and approximating the moments M + 4 and M − 3 by their thermal equilibrium values. In both panels the coupling terms to the higher-order moments are seen to be dwarfed by the leading terms over the entire ξ range, but especially at ξ = 0 (t = −1) where the resummed moments peak.
The suppression of coupling terms to higher order resummed moments is expected to partially persists in 3+1 dimensions. It is difficult to think of a realistic phase-space distribution that will not produce resummed moments that peak near ξ = 0 and decay for ξ larger than the inverse of the typical energy scale of the system (which is normally larger than the temperature to which the system settles after thermalization). For massive particles a momentum-independent suppression factor e −ξ 2 m 2 can be extracted from the defining integral (25). It is worth noting that the Gaussian suppression factors arise only in the resummed expansion; the ordinary expansion equations for the hydrodynamic moments are recovered from the resummed moments by performing a ξ integral for the scalar, vector and rank two tensor equations.
Still, looking at the general equation (32) one notices two terms in the second integral on the r.h.s. that couple to moments of higher tensor rank without being suppressed by factors of ξ 2 : the first one involving the acceleration and the last one which is a space-divergence. In the general case these terms have the potential of spoiling the rapid convergence of the resummed moment expansion and thus the good agreement of the leadingorder hydrodynamic approximation with the exact solution of the Boltzmann-Vlasov-Maxwell equations. Future research will settle this issue. There is, however, a somewhat handwaving argument suggesting that the re-summed moment expansion will indeed continue to converge in the general (3+1)-dimensional case. Looking again at Eq. (32) one sees that at any order only two higher order moments couple without ξ 2 suppression to the evolution of the lower ones (five if one includes the ξ 2suppressed couplings). On the other hand, the number of moments of the same or lower rank increases rapidly with the tensor rank s. Arguing that all the moments of a given tensor rank and order have the same physical dimensions and are therefore expected to be of similar order of magnitude, exact cancellations between many of them would be needed for their contribution to the evolution equations to be overwhelmed by the two terms coupling to higher ranks or orders. This argument works both with and without coupling to the electromagnetic fields.
We note that a similar argument cannot be made for the traditional moment expansion, even in the absence of electromagnetic fields where the coupled equations couple only to moments that remain well defined and well behaved in the limit m/T 1. Equation (10) for the evolution of the standard f-moments contains s+3 terms coupling to moments of different order and only s+1 moments of the same tensor rank s and energy index r. Moreover, two of the couplings to different moments are multiplied by r and r−1. Recalling that all members of the tower of moments coupled to the hydrodynamic moments feature the same sum r+s, one must therefore expect that at the truncation step the evolution of the last dynamically evolved moments is dominated by the r + r − 1 + 1 + s approximated moments, unless the approximation is fine-tuned to render their combined contribution subleading relative to the s+1 dynamically evolved moments. Such a cancellation does not even happen automatically in the highly symmetric Bjorken case; therefore, even there one must put some effort into designing a good enough approximation scheme for the nondynamical moments that does not ruin the convergence properties of the expansion.
We close this section with the following remark. In this work we used the approximation M ± n = M ± n,eq for the normalized moments of order n > l max where l max is the order of the last dynamically evolved moment. This corresponds in the general (3+1)-dimensional case to the approximation φ µ1···µs r = φ µ1···µs r | eq. . Due to its simplicity, this is a particularly convenient approximation if one wants to check the convergence properties up to very high orders. However, it does not make full use of all the information from the lower-order moments. This is illustrated in Fig. 9 where we compare the precision of the resummed moment method in reproducing the exact pressure anisotropy obtained from the Boltzmann equation (without electric fields) at different orders of truncation with that of two well-established second-order dissipative hydrodynamic approximations for uncharged fluids that were derived using different methods, namely the Müller-Israel-Stewart (MIS) [36,54] and Denicol-Niemi-Molnár-Rischke (DNMR) theories [20]. Fig. 9 shows that at leading order the resummed moment approximation reaches almost the same level of precision as the MIS theory, albeit with deviations from the exact result that have opposite signs. The DNMR approximation, on the other hand, is significantly closer to the exact solution; with the resummed moment expansion, using our approximation scheme, one must go to l max = 6 (i.e. NNLO) to reach comparable precision. At leading order without electric fields, the moment evolution equations for the resummed moment expansion and the traditional method of moments (DNMR) have the same form; they differ only in the approximation for the higher order moments. In particular, the DNMR approach does not assume thermal equilibrium values for the higher-order moments, but instead relates them to the dynamically evolved lower-order moments. One can imagine doing something similar in the resummed moment expansion. Other approximations can also be considered, such as anisotropic hydrodynamic [43] expectation values. We leave the problem of optimizing the approximation for the non-dynamical higher-order moments to future research.

VI. CONCLUSIONS
In this paper we studied the emergence of hydrodynamics from kinetic theory for an overall electrically neutral system of charged particles described by the Boltzmann-Vlasov-Maxwell equations. We found that for this problem a straightforward generalization of the standard method of moments becomes considerably more intricate, especially in the massless limit. We pursued a different approach which consists of introducing an expansion in terms of a set of resummed moments that converges rapidly for systems both with and without long-range mean fields. As was the case for the standard moments, the components of the particle (charge) current and energy-momentum tensor (i.e. the hydrodynamic degrees of freedom of the system) can be expressed through certain low-order resummed moments. Truncating the hierarchy of coupled evolution equations for the resummed moments at lowest order yields a set of hydrodynamic equations. The precision of the description of the microscopic dynamics via moment equations can be systematically improved by truncating the hierarchy at higher order, thereby including the effects of nonhydrodynamic moments of the distribution function into the evolution of the hydrodynamic degrees of freedom. Therefore, the resummed moments method provides a convenient and reliable basis for the hydrodynamic expansion of a plasma interacting with long-range mean fields.
We studied the convergence of this expansion, and the precision of the resulting macroscopic hydrodynamic description as a function of the truncation order, for a highly symmetric expansion scenario (Bjorken flow), where symmetries (longitudinal boost invariance and transverse homogeneity) simplify the system of equations to the extent that the microscopic kinetic evolution defined by the BVM equations can be solved exactly. Exploring a range of model parameters such as the initial expansion rate (controlled by the starting time τ 0 ), the initial value of the electric field, and the specific shear viscosity of the fluid (controlled by the relaxation time encoded in the RTA Boltzmann collision term), we found uniformly rapid convergence of the moment expansion to the exact solution.
As a figure of merit one should remember that a precision of better than 1/1000 for the hydrodynamic variables is typically reached by truncating the resummed moment expansion at order l max = 8 or higher. When performing the truncation, we here approximated the higher-order moments to which the last dynamically evolved moment couples by their thermal equilibrium values; though other approximation schemes are possible they were not investigated in this paper. The most efficient way to slow down the convergence of the expansion is to increase the viscosity to large values, i.e. to push the theory towards the free-streaming limit. Slower than typical convergence was also found for very large values of the initial electric field. 14 The first of these two observations can be explained by observing that in Eq. (51) the coupling of the moment M ± l to itself dominates over the coupling to the (approximated) higher moments as long as the inverse Knudsen number τ /τ r is small, but that the evolution becomes very sensitive to the approximation of the higher moments in the opposite limit, i.e. for large relaxation times τ r . The second observation is explained by the coupling of M ± l to M ∓ l+1 (which at truncation level is approximated non-dynamically) and the fact that this coupling is magnified by a large electric field.
In Section V C we explored the reasons for the observed fast convergence of the resummed moment expansion. We found that, differently from the standard method of moments, the Gaussian factor in the integrand of the definition of the resummed moments suppresses most couplings to higher-order moments, and thereby reduces sensitivity to the approximations made at the truncation level for the non-dynamical higher-order moments. It should be noted, however, that this argument fully holds only in a Bjorken expansion. A weaker argument is expected to apply also for the general (3+1)-dimensional case, but this needs to be validated by future numerical simulations.
On the other hand Therefore the projection of Eq. (A1) onto u α reads which simplifies to i.e. an exact equation for the divergence of F µ1···µ r . In particular, for = 1 and r = 0 the contraction of Eq. (8) with the four velocity corresponds to the conservation of energy and momentum: The same arguments can be used to extend the above result to the Boltzmann-Vlasov equation and to the case of the resummed moments Φ µ1···µs r defined in (25).
this particular case the relative error therefore reads (D4) In Fig. 10 we show the relative precision of the approximation for different truncation orders N (2) . It is necessary to consider 12 irreducible moments (N (2) = 12) of positive energy index to reproduce f zz −2 with 10 −2 accuracy (relative error smaller than 5%). One needs N (2) = 46 to reach 10 −3 precision, and N (2) = 152 for 10 −4 .
In the particular case of interaction with an electromagnetic field, the term proportional to f zz −2 can easily become the dominant one in the evolution of the electric current, and this indeed happens in the cases we tested numerically in this paper. If the dominant contribution on the r.h.s. of Eq. (20) is not reproduced at the desired precision, it is very unlikely that the desired precision can still be achieved for dynamically evolved moment on the l.h.s. This example shows that the treatment of the moments with negative energy index proposed in [20] can become numerically costly, because of the large number of degrees of freedom that must be taken into account. 15 Moments with higher tensor rank have the same problems. While this polynomial series has been very helpful for the lowest order of the expansion needed to obtain Israel-Stewart-like viscous hydrodynamic theories, and for situations very close to local thermal equilibrium, it is not convenient for the analysis done in this paper.
As a final remark we note that making use of the resummed moments presented in Sect. IV corresponds to having no truncation in Eq. (D4) at all, in the sense of taking N (l) → ∞. The numerical precision of the code, however, limits the accessible information on the f µ1···µs r moments for very large r values since those are represented by many consecutive derivatives of the resummed φ µ1···µs 1 moments. 15 In [18], on the other hand, a precision of the order of 10 −3 in the final solutions is obtained already with a handful of degrees of freedom, by using a modified set of moments. In this appendix we will consider the equation for the on-shell distribution functions of a multi-particle gas in Milne coordinates. We will show that it is necessary to have at least a two component system (particles and antiparticles) in order to have a longitudinally boost invariant and transverse homogeneous expansion and, at the same time, fulfill the Maxwell equations.
According to Eq. (18), making use of the simplifications between the connection coefficients terms, one finds in Milne coordinates (just like in the Cartesian ones) p · ∂f + qF iα p α ∂f ∂p i = −C.
We omitted the "on" in f on and the hat inp because in this section we will consider only the on-shell version of the Boltzmann-Vlasov equation. If one assumes that particles have a longitudinally boost invariant and transversely homogeneous distribution one has f = f (τ, p T , p η ). In particular, in order to have the so called Bjorken expansion, one must add the requirement f (τ, p T , p η ) = f (τ, p T , −p η ), i.e. Z 2 symmetry. Therefore, ∂f ∂f ∂p y =