Lattice Boltzmann method for warm fluid simulations of plasma wakefield acceleration

A comprehensive characterization of lattice Boltzmann (LB) schemes to perform warm fluid numerical simulations of particle wakefield acceleration (PWFA) processes is discussed in this paper. The LB schemes we develop hinge on the moment matching procedure, allowing the fluid description of a warm relativistic plasma wake generated by a driver pulse propagating in a neutral plasma. We focus on fluid models equations resulting from two popular closure assumptions of the relativistic kinetic equations, i.e., the local equilibrium and the warm plasma closure assumptions. The developed LB schemes can thus be used to disclose insights on the quantitative differences between the two closure approaches in the dynamics of PWFA processes. Comparisons between the proposed schemes and available analytical results are extensively addressed.


Introduction
The process of particle acceleration plays a role of primary importance at the interface between fundamental and applied physics [1].The growing costs and dimensions of conventional large-scale accelerators demand for new and more efficient technologies to push the energies reached by particle beams beyond the state of the art of modern day capabilities.In this context, plasma acceleration is a promising technique that would enable the construction of compact particle accelerators, while retaining the same (or superior) energy gains obtained with conventional methods [2][3][4].A ionized gas is perturbed via the injection of relativistic charged particles (particle wakefield acceleration, PWFA) [5] or intense lasers (laser wakefield acceleration, LWFA) [6], generally named driver: the interaction between the neutral plasma and the injected driver creates a wave like dynamics of positive and negative charges and hence strong accelerating fields (up to 100 GV/m) are developed; the interested reader might look into [7,8] to go into detail on the topic.The described process involves a large number of "actors": the injected driver components, whether particles moving near the speed of light -typically electrons -or laser fields, and both the ions and electrons that make up the plasma.All of them interact with each other via electromagnetic forces, thus making it really difficult to predict and control the final behavior of a plasma acceleration experiment.Theoretical modelling and numerical simulations are therefore a powerful tool to help guide the design of new experiments.On the side of numerical simulations, the most commonly used techniques in the field are represented by particle in cell (PIC) methods [9][10][11][12][13], that employ single particle dynamics to describe both the particles in the driver (in this paper we will focus on PWFA) and the plasma components.These schemes are deeply fine-grained and, while they can capture the most refined phenomena that happen at the microscopic scale, they bring in low-statistics numerical noise [9,14].An alternative numerical modelling can be proposed by recurring to continuum fluid descriptions of the system.Fluid models describe the plasma via macroscopic fields such as particle number density and fluid velocity, and they do so by solving the inviscid relativistic Euler equations that can be systematically derived from kinetic theory of charged gases via a coarse graining of the relativistic Maxwell-Vlasov system [15][16][17] (in their most general warm formulation these equations are not closed and hence require additional constrains -more on this later on).Despite losing the ability to describe some kinetically pertinent features, numerical methods relying on the fluid description are still able to capture non trivial features of the PWFA system, and are set up by construction not to show statistical noise.Some examples of fluid solvers used in the realm of PWFA are Architect [18], QFLUID [19], MARPLE [20], FLASH [21] and the code used for hydrodynamic optically-field-ionized plasma channels in [22].From the theoretical point of view, fluid models have traditionally been developed by neglecting thermal effects, i.e., by neglecting pressure terms in the Euler equations.This choice has been motivated first and foremost by the fact that initial electron thermal energy in the plasma is expected to be of the order of k b T i ∼ 10 eV [23,24], which is a small value when compared with the electron rest energy m e c 2 = 0.511 MeV (the initial thermal energy normalized to the electron rest mass, µ i = k b T i /(m e c 2 ), is the control parameter that is usually used to assess the importance of thermal effects), and secondly by the fact that these cold fluid models are easier to treat theoretically, providing even analytical results in some simplified cases [25][26][27].Nevertheless, there is a series of important reasons that drives the development of warm fluid models for PWFA.First, the wave-like solutions to cold fluid equations become singular in the presence of wide and highly charged driver pulses (Wave Breaking [28,29]).The presence of thermal effects may be one of the regularising mechanisms that mitigate the singularity [30][31][32][33], by altering the wakefield properties and allowing electron trapping in the wakefield [34].A significant heating is also expected in the post-wavebreaking dynamic [35].Second, studies targeted at the late stage dynamics of the process [36][37][38][39][40] point to the importance of electrons temperature (with particular emphasis on thermally driven ion-acoustic motion [39][40][41][42][43]) for the restoration of the equilibrium conditions after a driver pulse has passed in the plasma channel.Understanding the restoration conditions is pivotal to enable the possibility of having high repetition rates of driver pulses in order to create sustained accelerating fields.Third, although the aforementioned initial temperatures would not lead to meaningful divergences in behavior from the cold case (at least in the first wave periods [34,35]), a different situation is expected as many consecutive pulses are injected into the system.The energy deposited by every pulse would be partially transferred to the plasma [36,38], leading to significant increases in temperature (some estimates provide O(1) keV increases in post wave-breaking situations [39,40]).Lastly, we mention that temperature effects might be also relevant for the study of positron acceleration in quasi-hollow warm plasma channels [41,[44][45][46].If one wants to develop a warm fluid theory, often the only viable option when trying to derive analytical results, an immediate problem has to be tackled: additional fields are now present in the set of equations (namely the pressure tensor fields) and therefore a suitable closure has to be carefully selected.The closure problem has been studied in the community, and various models have been proposed [47][48][49][50][51][52] (see [53] for an outlook).In this paper we explore two closures.The first one [47] is based on the assumption that the underlying probability distribution solving for the Vlasov equation is a Maxwellian equilibrium, which is described in the relativistic framework by a Maxwell-Jüttner distribution [54].Hereafter, we will refer to it as the Local Equilibrium Closure (LEC).This choice leads to no entropy production and therefore grants the adoption of an isentropic equation of state to close for the pressure tensor field, that in this framework can be described as a single scalar quantity.The LEC is in principle not well suited for the description of early stage dynamics, as the restoring mechanism that drives the plasma back to an eventual initial equilibrium state (particle collisions) happens on longer time scales than the ones which are typical of the early phases of PWFA.Nevertheless, this closure would still retain a physical relevance when studying late stage dynamics (when particle collisions start to become relevant [40]).Furthermore, side-by-side comparisons against other closure schemes (or PIC solvers) might indeed reveal that the LEC is also helpful for qualitative assessments on early dynamics.The second closure (hereafter named warm closure -WARMC) is based on the idea proposed in [48,49,51] and later on reconsidered by [32,55,56]: the centralized moments equations obtained from the coarse-graining of the Vlasov equation are closed by neglecting the third order centralized moment, choice motivated by the assumption of weakly warm systems, and hence small momentum distribution variances.This leads to a closed set of equations that can be solved without having to make hypotheses on the underlying momentum probability distribution, except for the smallness of its variance.This gives also the possibility to evolve independently the various components of the pressure tensor, and hence to evaluate the expected momentum spread anisotropies: in fact, as the dynamics of the system is strongly focused in the acceleration direction, and no collisions are taken into account on these short timescales to regularize the process, momentum spread anisotropies are to be expected [57,58].The main target of this paper is to develop novel numerical schemes for the simulation of warm fluid models in the context of PWFA, for both the LEC and the WARMC.The schemes rely on the lattice Boltzmann (LB) method [59,60] to solve the fluid equations.LB is a popular numerical scheme commonly used in computational fluid dynamics as an alternate scheme to direct hydrodynamical solvers.Its formulation is rooted in the kinetic theory of gases, and this provides a strong physical basis to the method.Furthermore, the space locality of the calculations involved in the method makes this solver prone to multi CPUs and/or multi GPUs parallelization [59].In the past years LB has been generalized to work in many fields other than classical fluid dynamics [60], and it is now widely accepted as a numerical solver for many physics problems that rely on a set of continuum equations.The LB formulation used in this paper is the so-called moment matching LB, that solves systems described by advection diffusion equations [59,60].This is the very same formulation used in [61] to solve for the cold fluid models in PWFA.Here we extend such formulation to warm fluid models, hence we are pushing further the usage of the LB method in the context of PWFA.The LB schemes for warm fluids are coupled to a Finite Difference Time Domain (FDTD) scheme that solves for the electromagnetic fields [62,63]: we refer to [61] for a more in depth explanation of the coupling.The development of the LB schemes for two different warm fluids closures will enable us to look for thermal effects that depend on the adopted closure scheme, and furthermore to assess the importance of thermal spreads anisotropies showing side-byside comparisons between the two closures.This paper is organized as follows: in Sec. 2 we provide a basic introduction to the LB method, with particular emphasis on the moment matching procedure; in Sec. 3 basic equations for collisionless relativistic plasmas are reviewed and the two fluid closures LEC and WARMC are discussed in Sec. 4 and Sec. 5 respectively; numerical results will be presented in 6 and an outlook and a discussion on future perspectives is given in Sec. 7.

Lattice Boltzmann (LB) method
In this section, we introduce the Lattice Boltzmann (LB) method, which we use to solve the fluid equations both in the LEC and WARMC models.We first present the basics of the method, that are directly drawn from the kinetic theory of gases and in their original formulation tasked to the reproduction of classical (non relativistic) Navier-Stokes equations.We then illustrate how the whole procedure can be adapted to the simulation of generic advection equations (moment matching LB).We will see in the following sections how the relativistic warm fluid equations, in both the LEC and WARMC, can be recast into a set of advection equations, and hence can be numerically solved via moment matching LB.The interested reader might look into [59,60] for more detailed discussions on both LB in its original formulation and its moment matching variant.

Basic LB
As the theoretical cornerstone of LB is kinetic theory, the natural starting point in its algorithmic derivation is the Boltzmann transport equation which describes the evolution of the phase space density f (x, ξ ξ ξ, t) referring to the number of fluid particles with velocity ξ ξ ξ at position x at time t.On the RHS of the equation, Σ(x, ξ ξ ξ, t) represents a production term that is usually expressed as the sum of two components: S(x, ξ ξ ξ, t) representing the action of volume body forces F on the particles, and Ω(x, ξ ξ ξ, t) the action of binary collisions: where the latter has been here written by recurring to the customary Bhatnagar-Groos-Krook operator [64], which expresses the tendency of the distribution function f to relax towards an equilibrium f eq in a typical relaxation time τ .The moments of the distribution function f that solves Eq. ( 1) are then used to retrieve the hydrodynamic fields.In this basic formulation, the zeroth-and first-order moments (i.e., n = 0, 1) provide the mass density and the momentum density of the fluid, respectively [59].It can then be verified through the Chapman-Enskog expansion [59,60,65] that the moments obtained according to Eq. ( 5) verify the target continuum equations -again, the Navier-Stokes equations in this original formulation of the LB method -provided that f is sufficiently close to f eq .The most pivotal step in the development of the LB method is the realization that one needs only a truncated version of the distribution functions f and f eq to properly recover the desired moments M (n) (x, t) that solve the field equations [66,67].For this reason, it proves expedient to expand Figure 1: Lattice Boltzmann discretization.Through a pruning procedure, the continuum velocity space is discretized in a minimal set of N pop velocities.This discrete set (stencil) is selected in such a way as to guarantee the hydrodynamic equations for coarse grained fields (e.g.density, momentum).On the right we report one of the most popular LB stencils, with N pop = 19 velocities in 3D, which we have also used in this work.In the table we report the adimensional velocity components of the stencil, together with corresponding statistical weights.The different shades of blue group velocities with the same magnitude.
both f and f eq into series of orthogonal Hermite polynomials in the variable ξ ξ ξ, and then to truncate this expansion up to the point where one recovers the macroscopic observables of interest [68].
The second most important ingredient in the development of the LB is the velocity discretization procedure.Because of the truncated Hermite polynomials expansion, one usually adopts Gauss-Hermite quadrature rules to prune the continuum velocity space and select a discrete set of ξ ξ ξ i (i = 0, . . ., N pop − 1) velocities (stencil), each of them having a statistical weight w i [69,70].In Fig. 1 we report the discrete stencil adopted in this paper.This velocity discretization procedure splits Eq. ( 1) into a set of N pop equations, one for every discrete distribution function f i = f (x, ξ ξ ξ i , t), also called population.The careful selection of the velocity stencil is performed with the goal of exactly preserving the coarse-graining of the moments when passing from the continuous integrals of Eq. ( 5) to discrete summations: The performed velocity discretization, when joined with an explicit time marching discretization of step ∆t, finally delivers the Lattice Boltzmann equation: where space is discretized on a regular lattice of characteristic length ∆x = ξ ξ ξ i ∆t.The production term Σ i (x, ξ ξ ξ i , t) is a discretization of the continuous counterpart appearing in Eq. (1): with S i (x, t) being discretized according to one of the many forcing schemes available in the literature [59] to properly reproduce Eq. ( 3).The algorithmic steps of the LB scheme are now clearly outlined.A set of N pop versions of Eq. ( 7), one for every f i , is evolved on a regular spatial grid: the f i are updated at every node with the source term Eq. ( 8) (Source step), and then stream to neighbouring nodes with their corresponding velocity ξ ξ ξ i (Streaming step).At every iteration the hydrodynamic fields are obtained through the discrete summation appearing in Eq. ( 6).

Moment Matching procedure for Forced Advection Diffusion Equation
In this section, we will make use of the framework established in Sec.2.1 to explain how the moment matching LB, tasked for the solution of advection diffusion equations (ADEs), works.This variant was initially developed as an alternative, more straightforward formulation of thermal LBs [71][72][73] and at the same time applied to model different physics phenomena, such as diffusive chemical reactions [74], combustion problems [75,76], dissolution in porous media [77,78] (more use cases can be found in [59,60]).The target equation for a moment matching LB is a forced ADE for the scalar field A(x, t): where the quantity A is advected with velocity u(x, t) and diffused via the diffusion parameter κ.The previously established framework (Sec.2.1) can then be used by considering a discrete probability distribution function f i (x, t) whose zeroth-order discrete moment is exactly the field A(x, t): The velocity stencil described in Fig. 1 can then be used to evolve via Eq.( 7) the distribution functions f i , and the Chapman-Enskog analysis grants that the coarse graining Eq. ( 10) solves for Eq. ( 9).The last remaining ingredient to complete the algorithmic discussion of the method is the form of the Hermite truncated expansion f eq i .It has been shown that a suitable form that recovers the correct equations for A and minimizes the computational error is [79]: where ν = 1 √ 3 ∆x ∆t is a reference lattice velocity.This completes the description of the main algorithmic steps of the moment matching LB.The diffusion coefficient κ is shown via Chapman-Enskog procedure to be dependent on the relaxation time τ in the following way: While in this work we keep τ as small as possible (taking into account the known inferior limit τ > ∆t/2 imposed by LB bulk stability conditions [59]) in order to recover the physics of the problem, we note from Eq. ( 12) that by adjusting the parameter τ in the simulations one might control the diffusive effect in the ADE, obtaining therefore a tunable regularizing effect that might be helpful when dealing with stiff advection equations.Finally, we report a couple of additional details for our implementation of the method.First, the source term S i is chosen to be of the form in order to correctly reproduce the forcing appearing in Eq. ( 9) [59,60].Second, open boundary conditions are imposed at the extrema of the domain via the enforcing of the variables f i : where x b and x f are respectively the boundary position and the nearest bulk fluid node.Lastly, we remark that the whole procedure showed so far can be adapted to a 3D axisymmetric geometry: when doing so, after switching to cylindrical coordinates, we postulate that every quantity is independent from the azimuthal variable A = A(r, z).Furthermore, no angular motion is requested, hence the azimuthal component of the advection velocity u is zero.Following [80,81] one can develop an axisymmetric moment matching LB by considering a 2D Cartesian geometry (r, z) where all the differential operators -the divergences and laplacians appearing in Eq. ( 9) -are carefully adapted to the cylindrical geometry [82].To this extent, we remark that differently from the cited sources we will consider A to be either a scalar or the coordinate component of some hydrodynamic tensor field.Hence particular care will have to be taken when converting the differential operators to the cylindrical geometry.Lastly, when adopting the cylindrical geometry we have to modify the conditions applied at the r = 0 boundary node.If the LB advected quantity A is symmetric with respect to radial reflections, we keep employing Eq. ( 14).Instead, if A changes sign under radial reflections (the radial momentum component has this property) , we employ the following condition: with i being the mirrored direction to i with respect to the radial axis.

Relativistic Kinetic Equations for collisionless plasma
In this section we review the basic equations that can be used to build hydrodynamic models of warm plasmas starting from the kinetic theory of gases, all expressed within the framework of special relativity [15][16][17].In the next sections, we will see how the obtained set of equations can be closed via either a local equilibrium assumption or a warm closure, and then explain how to recast them into a set of advection equations.
In the following we will work in a flat space-time with Minkowski metric signature η αβ = diag(+, − − −).When expressing formulas in a manifestly covariant form, we will adopt Einstein's summation convention, with Greek indexes running from 0 to D, and Latin indexes from 1 to D. When not explicitly stated, we will use ∂ α ≡ ∂ ∂x α .The starting point in a relativistic fluid theory of warm plasmas is the Relativistic Vlasov equation: where x α = (ct, x) is the space-time coordinate vector, ρ α = (ρ 0 , ρ ρ ρ) is the relativistic kinetic momentum vector, F αβ is Maxwell's electromagnetic field tensor and f = f (x α , ρ α ) is the single-particle probability distribution function, describing the number of particles of mass m and charge q (m = m e and q = −e in case of electrons) in the Lorentz invariant momentum space volume dρ ρ ρ ρ 0 .The fields subject of hydrodynamic theories emerge from the coarse graining of the distribution function f in relativistic momentum space: with the first moments being assigned the following specific names: Invariant density h, Particle flow N α and Energymomentum tensor T αβ .Eq. ( 16) implies the following conservation equations for the moments [15,17,56]: Eqs. (18) to (20) are the constitutive equations of a relativistic fluid, which must be coupled to Maxwell equations [83] to provide a complete description of a plasma.Eqs.(18) to (20) express conservation of mass, energy and momentum once a proper decomposition for N α and T αβ is provided: as an example, the cold fluid equations can be obtained from Eqs. ( 18) and ( 19) by assuming the following: U β = γ(c, u) being the fluid four velocity and n 0 the particle number density as measured in the fluid rest frame (density transforms under Lorentz boosts as n = γn 0 ).γ is the Lorentz factor associated to u.At variance with the cold fluid treatment, any other fluid model involving a non-zero temperature requires additional closing equations.
In the next sections, we will show two possible closures: the first one (LEC) postulates that the underlying distribution function describing the warm plasma (and hence solving Eq. ( 16)) is a local equilibrium distribution.This grants the use of an isentropic equation of state to close for the system; the second one (WARMC) rephrases Eqs.(18) to (20) as centered moments equations, and then assumes the third order centered moment to be zero on the excuse of small thermal spreads, i.e. by assuming small values of the initial thermal energy (i subscripts stand for the initial values) to electron rest mass This is the control parameter that is usually used to assess the importance of thermal effects in warm fluid models, and will be extensively used in the next sections to label our results.

Local Equilibrium Closure (LEC)
The key assumption of the LEC hypotesis is that the distribution f solving for Eq. ( 16) is the relativistic version of the Maxwell-Boltzmann distribution, the Maxwell-Jüttner distribution [54,84].This has a number of important consequences, first and foremost the fact that under this assumption the plasma can be considered to be an ideal fluid.
It can be shown [15] that in this case the particle flow and the energy momentum tensor assume the following forms: where ϵ 0 and P 0 are respectively the plasma energy density and pressure, as measured in the fluid's rest frame (all quantities written with the 0 subscript have to be intended as they are measured in this frame).Consequently, the conservation equations Eqs. ( 18) and ( 19) can be re-expressed as the relativistic counterpart of Euler's equations: where h e = (P 0 + ϵ 0 )/n 0 is the relativistic enthalpy per particle and p = mγu is the relativistic fluid momentum.Eq. (27) are not yet closed, but another important feature of the LEC assumption is that there is no entropy production in the fluid.Therefore to close the set of equations instead of a statement of energy conservation we consider the entropy per particle here written in the small temperature limit [15], to be constant.This statement leads to the following scaling for temperature: Furthermore, this grants the use of the Synge equation of state [85], that we consider here for small temperatures: . This provides sufficient information for the closure of Eq. ( 27), and paves the way to a successful numerical treatment of this set via moment matching LB.

Local Equilibrium Closure Lattice Boltzmann (LEC-LB)
We can now express Eq. ( 27) as a series of forced advection equations: where This set can then be numerically solved via moment matching LB, using the techniques explained in Sec.2.2.As already stated, we can adapt the method to work in a 3D axisymmetric environment, and particular care has to be taken when translating the equations to said geometry.As the final expressions are rather bulky, we report them in full in Appendix A.
The two-way coupling of the fluid with the electromagnetic fields needs no particular explanation: the plasma hydrodynamic quantities (particle density n and the transport velocity u) are obtained from the LB evolved quantities (the various components of A) and then fed to the FDTD Maxwell solver together with the driving bunch terms.The evolved electromagnetic fields are then plugged into Eq.( 34) as source terms, and the next LB iteration can be performed.
The only detail worth of discussion is the determination of the transport velocity u from the LB-advected quantities A. In fact, due to the appearance of the γ Lorentz factor in Eq. ( 33), one cannot easily obtain the relationship between u and the second component of the vector A (the one containing momentum p), and therefore there is the need for a specific iterative algorithm to determine it.Initially set p (0) as the zero-th order moment of the distribution function coming out of the LB iteration, divided by n: Then: 3. Set p (k+1) as 4. Repeat until convergence is reached.

Warm Plasma Closure (WARMC)
The derivation given in this section closely follows [56].To explain how to derive the Warm Closure from the conservation equations Eqs.(18) to (20) one has to first define the thermal momentum w µ as The second and third order centered moments of the distribution function f , respectively θ µν and Q µνλ , are then defined as: Therefore one obtains: and due to the mass-shell condition ρ µ ρ µ = m 2 c 2 it follows that: The conservation equations Eqs. ( 18) to (20) become therefore The warm closure consists in taking Q αβγ = 0, based on the assumption of a small momentum spread: and this also provides an additional condition, obtained from Eq. ( 44): Note that the previous expressions have been written using Eq. ( 38) to express everything in terms of the more familiar quantities n 0 and U α .Eqs. ( 48) to (52) constitute a closed set of fluid equations, which may be solved for the plasma evolution through moment matching LB.

Warm Plasma Closure Lattice Boltzmann (WARMC-LB)
Eqs. ( 48) to ( 52) can be rephrased as advection equations by realizing that, for any given quantity G = 1, n 0 U β /h, θ βγ /h , one can write: Therefore one has: where The system of equations can now be simulated via moment matching LB, interpreting all terms on the RHS of equations as external forcings.The Fluid-Maxwell coupling discussed in Sec.4.1 applies also in this case, with the exception of the iterative scheme for the transport velocity u, since in this framework u can be naturally identified from the LB advected quantities.Also in this case we employ the axisymmetric description: the passages that are needed to adapt Eqs. ( 51), ( 52), ( 55) and (56) to this framework are rather lengthy but simple and we report the final full expressions in Appendix A. As a final remark, it can be worth mentioning how the tensor θ αβ can be decomposed in its transversal and longitudinal pressure components (w.r.t. the direction of the driving bunch, here chosen to be the z−axis), respectively P ⊥ and P ∥ = P ⊥ + ∆P , as to define a term of comparison against the isotropic pressure case provided by the LEC.We postulate [17,86] that in the fluid rest frame the anisotropic energy momentum tensor and the thermal momentum would read as T µν 0 = diag(ϵ 0 , P ⊥ , P ⊥ , P ∥ ) , It follows then from Eq. ( 41) that It is sufficient to apply a Lorentz boost Λ µ ν to derive this tensor in a generic reference frame: Other than showing how to obtain the rest frame values for both pressure terms, this form of the tensor is actually useful to determine that some of its components are null in the axisymmetric framework.This effectively reduces the number of advection equations that have to be solved via moment matching LB.

Numerical Results
In this section we present the current capabilities of the method by reproducing known analytical results and we also show side-by-side comparisons between the two closures, highlighting the emergence of momentum spreads anisotropies.We divide the showcase of our results in three different parts: in the first we start by considering temperature effects in a completely 1D scenario, where our equations of motions are made mono-dimensional by imposing translational symmetry along the radial directions (all derivatives w.r.t.transversal coordinates are therefore zero), imposing D = 1 in Eqs. ( 33) and ( 34) (this is equivalent to considering a 1D kinetic momentum space) and considering (1 + 1)−dimensional tensors in the WARMC case.In this 1D1V set-up, we compare our numerical result with known analytical solutions [30][31][32]55].Then, we move to the discussion of dispersion relations in a 1D3V setup [56,87]: there is still translational invariance along the transversal directions, but a 3D kinetic momentum space is considered (D = 3 in Eqs.(33) and (34) and θ µν tensors are (3+1) dimensional).Lastly, we consider full spatially resolved simulations in a 3D axisymmetric environment (3D3V).In this work, plasma ions are considered as an immobile background with constant plasma density n i = 10 16 cm −3 , as they are several order of magnitude more massive than plasma electrons.Furthermore, the driving electron pulse is modelled by a rigid bi-Gaussian density n b , with rms-sizes σ z = σ r = 25 µm and moving at the speed of light c (moving from right to left): The driving bunch is initially centered in z 0 and with a peak amplitude α such that the Normalized Charge Parameter representing the strength of the perturbation, is kept at a fixed desired value.Here k p is the plasma wave number, k p = ω p /c, with ω p = e 2 n i /(mϵ 0 ) the cold plasma frequency (m being the electron mass and ϵ 0 the vacuum permittivity).Unless explicitly stated, all the results presented from now on will consider a computational domain of size L r = 6/k p and L z = 30/k p , with cell resolutions ∆z = ∆r = 0.01/k p and computational time step ∆t = 0.001/ω p , for a simulated physical time of 30/ω p .The value of τ is tuned in every setup in order to gauge between the LB numerical stability conditions [59] and the smallness of the parameter κ introduced in Eq. ( 12), as to properly recover the physics of the problem: we have chosen a value of τ = 0.52∆t in the 1D setups, and τ = 0.53∆t in the 3D setups, and both these two values perfectly reproduce the available analytic solutions.As already mentioned, just like other common PWFA solvers [18] the electromagnetic fields are solved via an FDTD scheme [62,63] through numerical integration of the curl Maxwell equations (Faraday's and Ampere's laws in Eqs. ( 21) and ( 22)).It can be shown in fact [88] that when these equations are considered together with the continuity equation, the divergence Maxwell counterparts (Gauss laws in Eqs. ( 21) and ( 22)) are automatically satisfied, provided that the fields are correctly initialized.For this reason, we properly initialize the electromagnetic fields by solving analytically the Maxwell system with just a rigid Gaussian density (the driving bunch) in its rest-frame, and then boost the quantities to the lab-frame [89,90].KeV) 1D theory [31,32] (which coincide for this choice of the parameters and hence is represented here by a single solid black line).We also show for reference the analytic solution that can be derived from the cold fluid model (light purple).From top to bottom, we show results for the electron plasma density n, the electric field E (normalized w.r.t.E p = mcωp e ), and the plasma velocity u.All the curves are plotted with the co-moving variable ζ = z − ct on the x−axis.We show two significant regimes: Q = 0.01 on the left column and Q = 0.5 on the right column.The Gaussian driving bunch (green line) appearing in the top panels is vertically shifted by a unit factor for visualization purposes.

1D1V results
We adapt our scheme to work in a 1D environment by employing a value of σ r which is way bigger than the computational domain along the radial direction L r , so that the driving bunch of Eq. ( 61) reduces effectively to a single Gaussian profile along the z coordinate, and the system assumes a translational invariance along the r coordinate.Eq. ( 62) is therefore employed by imposing D = 1.Furthermore, the 1V condition is reached for the LEC by setting D = 1 in Eqs. ( 33) and (34), while in the WARMC case is sufficient to initialize to zero the transversal components of the θ µν tensor.We present in Fig. 2 the results for a numerical benchmark of the method obtained by comparing both the two fluid solvers against the analytic solutions that can be obtained by considering Eq. ( 27) and Eqs. ( 48) to (52) in a 1D1V environment.The corresponding equations are well studied in the literature: for the LEC see Sec.VI in [31] or equivalently [30]; for the WARMC see [32], and replace the laser driven case with a particle driven case.The theoretical solutions are obtained via finite differences integration of Eq. (6.4) in [31] (LEC) and Eq. ( 8) in [32] (WARMC).Although coming from different set of equations, these solutions are equivalent in the limit of small temperature.Being this the case, they are represented by a single solid black curve in Fig. 2. We observe that both the two methods are perfectly able to reproduce the theoretical results.In this simplified 1D scenario it is also possible to start appreciating some effects of temperature on the dynamics (the choice of µ i is done in order to evidence thermal effects in the first periods of the wave): by comparing against the analytic solutions of the cold fluid model (light purple in Fig. 2) it is possible to see that temperature leads to a decrease in the thermal wavelength of the plasma wave.As it will be seen in the next section, this effect is in stark contrast to what can be observed in a 3V environment, where the wavelength increases or decreases depending on the selected fluid closure.65) and (66).In both the two closures, our numerical codes well reproduce the theoretical results.We remark that theory predictions are operated in the linear regime (small Q perturbation) and in the assumption of small initial temperatures (small µ i values).

1D3V results
An increase in complexity with respect to the 1D1V case is represented by the 1D3V case, that is targeted specifically at the recovery of dispersion relations for the plasma wave, while still retaining a 1D simulation framework.In this context it is actually possible to extract for both closures a dispersion relation linking the frequency ω of the plasma wave to its wavenumber k.In the literature, such result has been obtained with various levels of approximations [87,91].In our case, a theoretical study can be conducted by considering Eq. ( 27) and Eqs. ( 48) to ( 52), linearly perturbing them and then obtaining the dispersion relation via Fourier analysis [56,92,93].The following formulas are obtained at order O(µ i ): Note that the WARMC result in the previous equation coincides with findings in [87] and represents the relativistic Bohm-Gross relationship for warm Langmuir waves, where the 5/2 correction term stands for thermal inertia i.e. hydrodynamic transport correction due to temperature.The 1D3V framework is such that only longitudinal waves are present in the system, and therefore one has phase velocity equal to the bunch velocity v ph = ω k = c.Therefore, in this context one obtains at order O(µ i ) At variance with what was observed in the previous section, Eq. ( 65) foresees a closure dependent behavior of the wave: with respect to the cold case (where ω = ω p ), LEC plasmas exhibit wavelengths decreasing with temperature, whereas WARMC plasmas wavelengths increase with temperature.This can be observed in Fig. 3, where we report a comparison between the theoretical predictions and the numerical results.In both the two closure schemes we observe a very good agreement.

3D3V results
We now move our analysis to fully spatially resolved simulations.In Fig. 4 we provide an initial comparison of the two fluid closures in a 3D3V axisymmetric environment.The temperature is set to an initial value of k b T i = 20 keV, and the chosen regime of the driving pulse is quasi-linear, Q = 0.5.In the figure, we show color plots for four significant quantities: the plasma number density n, the longitudinal fluid velocity u z , the longitudinal accelerating field E z and the transversal wakefield E r − cB ϕ .It is possible to appreciate that at variance with the usual dynamics of the cold case [61], the typical wave-like pattern of peaks and valleys is perturbed by the emergence of an acoustic behavior, that promotes electron motion out of the radial axis.The presence of pressure waves can be appreciated by the appearance of Mach cone structures [94], that can be easily spotted by looking at anyone of the panels of Fig. 4: as the driving bunch travels at the speed of light c, it perturbs the plasma medium.The perturbation propagates at a specific velocity c s (that depends on the initial temperature µ i ) smaller than c, and the wave-fronts of the perturbation, that in a sub-sonic flow would radially propagate all over the space, are instead constrained in a conical structure around the perturbation: a super-sonic Mach cone is observed.From a qualitative point of view, at this stage of the analysis there is no clear difference in behavior between the two fluid closures.A more quantitative characterization can be performed by analyzing these acoustic structures in the linear regime (i.e. for small values of the parameter Q), where it is possible to elaborate an analytical theory for both the two fluid models: in a nutshell, by using perturbation theory on respectively Eq. ( 27) and Eqs. ( 48) to (52), it is possible to derive a forced Klein-Gordon equation for the electron plasma density n, where the only difference between the fluid closures appear in the temperature dependent coefficients of the equation (see Appendix B for details).By looking at this equation the acoustic nature of the dynamics becomes clear, and it is possible to derive some analytical results to compare with the numerical simulations.In Fig. 5 we show the result of numerical simulations (top half of the panels) versus the analytical solutions (bottom half of the panel).There is good agreement between the two.By inspecting the two cited panels it is evident that the two fluid models, although exhibiting similar qualitative features, Figure 5: Comparison of simulations (top half of both panels) against analytic solutions of the warm linear theory (bottom half of both panels), for both the two closures (top panel LEC, while bottom panel WARMC).The quantity shown is the particle density n, normalized w.r.t. the initial plasma density n i .We also plot the Mach cone structure, together with its inclination Θ for reference.
differ in the inclination Θ (plotted for reference in the figure) of the conical structure.The theory of super-sonic flows [94] tells us that this inclination is linked to the intensity of the velocity c s via c s = c sin(Θ) , (67) which means that the wavefront of the perturbation propagates in the two models with different velocities.The warm linear theory provides in fact a dependency (look at Appendix B for more analytical details) of c s from the initial temperature: One can then run a set of simulations with different initial temperatures µ i , and study the inclination of the conical envelope (from which we extract c s ) as a function of the temperature.In Fig. 6 we show the results of this analysis, that again show good match between the simulations and expectations from theory.As a final element of discussion we present the results of Fig. 7, where we consider the rest frame pressures ratio P ∥ /P ⊥ for the WARMC model.This quantity is expected to be equal to one in the case of an isotropic pressure tensor, as it is the case in the LEC model, therefore in Fig. 7 we consider the following quantity:   68) and (69).In both the two closures, our numerical codes well reproduce the theoretical results.We remark that theory predictions are operated in the linear regime (small Q perturbation) and in the assumption of small initial temperatures (small µ i values).
The fluid rest frame pressures P ∥ and P ⊥ are obtained from the numerics by inverting Eq. ( 60) (the lab-frame components of the θ µν tensor are obtained from the simulations).We note thermal spread anisotropies are mostly relevant in the proximity of the bunch (located at ζk p = 3.0) where they become of the order of 10%, and become less and less important as one proceeds further away from the perturbation.To the best of our knowledge, this is the first time that such analysis is conducted on the WARMC model for a spatially resolved plasma (a study on a laser excited, 1D restricted plasma can be found in [95]) and this preliminary results indicate that the LEC model, that is missing this pressure anisotropy feature by construction, could still be used to characterize plasma behavior in late-stage dynamics studies.We reserve further investigation, in particular studies on the dependency of this anisotropic feature on driving bunch parameters and initial temperatures, for later research.

Outlook and perspectives
In the development of a fluid model for the simulations of PWFA processes, a multitude of physical ingredients have to be taken into account to provide a realistic description of the process.In an earlier paper [61], some of the authors started the exploration of lattice Boltzmann (LB) schemes for the construction of fluid models for PWFA and considered the simplified case of cold fluid models.This paper represents a step forward, in that we have explored Figure 7: Study of the thermal spread anisotropies in the WARMC model.The pressure ratio P ∥ /P ⊥ , whose value is obtained in the numerics via inversion of Eq. ( 60), is compared against its isotropic value 1. Differences of order 10% are observed in proximity of the bunch.The chosen initial temperature value is µ i = 0.04 → k b T i = 20 KeV (again, selected for enhancing thermal effects in the first wave periods) while Q = 0.5.
LB fluid schemes accounting for thermal effects [30-32, 34, 39, 41, 44-46].The inclusion of thermal effects is a rather non trivial task, due to the theoretical complication represented by the choice of a proper closure to the set of equations [53].We have handled this problem by selecting two of the most popular closure models that have been discussed in the literature: the first one relies on the assumption of a local equilibrium (LEC) [47], while the second one involves truncating the hierarchy of centered equations at an arbitrary order (WARMC) [56].We have then shown how to successfully adapt the LB schemes to both closures.If from one side the LEC is nominally not appropriate for a collisionless warm plasma, from the other side the WARMC is obtained under the assumption of an asymptotically small temperature.Any finite temperature, however small, can raise the question on what is the right closure scheme to obtain the correct fluid model for collisionless warm plasma dynamics.The preliminary comparisons shown in this work, although not presenting strong qualitative differences in the dynamics, give a clear indication that the selection of a closure scheme is pivotal for the quantitative assessment of PWFA experiments.To this aim, a oneto-one comparison between the predictions of fluid models and the PIC simulations (or a numerical solution of the Vlasov equations) could be helpful to shed lights on the matter.Work is in progress along this direction.The next logical step in the development of the method is the inclusion of ions dynamics.Since plasma ions are way more massive then electrons, their dynamics happens on longer time scales than the ones examined in PWFA studies that target the early stage evolution of system, unless a strongly dense driving bunch is considered [96].There is though a growing research interest for studies that inquire the late time evolution of the system where this kind of dynamics cannot be ignored anymore [36,38,40,97]: unlocking the movement of ions brings in a new wealth of physical processes such as soliton dynamic [40] and ion channels formation [39].Furthermore, most of these studies invoke thermal effects to explain the acoustic motion of the ions [38,40,41,43], and it is then clear how the development of a method that would be able of handling both temperature and ions dynamics is of the uttermost importance.The inclusion of ion motion in the numerical scheme presented in this paper is theoretically trivial, and only brings in a bigger computational effort (more equations need to be integrated): the aforementioned characterization is therefore completely within the reach of the LB method.We also would like to mention a couple of advancements that would make the method presented in this paper a more complete numerical tool for the simulation of wakefield acceleration processes: the first is represented by the possibility of including a laser source field as the plasma perturbating mechanism, as the LWFA has always been an alternative (w.r.t.PWFA) route to wakefield acceleration [8]; the second is the possibility of handling non-rigid particle bunches, as the tracking of the driving bunch properties is both an important element of diagnostic and also a feature to be kept under control in modern day PWFA experiments (see for example [98]).Again, work in this direction is in progress.On the computational side, we would like to mention that the code developed for this work is amenable to a number of numerical optimizations (such for example porting on GPUs) and it is a firm intention of the authors to realize these improvements in future works.As the most important design choice presented in this paper is the selection of the LB method as a fluid solver, a few comments are in place.In this paper we have shown that the LB method is well capable of recovering analytic results, and we would like to mention further advantages that could make LB a convenient and viable option in the context of PWFA simulations.Although PIC methods are able to model the dynamics at the particle level, and hence are able to model kinetic effects, they constantly have to keep at bay their inherent numerical noise [9,99].Therefore in some situations it would be preferable to dispose of alternative tools for prototyping, and to reserve more complete PIC simulations for later stages.Solvers based on a strict discretization of the Vlasov equations would not suffer of the noise problem and still be able to capture mesoscopic effects, but for any dimensionality more than 1D they would be too computationally demanding and hungry for memory resources.Fluid solvers (and the LB method among them) present then a viable option for the realization of quick PWFA simulations.They are coarse-grained and hence do not suffer of numerical noise, and still capture a good amount of the physical effects that are encountered in PWFA.Furthermore, the LB method lends itself exceptionally well to parallelization on both CPUs and GPUs [59].In our implementation, we achieved parallelization on multi CPUs using the Message Passing Interface (MPI).This approach involves dividing the computational domain into multiple rectangular sub-domains, corresponding to the number of processors.Leveraging the local nature of computations ( Eq. ( 7), r.h.s.), communications are solely required to exchange populations between neighboring processors during the "streaming' process ( Eq. ( 7), l.h.s.).By effectively disentangling "compute" and "communicate", this strategy greatly enhances the parallelization process [100][101][102][103]. Simulations in this study were conducted on an Intel Xeon E5-2695@2.40GHz processor.A representative simulation (like those shown in Fig. 4-5) run on 96 processors requires approximately 182 minutes for the LEC-LB model and 368 minutes for the WARMC-LB model for 3 • 10 4 time steps.Memory requirements for such simulations are ∼ 1.6 GB for the LEC-LB model and ∼ 4.3 GB for the WARMC-LB model, although these numbers could be decreased by further optimizations of the code.As a final note, we would like to remark that LB takes roughly 50% of the compute time in the simulation, whereas the Maxwell solver for the electromagnetic fields takes 1%.This is to be expected as the amount of computation appearing in LB is significantly higher.Fig. 8 presents the execution time per iteration and speedup data for varying numbers of processors (strong scaling).Furthermore, an advantage of the LB method over Vlasov solvers is that the first adopts a smart pruning of the velocity space [59,60], thus improving the computational efficiency.It is possible to increase the number of discrete velocities N pop (see Sec. 2.1).However, this leads to a proportional increase in both the computational cost and memory requirements, making it a crucial factor to consider when choosing the LB stencil.Common practice is indeed to choose the minimum value N pop that ensures the recovery of the hydrodynamic properties [59].We hasten to remark, however, that at variance with usual hydrodynamic solvers, LB's theoretical formulation is strongly grounded in kinetic theory, and its fluid behavior is obtained via the smart discretization of the velocity space cited above.This is a remarkable feature that makes the method a strong candidate for the inclusion of kinetic effects into PWFA fluid solvers.In fact, recent studies [104][105][106] in other research fields show that LB is indeed capable of capturing behaviors beyond hydrodynamics just by increasing the number of discrete kinetic velocities.This could open novel perspectives for developing more refined numerical schemes for simulations of PWFA processes, with the obvious need of a precise comparison/benchmark against some reference PIC/Vlasov simulations.Further investigation on this is in progress.We finally would like to mention that an alternative route to LB wakefield simulation could be represented by the Relativistic Lattice Boltzmann method [107].This is an extension of LB hydrodynamic schemes (like the ones exposed in Sec.2.1) to the theory of special relativity, originally created for simulation in astrophysical [108] and condensed matter [109,110] contexts, and therefore constitutes a promising tool for the simulation of warm plasmas within the LEC assumption.Adapting it to a PWFA framework is though more technical and less immediate w.r.t. the moment matching LB used in this paper, therefore the authors reserve further development on this line for future works.

B Warm linear theory
The fluid equations, respectively Eq. ( 27) and Eqs. ( 48) to (52) for the two closures, can be linearly perturbed w.r.t. the initial rest state, when coupled with Maxwell equations Eqs. ( 21) and (22).One then obtains a forced Klein-Gordon equation for the density perturbation n 1 = n − n i : Where the coefficients c s , a e and a b depend on the initial temperature in a way that's dictated by the selected closure scheme.At order O(µ i ), one has: This result, aside from giving an explicit dependency of the acoustic velocity c s w.r.t.temperature (result Eqs. ( 68) and ( 69)), tells us that the behavior of the two fluid schemes is qualitatively the same.One can further inspect Eq. ( 92) to determine a dispersion relation for the Langmuir wave.Assuming to be far away from the driving bunch (so that n b can be ignored), one derives the following by recurring to Fourier transformation (∂ t , ∇ ∇ ∇) → (−iω, ik): which reduces to the result Eqs. ( 65) and ( 66) when considered at order O(µ i ).Last thing to discuss is the strategy to the solution of Eq. ( 92).After performing the traveling wave ansatz, which states that every field is dependent on z and t only through the co-moving coordinate ζ = z − ct, we obtain the following equation: At this point, we introduce the Hankel Tranformation of order 0 [111]: where J 0 (wr) is the Bessel Function of first kind and order 0. This transform behaves nicely in the presence of transversal laplacians of radial functions (∇ ⊥ ) → (iw).Therefore, if one Hankel transforms Eq. ( 96), a Forced Harmonic oscillator in the variable ζ is obtained: which can be solved via the Gauss-Green method.The final anti-transformed solution, to be numerically integrated, reads therefore as:

Figure 2 :
Figure 2: Comparison between the numerical results obtained from the moment matching LB with the LEC (LEC-LB), the moment matching LB with the WARMC (WARMC-LB) and their respective warm (µ i = 0.04 → k b T i = 20 KeV) 1D theory[31,32] (which coincide for this choice of the parameters and hence is represented here by a single solid black line).We also show for reference the analytic solution that can be derived from the cold fluid model (light purple).From top to bottom, we show results for the electron plasma density n, the electric field E (normalized w.r.t.E p =

Figure 3 :
Figure 3: Dispersion relation dependency on temperature.Comparison of the numerics obtained from the moment matching LB simulations versus the theory predictions Eqs.(65) and(66).In both the two closures, our numerical codes well reproduce the theoretical results.We remark that theory predictions are operated in the linear regime (small Q perturbation) and in the assumption of small initial temperatures (small µ i values).

Figure 4 :
Figure 4: Qualitative comparison of 3D3V LEC-LB simulation against the WARMC-LB counterpart.We show four significant quantities, all properly adimensionalized w.r.t.previously defined values: the particle density n, the longitudinal fluid velocity u z and both the two components of the wakefields, E z and E r − cB ϕ .The chosen initial temperature value is µ i = 0.04 → k b T i = 20 keV (again, selected for enhancing thermal effects in the first wave periods) while Q = 0.5.

Figure 6 :
Figure 6: Acoustic velocity c s dependency on the initial plasma temperature µ i .Comparison of the numerics obtained from the moment matching LB simulations versus the warm linear theory predictions Eqs.(68) and(69).In both the two closures, our numerical codes well reproduce the theoretical results.We remark that theory predictions are operated in the linear regime (small Q perturbation) and in the assumption of small initial temperatures (small µ i values).

Figure 8 :
Figure 8: Top panel: execution time per iteration (in [s]) as a function of the number of processors for both LEC-LB (blue dotted) and WARMC-LB (orange solid).The scale is logarithmic.Bottom panel: corresponding speedups (data from the two closures are mostly overlapping).The simulation setup is the same as for results in Fig. 4-5.