Physics of relativistic collisionless shocks: The scattering center frame
PPhysics of relativistic collisionless shocks:The scattering center frame
Guy Pelletier, Laurent Gremillet, Arno Vanthieghem,
3, 4 and Martin Lemoine Universit´e Grenoble Alpes, CNRS-INSU, Institut de Plan´etologieet d’Astrophysique de Grenoble (IPAG), F-38041 Grenoble, France CEA, DAM, DIF, F-91297 Arpajon, France Institut d’Astrophysique de Paris, CNRS – Sorbonne Universit´e, 98 bis boulevard Arago, F-75014 Paris Sorbonne Universit´e, Institut Lagrange de Paris (ILP), 98 bis bvd Arago, F-75014 Paris, France (Dated: July 26, 2019)In this first paper of a series dedicated to the microphysics of unmagnetized, relativistic col-lisionless pair shocks, we discuss the physics of the Weibel-type transverse current filamentationinstability (CFI) that develops in the shock precursor, through the interaction of an ultrarelativisticsuprathermal particle beam with the background plasma. We introduce in particular the notionof “Weibel frame”, or scattering center frame, in which the microturbulence is of mostly magneticnature. We calculate the properties of this frame, using first a kinetic formulation of the linear phaseof the instability, relying on Maxwell-J¨uttner distribution functions, then using a quasistatic modelof the nonlinear stage of the instability. Both methods show that: (i) the “Weibel frame” moves atsubrelativistic velocities relative to the background plasma, therefore at relativistic velocities relativeto the shock front; (ii) the velocity of the “Weibel frame” relative to the background plasma scaleswith ξ b , i.e. , the pressure of the suprathermal particle beam in units of the momentum flux densityincoming into the shock; and (iii), the “Weibel frame” moves slightly less fast than the backgroundplasma relative to the shock front. Our theoretical results are found to be in satisfactory agreementwith the measurements carried out in dedicated large-scale 2D3V PIC simulations. I. INTRODUCTIONA. Motivations and objectives
The Weibel-type current filamentation instability(CFI) that develops in anisotropic plasma flows [1–5], hasgained a lot of attention in past decades because of itsrelevance to various fields of physics. Through the gener-ation of skin-depth-scale current filaments surrounded bytoroidal magnetic fields, which gradually isotropize thecounterstreaming plasmas, this instability is a key pro-cess in the formation of astrophysical collisionless shockwaves [6]. Because it also arises in the precursor ofweakly magnetized shock waves, i.e. , the upstream re-gion where the unshocked plasma flows against a beamof suprathermal particles, it further sustains the shocktransition and likely generates the magnetized turbulencethat is required for particle acceleration and radiation.Such microphysics may well underpin the outstandingphenomenon of gamma-ray burst afterglows, in which asubstantial fraction of the ∼ ergs liberated by thecataclysmic event is radiated away by shock-acceleratedelectrons, on hours to months timescales [7, 8].In laser-driven high-energy-density physics and labo-ratory astrophysics, the CFI also plays a central role [9–17]. In particular, it has been recently observed to growin the interpenetration of laser-ablated plasmas [18, 19],marking a first step towards the generation of collision-less shocks in the laboratory [20–23].The development of the CFI in symmetric collidingplasmas has been studied theoretically in both the linear[5, 7, 24, 25] and nonlinear [26–31] regimes, and numer-ically using particle-in-cell (PIC) simulations, [3, 32–38]. In such a symmetric configuration, the lab frame appearsas a privileged frame in which to discuss the physics ofthe CFI and it is found, indeed, that the CFI generatesan essentially magnetic structure in that frame.However, in the precursor of relativistic collisionlessshocks, the CFI develops in a highly asymmetric config-uration: in the reference frame in which the shock frontlies at rest (the shock rest frame R s ), the backgroundplasma appears cold, and streams at relativistic veloci-ties through a quasi-isotropic, ultrarelativistically hot gasof suprathermal particles. Conversely, in the backgroundplasma rest frame R p , the suprathermal particles form adense beam of angular dispersion (cid:28)
1, carrying high iner-tia, and moving through a tenuous gas of sub- or mildly-relativistic temperature. The growth of the CFI in theprecursor of relativistic collisionless shocks has been dis-cussed in a number of numerical studies [39–46], but onlya few theoretical studies have discussed its properties insuch conditions [47–50]. In particular, the notion of aframe in which this CFI is essentially of magnetic naturehas not, to the best of our knowledge, received attentionso far; the only exception being Ref. [23] where this framewas introduced to simplify the calculation of the CFI inthe precursor of a nonrelativistic electron-ion shock.The present paper is the first of a series in whichwe discuss the microphysics of unmagnetized, relativis-tic collisionless pair shocks. Here, we address in detailthis notion of a “Weibel frame”, in which the electro-magnetic configuration is essentially magnetic in nature.As shown in accompanying papers of this series, and inparticular [51], this reference frame plays a fundamentalrole in the physics of collisionless shock waves. In Pa-per II [52], it is shown how the noninertial character of a r X i v : . [ a s t r o - ph . H E ] J u l the “Weibel frame” controls the heating and slowdownof the background plasma. In Paper III [53], this ref-erence frame is used to calculate the scattering rate ofsuprathermal particles; the latter quantity controls theresidence time of suprathermal particles in the upstream,hence the acceleration timescale, thus the maximum en-ergy etc. Therefore, this notion of a “Weibel frame” hasdirect phenomenological consequences and potential ob-servable radiative signatures. Finally, in Paper IV [54],we discuss the growth of the filamentary turbulence in theprecursor of a collisionless shock. In the following, as inall accompanying papers, we support our present theoret-ical findings through detailed comparisons to dedicatedlarge-scale 2D3V PIC simulations of relativistic collision-less shocks.The present paper is laid out as follows. Section IIintroduces the notion of the “Weibel frame” in a fluidmodel of the linear phase of the current filamentation in-stability. Section III then extends those results to a fullykinetic model of the CFI; the details of these calculationsare provided in Apps. A, B and C. Section IV introducesthe notion of the “Weibel frame” in a nonlinear modelof the filamentation phase, in which the filaments aremodeled in quasistatic equilibrium at all points in theprecursor. Finally, Sec. III D and IV B provide detailedcomparisons of those various models to dedicated PICsimulations. A summary and conclusions are provided inSec. V. We use units such that k B = c = 1. The metricsignature is ( − , + , + , +). B. Setup
In this work, we consider an unmagnetized, relativisticcollisionless shock wave propagating through an electron-positron plasma. In the shock rest frame R s , the back-ground plasma is initially incoming at relativistic velocity β ∞ (cid:39) − γ ∞ (cid:29)
1) along the x − axis from+ ∞ . Inside the precursor, its mean velocity is written β p (Lorentz factor γ p ). The precursor is defined as theregion permeated by a beam of suprathermal particles,which were reflected off the shock surface or acceleratedthrough a Fermi-like process, possibly up to high ener-gies.Our conventions are as follows: the index p refersto the background plasma, while the index b refers tothe suprathermal particle population. For species α , n α represents the density, w α the enthalpy density, p α the pressure, and T α the temperature, all measured inthe species rest frame, while u αµ = γ α (1 , β α ) repre-sents its four-velocity. All throughout, α will denote aspecies (‘p’ or ‘b’) and not a space-time index. Boththe background plasma and the suprathermal beam arecomposed of electrons and positrons, so that there are ac-tually four species. However, unless otherwise noted, wedo not make any particular distinction between electronsand positrons in a given population, and therefore dropany reference to the charge. The above proper hydrody- namic quantities nevertheless correspond to one chargedspecies (electrons or positrons) of one population (eitherthe background plasma or the suprathermal beam). Asin other papers of this series, we rely on the followingframes: the shock front rest frame R s , the backgroundplasma rest frame R p , the downstream rest frame R d and the “Weibel frame” R w , to be determined hereafter.Quantities evaluated in one or the other reference frameare indicated with the respective subscripts | s , | p , | d and | w . By default, however, frame-dependent quantities thatlack subscripts are defined in the lab frame, which coin-cides with the shock rest frame R s . Finally, thermody-namic moments n α , p α , T α and w α are always defined asproper, unless otherwise stated; they are thus defined inthe rest frame of species α .According to the fluid shock jump conditions [55],the typical temperature of shocked particles is T b = κ T b γ ∞ m e , where κ T b = ( (cid:98) Γ b − − (cid:98) Γ b ) / (cid:98) Γ − / , in termsof the adiabatic index (cid:98) Γ b of the shocked gas; the latter isrelativistically hot, so (cid:98) Γ b = 4 / (cid:98) Γ b = 3 / κ T b = 1 / (3 √
2) in 3Dand κ T b = 1 / (2 √
3) in 2D. Strictly speaking, the tem-perature of suprathermal particles is by definition largerthan that of the shocked plasma, defined above. Yet, asthis temperature of suprathermal particles always scaleswith γ ∞ , we retain the above definition, emphasizing that κ T b is generically larger than the above, by about anorder of magnitude, so that κ T b ∼ a few. This is illus-trated in Fig. 1, which shows the profiles of the mainhydrodynamic properties of the background plasma andsuprathermal beam for two reference PIC simulations,with respective Lorentz factors γ ∞ = 17 and γ ∞ = 173.Our PIC simulation are described further on inSec. III D. Let us note here, however, that the refer-ence frame of these 2D3V simulations (2D in configura-tion space, 3D in momentum space) coincide with thedownstream rest frame R d . Hence, a Lorentz factor γ ∞ = 17 (resp. γ ∞ = 173) in the R s frame correspondsto a simulation frame Lorentz factor γ ∞| d = 10 (resp. γ ∞| d = 100). Simulations are conducted in the x − y plane, with (cid:98) x oriented along the shock normal, and (cid:98) y defining the transverse dimension.The ultrarelativistic suprathermal particle gas is nearlyisotropic in the shock rest frame, meaning a drift Lorentzfactor γ b ∼ p b ,written ξ b in units of the incoming momentum flux atinfinity F ∞ = γ ∞ β ∞ n ∞ m e : p b = ξ b F ∞ . (1)The leading instability that develops in the precur-sor of weakly magnetized ultrarelativistic shocks is theWeibel-type CFI described above. In principle, thisinstability is defined in all momentum space ( k (cid:107) , k ⊥ )where k (cid:107) = k · (cid:98) x and k ⊥ = k ⊥ (cid:98) y in our PIC sim-ulations. The purely transverse modes correspond to - - - - - - - - - - - - - FIG. 1. Spatial profiles of the main hydrodynamic quantities characterizing the background plasma and the beam of suprather-mal particles in the shock precursor, as a function of distance to the shock front x | d (simulation frame). Upper row, panels (a),(b) and (c): profiles extracted from a PIC simulation of shock Lorentz factor γ ∞ = 17 (corresponding to a relative upstream-downstream Lorentz factor γ ∞| d = 10) at time t = 3600 ω − ; lower row, panels (d), (e) and (f): from a PIC simulation of shockLorentz factor γ ∞ = 173 ( γ ∞| d = 100) at time t = 6900 ω − . Panels (a) and (d): Lorentz factor of the background plasma inthe simulation (downstream) rest frame (dark gray), of the beam in the simulation frame (light red) and the relative Lorentzfactor between the beam and the background plasma (light blue); panels (b) and (e): proper temperature of the backgroundplasma (dark gray) and of the beam (light red); panels (c) and (f): proper density of the background plasma (dark gray) andof the beam (light red). The asymmetry of the beam-plasma configuration is manifest in those figures. k (cid:107) (cid:28) k ⊥ , while the so-called oblique modes correspond tothe limit k (cid:107) (cid:39) k ⊥ (cid:29) ω p , with ω p = (cid:0) πn ∞ e /m e (cid:1) / theplasma frequency of (one charged species of) the unper-turbed far-upstream background plasma. Oblique modesare likely relevant far in the precursor, but most likelyLandau damped once the background plasma heats up[38, 49, 50]. We thus assume that the transverse modesdominate in most of the precursor and restrict ourself tothe limit k (cid:107) → δE x , δE y and δB z in our two reference PIC simulations,normalized to the incoming momentum flux at infinity F ∞ . The purely transverse CFI generates a transverse δB z and its associated electrostatic δE y component; bydefinition of the “Weibel frame” R w (see further on), δE y | w = ,
0, and hence δE y = β w | d δB z in the R d simula-tion frame. The growth of this instability also generatesan inductive longitudinal electric field δE x . Nonstrictlytransverse CFI modes (with k (cid:107) (cid:54) = 0) further generate anelectrostatic longitudinal δE x component. Figure 2 re-veals that the energy density in the transverse field com- ponents dominates that in the longitudinal electric fieldeverywhere in the precursor in the R d frame, and thatthe electric field is smaller than the magnetic field in thenear precursor where the ratio can be measured accu-rately (up to x of the order of a few hundred to a thou-sand c/ω p ). Finally, the transverse magnetic field energyremains larger than the longitudinal electric field compo-nent when deboosted back to the R w frame. Therefore,this “Weibel frame” appears well defined in the 2D3Vsimulations.Let us emphasize that the limit k (cid:107) (cid:54) = 0 certainlyremains of interest; as a matter of fact, a nontrivialstructure along the streaming axis turns out to be amandatory requirement to achieve pitch-angle diffusionof suprathermal particles [5, 53]. Yet, we expect thatthe main features of the (already nontrivial) calculationsthat follow will remain valid in a limit 0 < k (cid:107) (cid:28) k ⊥ . Wethus consider purely transverse electromagnetic pertur-bations, i.e. , δA µ = ( δ Φ , δA x , , ∂ x = 0 for allquantities.Figure 2 suggests that the ratio (cid:10) δE y (cid:11) / / (cid:10) δB z (cid:11) / ,namely β w | d , depends on x , which implies that the“Weibel frame” is not globally inertial. This frame ac-tually decelerates from the far to the near precursor, as - - - - - - - - - - - - - - - - - - - - FIG. 2. Spatial profiles of normalized electromagnetic field en-ergy densities in transverse magnetic field δB z (black), trans-verse electric field δE y (blue) and longitudinal electric field δE x (red), as a function of distance to the shock. Panels(a) and (b): extracted from simulation γ ∞ = 17 (simulationframe γ ∞| d = 10). Panels (c) and (d): simulation γ ∞ = 173(simulation frame γ ∞| d = 100). Panels (a) and (c) show theenergy densities measured in the simulation frame; panels (b)and (d) show the corresponding energy densities deboostedto the “Weibel frame” R w (where δE y | w vanishes by defini-tion). The velocity of R w is measured in the simulation as β w | d = (cid:104) δE y (cid:105) / / (cid:104) δB z (cid:105) / , where the average is taken overthe transverse dimension of the simulation box. See text fordetails. discussed in [51]. The noninertial nature of R w has im-portant consequences for the physics of the shock, mostnotably the deceleration and heating of the backgroundplasma, which form the focus of a subsequent paper inthis series [52].In the present paper, we characterize the velocity β w | d (more specifically, β w | p ) at each point x of the precur-sor, given the physical conditions at this point. The dis-cussion that follows thus relies on an implicit WKB-likeapproximation, which stipulates that the ‘Weibel frame”has time to adjust at each point to the local physical con-ditions on a deceleration length scale | u w d x/ d u w | . Im-proving on this assumption would necessitate a properinclusion of noninertial effects, which are characterizedby the velocity profile d u w / d x , in the calculations thatfollow. However, this velocity profile is itself deter-mined by the response of the background plasma andthe suprathermal beam to the microturbulence. Such anendeavor thus represents a rather formidable task, wellbeyond the scope of any current study. II. THE “WEIBEL FRAME” IN A FLUIDMODEL
Here we present a fluid derivation of the CFI and ofits associated “Weibel frame”, in the context of the pre-cursor of an unmagnetized, relativistic electron-positronshock. Although Sec. III D will demonstrate that, in ac-tual relativistic shock precursors, kinetic corrections aremandatory to describe the physics of the CFI, the fol-lowing fluid model retains the advantage of simplicity aswell as a pedagogical virtue. A fully kinetic descriptionof the CFI and its “Weibel frame” in the case of Maxwell-J¨uttner plasma distribution functions will be provided inthe forthcoming Sec. III.To preserve covariance, we use a relativistic fluid for-malism. The conservation of the total energy-momentumtensor can be written in the compact way [5]: w α u αµ ∂ µ u αν + h αµν ∂ µ p α = q α n α u αµ F νµ , (2)introducing h αµν = η µν + u αµ u αν , which projects or-thogonally to u αµ (since u αµ u αµ = − w α u αµ ∂ µ δu αν + h αµν ∂ µ δp α = q α n α u αµ δF νµ . (3)In the following, the four-velocity perturbation is de-composed as δu αµ = (cid:0) δγ α , γ α δβ αx , γ α δβ αy , γ α δβ αz (cid:1) ,with δγ α = γ α β α δβ αx . Using the short-hand notation δ β ⊥ α = (0 , δβ αy , δβ αz ) and ∇ ⊥ = (0 , ∂ y , ∂ z ), the sys-tem can be rewritten explicitly as γ α w α ∂ t δβ αx + β α ∂ t δp α = − q α n α γ α ∂ t δA x , (4) γ α w α ∂ t δ β ⊥ α + ∇ ⊥ δp α = − q α n α γ α ∇ ⊥ ( δ Φ − β α δA x ) . (5)Current conservation written to first order also yields γ α ∂ t δn α + n α ∂ t δγ α + γ α n α ∇ ⊥ · δ β ⊥ α = 0 . (6)The pressure perturbation can be related to the densityperturbation through the adiabatic index (cid:98) Γ α : δp α = (cid:98) Γ α p α n α δn α . (7)In the following, we use the isentropic sound speedsquared c α ≡ (cid:98) Γ α p α w α , (8)so that δp α = c α ( w α /n α ) δn α . Given the relation be-tween δγ α and δβ αx , Eq. (4) can then be used to expressthe pressure perturbation in terms of the perturbed ap-parent density δN α ≡ δ ( γ α n α ): ∂ t δp α = c α γ α (1 − c α β α ) (cid:20) q α γ α n α β α ∂ t δA x + γ α w α n α ∂ t δN α (cid:21) . (9)This equation involves the effective sound speed squaredof the streaming plasma species, c α = c α γ α (1 − c α β α ) . (10)We have c α (cid:39) T α / (3 γ α ) and c α (cid:39) / (2 γ α ) in thenonrelativistic ( T α (cid:28) m e , (cid:98) Γ α (cid:39) /
3) and ultrarelativis-tic ( T α (cid:29) m e , (cid:98) Γ α (cid:39) /
3) thermal limits of a 3D gas,respectively. After integration, Eq. (9) provides a directrelationship between δp α and δA x and δN α , which canbe inserted in Eq. (5) to yield γ α w α ∂ t δ β ⊥ α + c α γ α w α n α ∇ ⊥ δN α = − q α γ α n α ( ∇ ⊥ δ Φ − β α ∇ ⊥ δA x ) . (11)Finally, combining this equation with Eq (6), one obtains ∂ t δN α − c α (cid:52) δN α = q α n α w α (cid:2) (cid:52) δ Φ − β α (cid:0) − c α (cid:1) (cid:52) δA x (cid:3) , (12)which determines the response of the apparent chargedensity to the electromagnetic perturbation: in Fourierspace, with δρ α = q α δN α , one finds δρ α = Ω α π k ω − c α k (cid:2) δ Φ − β α (cid:0) − c α (cid:1) δA x (cid:3) , (13)where the relativistic (proper-frame) plasma frequency,Ω p α , is defined by Ω α ≡ πn α q α w α /n α . (14) Note that the far-upstream plasma frequency Ω pp coin-cides with ω p defined earlier.The above also allows us to determine the response ofthe current density, δj αx , which comprises both a per-turbed conduction current density as well as a perturbedadvection current density: δj αx = ρ α δβ αx + β α δρ α . (15)Equations (4), (9) and (13) can then be combined toderive the response in Fourier space: δj αx = Ω α π (cid:26) β α (cid:0) − c α (cid:1) k ω − c α k δ Φ − (cid:20) − β α (cid:0) − c α (cid:1) ω − k ω − c α k (cid:21) δA x (cid:27) . (16)We define the “Weibel frame” as that in which thelinear instability becomes purely magnetic, i.e. , the elec-trostatic potential and the total electric charge densityvanish. From the response of the charge density, onefinds that in this frame the following relation must befulfilled: (cid:88) α Ω α β α | w − c α ζ − c α = 0 , (17)where ζ = ω/k is the (complex) phase velocity . In theshock precursor, the background plasma, of proper den-sity n p and nonrelativistic temperature T p (cid:46)
1, flows ata velocity β p | s (cid:39) − n b and relativistic temperature T b (cid:29) n p β p | w (cid:32) (cid:98) Γ p T p γ | w ζ (cid:33) + n b β b | w T b (cid:18) − (cid:98) Γ b (cid:19) (cid:39) , (18)where we have assumed | ζ | (cid:29) c α for both popula-tions, as is consistent within a hydrodynamic model (seeSec. III). This relation indicates that the velocity of R w is, in principle, mode-dependent. Yet the response of thebeam proves to be more sensitive to thermal effects thanthat of the plasma since the inequality ζ (cid:29) c α implies T p /γ | w ζ (cid:28)
1. As a consequence, β w | p ≡ − β p | w (cid:39) n b β b | w n p T b (cid:18) − (cid:98) Γ b (cid:19) (19)is a good (mode-independent) approximation of the ve-locity of R w relative to the background plasma frame.The beam is more conveniently characterized by itsnormalized pressure ξ b and temperature T b = κ T b γ ∞ m e ,so that n b = κ − T b γ ∞ ξ b n ∞ , see Sec. I B. Current con-servation further implies n ∞ /n p = γ p β p / ( γ ∞ β ∞ ), with β p (cid:39) β ∞ (cid:39) −
1, so that β w | p (cid:39) β b | w ξ b γ p γ ∞ (cid:98) Γ b − κ T b (cid:98) Γ b . (20)The beam moves relativistically with respect to the back-ground plasma, therefore either β b | w ∼ β p | w ∼ − ξ b <
1, however, the former must hold, which providesthe final result: β w | p (cid:39) ξ b γ p γ ∞ (cid:98) Γ b − κ T b (cid:98) Γ b . (21)In the case of negligible deceleration of the incomingplasma ( γ p (cid:39) γ ∞ ) and for a 3D adiabatic index (cid:98) Γ b = 4 / β w | p (cid:39) κ T b ξ b . (22)One can also directly calculate β w : β w = β w | p + β p β w | p β p (cid:39) β p (cid:32) − κ T b ξ b γ (cid:33) . (23)The above indicates that: (i) the “Weibel frame” R w moves at a sub-relativistic velocity relative to the back-ground plasma [ i.e. , u w | p = β w | p (1 − β | p ) − / < (cid:39) β p relative to theshock front; (ii) β w | p is of opposite sign to β p , whichimplies that, in magnitude, the “Weibel frame” movesslightly less fast than the background plasma relative tothe shock front; (iii) the relative velocity β w | p scales with ξ b .Finally, because ξ b is a function of distance x to theshock, Eq. (23) indicates that β w itself depends on x , i.e. ,the “Weibel frame” is not globally inertial, as alreadyobserved in Fig. 2. This observation has important im-plications with respect to the physics of the shock, whichare addressed in detail in the follow-up paper [52].In Sec. III, we extend the above calculations to a ki-netic description. Although the expression for β w | p willbe found to take somewhat different values, the abovethree features will remain valid.Using the above equations, it becomes straightforwardto derive the dispersion relation of the CFI in this warmfluid description. Consider the CFI in R w , where δ Φ | w =0. The response δj p x | w of the background plasma can bewritten δj p x | w (cid:39) ω π δA x | w (24)to first order in ξ b , since the term that has been neglectedhere with respect to Eq. (16) is of the order of β | w ∼ ξ .Note also that we have assumed that the backgroundplasma temperature remains sub-relativistic over mostof the precursor, as discussed in Paper II [52], so thatΩ pp (cid:39) ω p . From the Maxwell equations written in theLorentz gauge ( k µ δA µ = 0), i.e. , (cid:50) δA µ = − π (cid:80) α δj αµ ,one derives the dispersion relation in the “Weibel frame”(subscript w omitted for clarity): (cid:0) ω − k − ω (cid:1) = Ω (cid:20) − β | w (cid:0) − c (cid:1) ω − k ω − k c (cid:21) , (25) which, assuming | ω | (cid:28) k and c (cid:28)
1, can be ap-proximated by (cid:0) ω − k c (cid:1) (cid:0) ω − k − ω (cid:1) (cid:39) Ω β | w k , (26)with solution ω (cid:39) k c − Ω β | w k k + ω . (27)The stabilizing effect of the beam dispersion is manifestin this equation; this effect has been noted before in [48,49]. III. RELATIVISTIC KINETIC MODEL
In this Section, we evaluate the “Weibel frame” veloc-ity and the local instability growth rate within a rigorouskinetic formalism. As discussed further below, kinetic ef-fects must indeed be taken into account when consideringthe development of the CFI in the shock precursor. Thederivation of the kinetic dielectric tensor involves ratherheavy calculations, which are relegated to Apps. A, Band C. Below, we describe the general method, summa-rize the approximate expressions of the relevant dielectrictensor elements, and provide the estimates of the growthrate of the CFI and the associated estimate of β w in var-ious limits. The latter results, in particular, are given inSec. III C.For the time being, the reference frame in which wework is left unspecified. We first recall the linear disper-sion relation fulfilled by the CFI modes [56]: ε yy ( ε xx − /ζ ) = ε xy , (28)with ζ = ω/k as the complex phase velocity. In a fullyrelativistic framework, the dielectric tensor elements read( i, j = 1 , ,
3) [57] ε ij ( ω, k ) = δ ij + (cid:88) α γ α ω α ζ k (cid:90) u i γ ∂f (0) α ∂u j d u + (cid:88) α γ α ω α ζ k (cid:90) u i u j γ m e k∂f (0) α /∂u y ζ − v y d u , (29)where v y = u y / ( γm e ), ω α = 4 πn α e /m e representsthe nonrelativistic plasma frequency squared of species α ( n α represents as before the proper density) and f (0) α ( u )the corresponding unperturbed momentum distributionfunction, normalized such that (cid:82) d u f (0) α = 1. If the non-diagonal tensor element ε xy happens to vanish, Eq. (28)implies either ε yy = 0 or ζ ε xx − δ E (cid:107) (cid:98) y ) and purely electromagnetic (or in-ductive) modes (with δ E (cid:107) (cid:98) x ). Assuming that f (0) α ( u ) iseven in p y , ε xy reduces to ε xy ( ω, k ) = (cid:88) α γ α ω α ω (cid:90) u x u y m e γ ∂f (0) α /∂u y ζ − v y d u . (30)This expression is generally nonzero, and hence the CFIexcites mixed electromagnetic/electrostatic fluctuations.This feature has been often overlooked in the past, δE y =0 being assumed from the outset in a number of calcu-lations [58–61]. The electromagnetic ( δE x ) and electro-static ( δE y ) components of the solutions to Eq. (28) ver-ify [62, 63]: δE y δE x = − ε xy ε yy . (31)In Sect. II, β p | w was determined by solving ε xy = 0 inthe fluid limit, exploiting the fact that, to leading order,this equation is independent of the complex frequency ζ . In the general kinetic case, however, ε xy depends on ζ , whose knowledge involves solving Eq. (28). In prac-tice, we are interested in the fastest-growing mode ( ζ max ),which should be calculated in the (unknown) “Weibelframe”. Instead of solving simultaneously Eq. (28) and ε xy | w = 0 for ζ max | w and β p | w , we follow a different ap-proach, noting that the ( δE y , δB z ) fluctuations pertain-ing to a given mode in the plasma and “Weibel frames”are related through δE y | p = γ w | p β w | p δB z | w , (32) δB z | p = γ w | p δB z | w , (33)since δE y | w = 0. The velocity of the “Weibel frame”relative to the plasma rest frame is therefore given by β w | p = δE y | p δB z | p . (34)Making use of Eq. (31) and of δB z | p = − δE x | p /ζ max | p ,one obtains β w | p = ζ max | p ε xy | p ε yy | p . (35)This formula has the advantage of involving only quan-tities measured in the plasma frame.We apply the above formalism to the case of Maxwell-J¨uttner momentum distribution functions [64]: f (0) α ( u ) = µ α πm e γ α K ( µ α ) exp [ − γ α µ α ( γ − β α u x /m e )] , (36)where β α ≡ (cid:104) u x /γm e (cid:105) is the normalized mean drift veloc-ity of species α (corresponding drift Lorentz factor γ α ), µ α ≡ m e /T α and K denotes a modified Bessel functionof the second kind.Compact expressions of the tensor elements ε lm can beobtained in terms of one-dimensional integrals over thevelocity parallel to the wave vector ( v (cid:107) ) [63, 65]. Thesecalculations are detailed in Appendix A. In Sec. III D,such expressions will be used for the numerical resolu-tion of Eq. (35) along with Eq. (28). The parameters ofthe background plasma and suprathermal beam will bethen extracted from PIC simulations of relativistic colli-sionless shocks. In the remainder of this section, we will derive analytic approximations of ζ max | p and β w | p , validin distinct instability regimes for the plasma and beamparticles.The starting points of these calculations are the alter-native expressions (A30), (A31) and (A32) of the dielec-tric tensor. For instance, ε xx can be rewritten as (A30) ε xx = 1 + (cid:88) α ω α k ζ µ α γ α β α − ω α k ζ µ α γ α K ( µ α ) A αxx , (37)where A αxx = 2 µ α (cid:112) − ζ (cid:90) ∞−∞ d s χ α − s × (cid:26) γ α β α ( s + 1) / + 1 µ α (cid:20) s + 1 + 3 β α γ α ( s + 1) (cid:21) + 1 µ α (cid:20) s + 1) / + 3 β α γ α ( s + 1) / (cid:21)(cid:27) e − µ α √ s +1 , (38)and χ α = γ α ζ/ (cid:112) − ζ . The integrals involved in ε yy and ε xy can be put in a similar form, see (A34) and (A35).Introducing ∆ s α the characteristic width of the integrandin Eq. (38) [except for the denominator ( χ α − s ) − ], twolimiting cases can be considered for each plasma species: • the ‘hydrodynamic’ limit: | χ α | (cid:29) ∆ s α ; • the ‘kinetic’ limit: | χ α | (cid:28) ∆ s α .The dimensionless variable s is introduced immediatelybefore Eq. (A30); it corresponds to γ α β (cid:107) γ (cid:107) , with β (cid:107) thecomponent of the particle velocity along the wavevector.Since the wavevector is transverse to the streaming di-rection, the typical extent of s in the above integral is,up to the resonant factor, controlled by the proper tem-perature 1 /µ α ; the parameter χ α itself corresponds tothe apparent phase four-velocity of the mode. There-fore, the meaning of the hydrodynamic limit is that theapparent typical transverse momentum (normalized to m e ) exceeds the apparent phase four-velocity, while thekinetic limit corresponds to the opposite case. In the fol-lowing, approximate formulas of the dielectric tensor willbe derived in these two limits. A. Evaluation of the dielectric tensor for thebackground plasma
1. Hydrodynamic limit
In the outermost part of the precursor, the backgroundplasma is characterized by a nonrelativistic proper tem-perature, µ p (cid:29) s p (cid:39) (cid:112) /µ p ,and hence the hydrodynamic response of the backgroundplasma implies (cid:112) µ p / | χ p | (cid:29)
1. This condition coincideswith the large-argument limit ( ˜ χ p (cid:39) | χ p | (cid:112) µ p / (cid:29)
1) ofthe Z and Z (cid:48) functions involved in the low-temperatureexpressions derived in Appendix B. These formulas canbe further expanded to first order in 1 /µ p by making useof the asymptotic series Z ( η ) (cid:39) − /η − / η − / η · · · [66]: ε p xx (cid:39) − ω k ζ (cid:26) (cid:20) β + 1 µ p (cid:18) − β (cid:19) (cid:18) ζ − (cid:19)(cid:21)(cid:27) , (39) ε p yy (cid:39) − ω k ζ (cid:20) γ p µ p (cid:18) ζ − (cid:19)(cid:21) , (40) ε p xy (cid:39) − ω β p k ζ (cid:20) − µ p (cid:18) − γ (cid:19) (cid:0) − ζ (cid:1)(cid:21) . (41)In the rest frame of the background plasma, β p | p = 0;further assuming the weak-growth limit, | ζ | (cid:28)
1, theabove relations simplify to ε p xx (cid:39) − ω k ζ (cid:18) µ p ζ (cid:19) , (42) ε p yy (cid:39) − ω k ζ (cid:18) µ p ζ (cid:19) , (43) ε p xy = 0 . (44)
2. Kinetic limit
We now consider the limit ˜ χ p (cid:28) Z ( η ) (cid:39) i √ π exp( − η ) − η · · · [66] and assuming | ζ | (cid:28)
1, thisyields, in the background plasma rest frame ε p xx (cid:39) i (cid:114) πµ p ω k ζ , (45) ε p yy (cid:39) ω µ p k (cid:18) i (cid:114) πµ p ζ (cid:19) , (46) ε p xy = 0 . (47) B. Evaluation of the dielectric tensor for thesuprathermal particles
1. Hydrodynamic limit
In contrast to the background plasma, the beam parti-cles are characterized by an ultrarelativistic drift velocityin the background plasma rest frame ( γ b ≡ γ b | p (cid:29) µ b (cid:28) s b (cid:39)
1, so that the hydrodynamic response ofthe suprathermal particles requires | χ b | (cid:29)
1. The corre-sponding dielectric tensor, ε lm, b , is obtained by expand-ing ( χ b − s ) − (cid:39) χ − (cid:2) s/χ b ) (cid:3) in Eqs. (A33)-(A35),and evaluating the various resulting integrals. For ε xx, b , this gives A b xx (cid:39) µ γ b ζ (cid:20) γ β ∂ ∂b J (0 , µ b , − ∂∂b J (0 , µ b , (cid:21) + 4(1 − ζ ) µ γ ζ (cid:20) γ β ∂ ∂b I (0 , µ b , − ∂∂b I (0 , µ b , (cid:21) , (48)where the functions I ( t, λ, b ) and J ( t, λ, b ) are definedby Eqs. (C11) and (C4), respectively. Working out thederivatives, we find A b xx (cid:39) µ γ b ζ (cid:2) γ β µ b K ( µ b ) + K ( µ b ) (cid:3) + 4(1 − ζ ) µ γ ζ (cid:2) γ β µ b K ( µ b ) + K ( µ b ) (cid:3) . (49)Inserting this expression into Eq. (37) yields, to leadingorder, ε b xx (cid:39) − ω µ b β k ζ , (50)assuming | ζ | (cid:28) γ − (cid:28) µ b (cid:28)
1, and expandingthe Bessel functions accordingly.Applying the same procedure to Eqs. (A34) and (A35)leads to the hydrodynamic expressions ε b yy (cid:39) − ω µ b k ζ (cid:26) µ γ ζ (cid:20) − γ µ ln µ b (cid:21)(cid:27) , (51) ε b xy (cid:39) − ω µ b β b k ζ . (52)
2. Kinetic limit
The kinetic response of the beam particles can be read-ily obtained, to leading order in | χ b | , from the expansion( χ b − s ) − (cid:39) − iπδ ( s ) − P (1 /s ) in Eq. (A30), where δ ( s )is the Dirac delta function and P denotes the Cauchyprincipal value, which here vanishes. In general, how-ever, the beam particles appear to be only marginallykinetic in PIC shock simulations, so it could be useful togo to the next order. The series expansions derived inAppendix C are well suited to this purpose. In the limit | χ b | (cid:28)
1, Eqs. (C8), (C14) and (C17) reduce to A b xx (cid:39) − γ b (cid:26) i (cid:114) π µ b (cid:20) µ b K / ( µ b ) + β γ K / ( µ b ) (cid:21) − γ b ζ (cid:20) µ b K ( µ b ) + β γ K ( µ b ) (cid:21)(cid:27) , (53) A b yy (cid:39) − γ b ζµ b (cid:26) K ( µ b ) + 2 µ b K ( µ b )+ iγ b (cid:114) πµ b ζ (cid:20) K / ( µ b ) + 2 + β µ b K / ( µ b ) (cid:21)(cid:27) , (54) A b xy (cid:39) − γ b β b (cid:26) i (cid:114) π µ b K / ( µ b ) − ζ (cid:20) γ K ( µ b ) − µ b K ( µ b ) (cid:21)(cid:27) . (55)Combining these approximate expressions withEqs. (A30)-(A32) and expanding the Bessel functions inthe small-argument limit gives ε b xx (cid:39) ω µ b γ β k ζ + ω µ b γ k ζ [3 iπ − γ b ζ ] , (56) ε b yy (cid:39) ω µ b γ k (cid:20) iπ γ b ζ (cid:21) , (57) ε b xy (cid:39) ω µ b β b γ k ζ (cid:20) iπ γ b ζ − γ ζ (cid:21) . (58)where we have further assumed γ b (cid:28) C. CFI growth rates and frame velocities invarious plasma response regimes
The previous formulas may now be applied to thecase of the precursor of a relativistic shock to derive thegrowth rate of the purely transverse CFI and the veloc-ity of the “Weibel frame” R w . We consider the variouslimits in which the plasma and/or the beam can be de-scribed in a fluid-like or kinetic approximation, keepingin mind that the most relevant limit for the precursor isthat of both kinetic beam and background plasma.For reference, let us recall that in the limit | ζ | = | ω/k | (cid:28)
1, which is applicable here, the plasma can bedescribed in the hydrodynamic regime iff γ p | ζ | (cid:28) (cid:112) /µ p (with µ p = m/T p assumed greater than unity). As for thebeam, it can be described in the hydrodynamic regimeiff γ b | ζ | (cid:28) β w | p according to Eq. (35). We also assume that theplasma frame is close to the “Weibel frame”, so that theoff-diagonal term (cid:15) xy can be neglected in the dispersionrelation as written in the background plasma frame. Thedispersion relation may then be approximated as ζ ε xx − (cid:39) . (59) Finally, in order to make contact with our previ-ous notations, we will repeatedly use the substitution: ω µ b /ω = ξ b ( n ∞ /n p ) /κ T b . This notably implies ω µ b /ω (cid:28)
1. We also have γ b | p (cid:39) γ p | s and β b | p (cid:39)
1. Hydrodynamic plasma and beam
In the hydrodynamic regime (and in the backgroundplasma rest frame), | ζ | (cid:28) µ − / and | ζ | (cid:28) /γ b | p . Hence,the dispersion relation gives to leading order: ζ (cid:39) − ω µ b β | p k + 2 ω , (60)and so the growth rate saturates at Γ max (cid:39) ω µ b for k (cid:29) √ ω p .Adding up the hydrodynamic plasma and beam con-tributions into Eq. (35) and retaining only leading orderterms yields the “Weibel frame” velocity β w | p (cid:39) ω µ b β b | p − k ζ + 2 ω + ω µ b . (61)As | ζ | (cid:28) ω p /k according to Eq. (60), the expression for β w | p boils down to β w | p (cid:39) ω µ b β b | p ω (cid:39) κ T b ξ b n ∞ n p . (62)As β b | p (cid:39) β b | w (cid:39)
1, we recover the formula derivedwithin a fluid approach, Eq. (20), provided one sets inthe latter the adiabatic index (cid:98) Γ b = 2. This is the valueexpected for a gas with one degree of freedom in the rela-tivistic limit; the reduced effective dimensionality for thebeam response results from the assumption of a purely1D transverse fluctuation spectrum.Finally, the fully hydrodynamic regime holds as long asmin( (cid:112) µ p / , γ b | p ) | ζ max | (cid:29)
1. Now, expressing Eq. (60)at k max ≡ √ ω p gives ζ max (cid:39) i (cid:32) ω µ b ω (cid:33) / , (63)so that another way of expressing the validity of the hy-drodynamic regime is:min (cid:16) µ p , γ | p (cid:17) > ω ω µ b (cid:39) κ T b ξ b n p n ∞ . (64)
2. Kinetic plasma and hydrodynamic beam
In this limit, 1 /γ b | p (cid:29) | ζ | (cid:29) µ − / ; to leading order,the dispersion relation takes the form i (cid:112) πµ p ω ζ − ω µ b β | p ζ − k (cid:39) . (65)0The growth rate reaches its maximum value Γ max = (cid:113) ω µ b for k (cid:46) k max (cid:39) (2 πω µ p µ b ) / ω / . We de-fine | ζ max | = Γ max /k max (cid:39) (cid:16) ω µ b / (cid:112) πµ p ω (cid:17) / .Making use of Eqs. (46), (47), (51), and (52), we ex-press Eq. (35) as β w | p (cid:39) − ω µ b β b ω µ p ζ , (66)which gives β w | p (cid:39) (cid:32) π ω µ b ω µ (cid:33) / (cid:39) (cid:32) π κ T b (cid:33) / ξ / (cid:32) n ∞ n p T m (cid:33) / . (67)
3. Hydrodynamic plasma and kinetic beam
In this limit, we now have 1 /γ b | p (cid:28) | ζ | (cid:28) µ − / . UsingEqs. (42) and (56), the dispersion relation reduces to2 ω µ b γ | p (cid:20) i π γ b | p ζ (cid:21) − ω (cid:20) µ p ζ (cid:21) − k (cid:39) . (68)Let us first assume that ω / ( µ p ζ ) can be neglected infront of 2 ω µ b γ | p ζ . An unstable solution then existsprovided γ | p > ω / (cid:16) ω µ b (cid:17) ∼ /ξ b . This is the samecondition as encountered in the 4 − fluid system of Sec. II,see Eq. (27) where c , b ∼ /γ | p and Ω ∼ ω µ b .If this condition is fulfilled, the fastest-growing solu-tion corresponds to Γ max (cid:39) π (cid:113)
23 [ ω µ b γ | p − ω ] / ω µ b γ | p and ζ max (cid:39) i πγ b | p (cid:18) − ω ω µ b γ | p (cid:19) . Note that the condi-tions for a hydrodynamic plasma ( (cid:113) µ p | ζ max | >
1) anda kinetic beam ( γ b | p | ζ max | <
1) are then verified, albeitmarginally.Moreover, combining Eqs. (43), (44), (57), and (58)gives the “Weibel frame” velocity β w | p (cid:39) ω µ b γ | p ζ k ζ − ω + 2 ω µ b γ | p ζ . (69)Inserting the above expressions of Γ max and ζ max , therefollows β w | p (cid:39) (cid:18) π (cid:19) ω µ b β b | p ω (cid:32) − ω ω µ b γ | p (cid:33) (70) (cid:39) (cid:18) πκ T b (cid:19) ξ b n ∞ n p (cid:32) − n p n ∞ κ T b γ | p ξ b (cid:33) . (71)In the opposite limit, in which ω / ( µ p ζ ) domi-nates over 2 ω µ b γ | p ζ in Eq. (68), the dominant mode satisfies Γ max (cid:39) (cid:113) ω /µ p and ζ max (cid:39) i [2 µ p (1 − ω µ b γ | p /ω )] − / for k max ≡ (cid:113) ω − ω µ b γ | p . Toleading order, we thus derive β w | p (cid:39) ω µ b γ | p ω µ p (cid:16) − ω µ b γ | p /ω (cid:17) (72) (cid:39) κ T b γ | p ξ b T p m n ∞ n p (cid:32) − κ T b γ | p ξ b n ∞ n p (cid:33) − (73)
4. Kinetic plasma and beam
Finally, we consider the case of a fully kinetic beam-plasma system. This regime is of particular importancesince it is found to hold in most of the precursor regionin long-time shock simulations (see Sec. III D). Using theexpressions (45) and (56), the dispersion relation writes i (cid:112) πµ p ω ζ + 2 ω µ b γ | p (cid:18) i π γ b | p ζ (cid:19) − k (cid:39) . (74)The dominant CFI mode is then defined byΓ max (cid:39) (cid:114)
23 ( ω µ b ) / γ | p β | p (cid:112) πµ p ω + π ω µ b γ | p , (75) ζ max (cid:39) i ω µ b γ | p (cid:112) πµ p ω + π ω µ b γ | p , (76)and k max (cid:39) (cid:113) µ b ω pb γ b | p .The corresponding expression for the “Weibel frame”velocity is obtained by combining Eqs. (46), (47), (57),and (58). After some algebra, one finds β w | p (cid:39) ω µ b γ b | p β b | p i π γ b | p ζ max − γ | p ζ ω µ p + ω µ b γ | p (77)or, to leading order and in terms of our usual parameters, β w | p (cid:39) κ T b ξ b n ∞ n p γ | p T p m (78) D. Comparison to PIC simulations
In this Section, we compare the relative velocity β w | p estimated using the kinetic model of the CFI developedin Sec. III C with that extracted from 2D3V PIC simula-tions.These simulations have been performed using the mas-sively parallel, relativistic PIC code calder [67]. Theshock is generated by means of the standard mirror tech-nique [40]. The background pair plasma is continuously1injected into the domain from the right-hand boundary,and is made to reflect specularly off the left-hand bound-ary ( x | d = 0). Electrons and positrons are injected witha Maxwell-J¨uttner momentum distribution [Eq. (36)] ofproper temperature T ∞ = 0 . m e and mean drift veloc-ity β ∞| d = − (1 − /γ ∞| d ) / in the simulation frame. Asmentioned earlier, this simulation reference frame coin-cides with the downstream plasma rest frame (as a conse-quence of the use of the mirror technique). Our two ref-erence PIC simulations use γ ∞| d = 10 and γ ∞| d = 100,which respectively correspond to shock Lorentz factors(with respect to the upstream) of γ ∞ = 17 and γ ∞ = 173.As the simulation proceeds, the right-hand boundary(injector) is progressively displaced towards the right soas to keep the reflected ballistic particles inside the (in-creasingly large) domain, while speeding up the calcu-lation at early times [40]. In order to quench the nu-merical Cherenkov instability, which notoriously ham-pers simulations of relativistically drifting plasmas, wemake use of the Godfrey-Vay filtering scheme, com-bined with the Cole-Karkkainnen finite difference fieldsolver [68]. We use a mesh size ∆ x = ∆ y = 0 . c/ω p and a time step ∆ t = 0 . x/c . Periodic bound-ary conditions are employed in the tranverse directionfor both particles and fields. The initial domain sizeis L x × L y = 2700 ×
340 ( c/ω p ) for γ ∞| d = 10 and L x × L y = 2000 ×
200 ( c/ω p ) for γ ∞| d = 100. Eachcell is initially filled with 10 macro-particles per species(electrons or positrons). The simulation is run until t max = 3600 ω − (resp. t max = 6900 ω − ) for our sim-ulation with γ ∞| d = 10 (resp. γ ∞| d = 100).In order to test our model of the CFI developing in theprecursor through the interpenetration of suprathermaland background plasma populations, we need to care-fully distinguish these two in the simulations. In orderto do so, we track the particles according to the signof their x − momentum and how many times this signhas changed, due only to interactions with the electro-magnetic turbulence. We then define the backgroundplasma as those particles that move toward negative val-ues of x and that do not have undergone turnarounds, i.e. , any change of sign of u x . We define the beam par-ticles as those that move towards positive values of x ,independently of their number of turnarounds. This def-inition leaves a minority of particles: those that movewith u x < ξ b (cid:39) / / F ∞ .We extract hydrodynamic moments n α , T α , u α foreach of the beam and background plasma population,assuming that they obey Maxwell-J¨uttner momentumdistributions. The spatial profiles of these various hy-drodynamic quantities have already been presented inFig. 1, as extracted from the simulations with γ ∞| d = 10and γ ∞| d = 100 at respective times t (cid:39) ω − and t (cid:39) ω − . One can see that n p and γ p | d vary weaklyacross the precursor region, except near the shock frontwhere the incoming plasma slows down significantly andexperiences compression. By contrast, the plasma tem-perature steadily increases from its far-upstream value( T p = 0 . m e ) to unity and beyond when approachingthe shock front. This heating results from the interactionwith the beam particles, whose density rises by ∼ − max ), wave number ( k max ) and wave phase veloc-ity ( ξ max ) of the fastest-growing mode, and then useEq. (35) to evaluate the corresponding value of the“Weibel frame” velocity ( β w | p ).Figures 3 display the spatial profiles of the χ b and˜ χ p parameters defined by Eqs. (A36) and (B5), respec-tively. Interestingly, both the background plasma andbeam populations appear to lie in the kinetic CFI regime( ˜ χ p < χ b <
1) throughout the precursor region.However, whereas ˜ χ p shows weak variations around rel-atively low values ( ˜ χ p ∼ . χ b values are larger byabout an order of magnitude and present stronger varia-tions. Where ξ b becomes of the order of unity, the beamresponse is then only marginally kinetic; consequently,analytical approximations present an error of about afactor 2 with respect to the full numerical calculation of β w | p .As discussed in Sec. I B, we estimate the 3-velocityof the “Weibel frame” in PIC simulations as the ra-tio β w | d = (cid:10) δE y (cid:11) / / (cid:10) δB z (cid:11) / , where averaging is per-formed over the transverse dimension. Given the lo-cally measured value of the background plasma ve-locity β p | d , we convert it to the instantaneous lo-cal plasma rest frame through the standard transform β w | p = (cid:0) β w | d − β p | d (cid:1) / (cid:0) − β w | d β p | d (cid:1) . We then comparein Fig. 4 this measurement with our theoretical estimatesof the relative velocity, β w | p , between the “Weibel frame”and the background plasma, using both the numericalcalculation (red circle/dashed line) and the analytical ap-proximation (purple square/dashed). Both for γ ∞ = 17and γ ∞ = 173, it is seen that, as predicted, β w | p re-mains subrelativistic throughout the precursor, increas-2 - - - - FIG. 3. Parameters χ b and ˜ χ p of the beam and the back-ground plasma as a function of distance to the shock front,extracted from PIC simulations with γ ∞ = 17 [panel (a), top]and γ ∞ = 173 [panel (b), bottom]. We recall that the hy-drodynamic (kinetic) regime for the beam and/or the plasmacorresponds to χ b (cid:29) χ b (cid:28)
1) and/or ˜ χ p (cid:29) χ p (cid:28) ing from β w | p ∼ − at the far end of the precursor upto β w | p (cid:39) . β w | p can be measured accurately.For both reference simulations, Fig. 4 reveals signif-icant fluctuations in the measured values of β w | p , withincreasing amplitude at large x , which deserves some dis-cussion. Far from the shock front, the magnitude of ourobservable (cid:10) δE y (cid:11) / / (cid:10) δB z (cid:11) / is close to unity (see e.g. ,Fig. 2), which implies that any amount of numerical noisecan artificially push | β w | d | to values larger than unity,even though its true value might be <
1. Furthermore,when transforming values to the background plasma restframe, any error in β w | d is amplified by ∼ γ | d , which is ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ● ● ◼ ◼ - - - ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ● ● ◼ ◼ - - - FIG. 4. Theoretical estimates of the relative velocity betweenthe “Weibel frame” R w and the background plasma, β w | p , asa function of distance to the shock front, compared with thevelocity extracted from our reference PIC simulations throughthe ratio (cid:10) δE y (cid:11) / / (cid:10) δB z (cid:11) / (solid green); where this valuecannot be measured accurately from the simulation, the dataare colored in light green (see text for details). The red circlesymbols/dashed curve uses the numerical solution to the gen-eral kinetic dispersion relation (28) to derive ζ max , while thepurple square/dashed curve plots the analytic approximationsderived in Sec. III C. Panel (a), top: γ ∞ = 17 ( γ ∞| d = 10);panel (b), bottom: γ ∞ = 173 ( γ ∞| d = 100); large far from the shock.In a given portion of the precursor, the “Weibel frame”can be considered well defined where | β w | d | < | β w | p | < | β w | d | <
1, but the use of raw numerical data,binned linearly, and plotted on a log-scale, precludes aclear identification of the region where this frame is well-defined, at least by eye. By rebinning the data, however,we infer that the “Weibel frame” is well defined in therange x (cid:46) c/ω p , for both PIC simulations, as claimedin Sec. I B. At larger distances from the shock, | β w | p | fluc-3tuates too often on both sides of unity, either because ofnumerical noise, or because of the contribution of electro-static modes. Given the magnitude of the fluctuations,we have chosen to plot the data in full color only in theregion x ≤ c/ω p , for the sake of being conservative. IV. THE “WEIBEL FRAME” IN THENONLINEAR REGIMEA. Theoretical model
We adopt here another perspective on the “Weibelframe”. Specifically, we consider the nonlinear evolutionof the CFI, once the current filaments have formed, bor-rowing on the work of Ref. [31]. We assume that, at eachpoint in the shock precursor, the CFI has developed aquasistatic, transversally periodic system of current fil-aments. By quasistatic, it is meant that an equilibriumapproximately holds between magnetic and thermal pres-sures in the filaments, according to the physical condi-tions at the point considered, and that these conditionsevolve slowly enough that this equilibrium has time toadapt from one point to another. We also consider a2D configuration, without loss of generality, so that theplasma is periodic along the y -axis and the species driftalong the x -axis; the magnetic field is then transverse tothe ( x, y )-plane. The dominant components of the four-potential now explicitly depend on y and are given by A µ = [ φ ( y ) , A x ( y ) , , n α = n α exp (cid:20) − γ α q α T α ( φ − β α A x ) (cid:21) , (79)where q α denotes the charge of the species and n α repre-sents a normalization prefactor. We can inject these den-sity profiles in the potential formulation of the Amp`ere-Maxwell and Gauss-Maxwell equations, leading to ∂ y A x = (cid:88) α ω α c γ α β α sinh (cid:20) γ α eT α ( φ − β α A x ) (cid:21) m e c e , (80) ∂ y φ = (cid:88) α ω α c γ α sinh (cid:20) γ α eT α ( φ − β α A x ) (cid:21) m e c e , (81)where ω α = (cid:112) πn α e /m e scales to the plasma fre-quency of species α . It is worth noting that this systemcan also be obtained in a 4-fluid framework with isother-mal closure condition [31].Following [31], we introduce the plasma nonlinearityparameter: Ξ p = (cid:12)(cid:12)(cid:12)(cid:12) γ p β p eT p max y A x ( y ) (cid:12)(cid:12)(cid:12)(cid:12) . (82) In the weakly nonlinear limit, Ξ p (cid:28)
1, the sinh functionsin the above equations can be approximated to unity, sothat the vanishing of the electrostatic component entails: n b γ | w β b | w T b + n p γ | w β p | w T p = 0 (83)where the subscript w has been introduced, because theabove quantities are now defined in the “Weibel frame”, R w , in which φ | w = 0.In the weakly nonlinear limit, the velocity of R w canbe computed exactly using relation (83). Writing β b | w =( β b − β w ) / (1 − β b β w ), γ b | w = γ b γ w (1 − β b β w ) etc., in anygiven frame, one finds that β w is solution to the equation: β − Q w β w + 1 = 0 (84)with Q w = n b T b γ (cid:0) β (cid:1) + n p T p γ (cid:0) β (cid:1) n b T b γ β b + n p T p γ β p . (85)Writing n b in terms of ξ b as before, we expand the abovesolution to first order in ξ b to obtain the relative velocity β w | p : β w | p (cid:39) + γ | p ξ b T p m e n ∞ n p β ∞ β b | p κ T b (86)where γ b | p ∼ γ p represents the relative Lorentz factorbetween the beam and the background plasma. Inter-estingly, Eq. (86) corresponds to our earlier expressionEq. (78) obtained from the linear dispersion relation ofthe CFI in the kinetic plasma – kinetic beam limit.In Paper II [52], it is argued that γ ξ b is much smallerthan unity in the far precursor, where the incomingbackground plasma maintains its initial velocity, i.e. , γ p (cid:39) γ ∞ , and of the order of unity in the near pre-cursor, where γ p < γ ∞ due to deceleration. The abovethus implies that the “Weibel frame”, in this nonlineardescription, moves at subrelativistic velocities relative tothe background plasma, as in the linear limit studied ear-lier. That β w | p is positive means that the “Weibel frame”moves at slightly smaller absolute velocity towards theshock front than the background plasma.In the near precursor, as the background plasma isheated up to relativistic temperatures, β w | p increases inmagnitude; this implies that the background plasma de-couples from the “Weibel frame”, hence increasing theheating rate and leading to the shock transition.One can also compute the first nonlinear correction inΞ p to the above velocity. To this effect, we recast Eq. (81)in terms of the nonlinearity parameters of the beam (Ξ b )and the plasma (Ξ p ) in relation (83): n b γ | w β b | w T b sinh Ξ b | w Ξ b | w + n p γ | w β p | w T p sinh Ξ p | w Ξ p | w = 0 . (87)4As discussed in Paper III [53], the beam particles carrysuch inertia that they hardly participate in the filamen-tation, meaning Ξ b | w (cid:28)
1. In R w , Ξ b | w indeed repre-sents the ratio of the electromagnetic component eA x | w to the typical momentum T b | w = T b /γ b | w of the parti-cles, and for suprathermal particles, this ratio is muchsmaller than unity. Assuming Ξ p | w (cid:46)
1, we then obtainto lowest order n b γ | w β b | w T b + n p γ | w β p | w T p (cid:32) | w (cid:33) (cid:39) . (88)In this configuration, the nonlinearity appears as asecond-order correction to the “Weibel frame” speed. Indetail, one obtains β (n − lin)w | p (cid:39) β w | p (cid:32) − Ξ (cid:33) . (89)Note that the “Weibel frame” is not always well-defined in the strongly nonlinear limit: if Ξ p , b (cid:29) γ b | w n b = γ p | w n p , (90) γ b | w β b | w T b = − γ p | w β p | w T p , (91)which overdetermine the system for a given set of param-eters. From Eq. (81), we see that the error on the elec-trostatic fields at leading order in Ξ p | w if we evaluate the“Weibel frame” from relation (83) evolves as δφ ∝ Ξ | w .In the case of interest, however, PIC simulations indi-cate that the weakly nonlinear limit represents a goodapproximation in the precursor of relativistic shocks, seeSec. IV B. B. Comparison to PIC simulations
In this section, we confront the result of Sec. (IV) withthe PIC simulations presented in the Sec. III D. A stronghypothesis underlying the above formulae is that of amarginally nonlinear plasma response Ξ p | w <
1. We wishto motivate this hypothesis by measuring this nonlinear-ity parameter. Assuming Eq. (79) holds, one expects n p ( y ) = n p0 exp (cid:20) − γ p | s T p (cid:0) φ − β p | s A x (cid:1)(cid:21) = n p0 exp (cid:20) γ p | w β p | w T p A x (cid:21) . (92)We can then estimate the nonlinearity parameter in thesimulation through the following relationΞ p | w (cid:39) √ (cid:42)(cid:20) log (cid:18) n p (cid:104) n p (cid:105) (cid:19) − (cid:28) log (cid:18) n p (cid:104) n p (cid:105) (cid:19)(cid:29)(cid:21) (cid:43) / , (93) - - - - FIG. 5. Evolution of the nonlinearity parameter of the back-ground plasma Ξ p | w . Panel (a), top: γ ∞| d = 10 at time t = 3600 ω − . Panel (b), bottom: γ ∞| d = 100 at time t = 6900 ω − . where (cid:104)·(cid:105) denotes the mean value along the directiontransverse to the drift.Using relation (93), we present in Figs. 5 our estimateof Ξ p | w for the two reference simulations ( γ ∞ = 17 in thetop panel, γ ∞ = 173 in the bottom panel). In both cases,the nonlinearity parameter tends to increase from valueswell below unity in the far precursor, up to near unitywithin a few hundred c/ω p to the shock front. This be-havior suggests that one should use the present nonlinearequilibrium model to estimate β w | p in the near precur-sor, say x (cid:46) c/ω p , and the linear estimate obtainedin Sec. III at larger distances, where Ξ p | w falls to val-ues small compared to unity. One must, however, keepin mind that the present approach is based on a fluidmodel, while the former is fully kinetic, and that theanalysis of Sec. III indicates that, from the point of viewof the instability, both the beam and the backgroundplasma should be treated kinetically. Fortunately, boththe present nonlinear estimate of β w | p , Eq. (86), and its5 ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆◆ ◆ - - - ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆◆ ◆ - - - FIG. 6. Similar to Fig. 4, for our theoretical estimate of β w | p obtained assuming nonlinear equilibrium filaments, as devel-oped in Sec. IV. The purple diamond symbols/dashed curveuses (86) to compute the β w | p . Panel (a), top: γ ∞ = 17( γ ∞| d = 10). Panel (b), bottom: γ ∞ = 173 ( γ ∞| d = 100). linear kinetic counterpart, Eq. (78), match one another.Hence, our theoretical estimate of β w | p can be consideredas rather well established throughout the precursor.This is supported by Figs. 6, which compare the valuesof β w | p extracted from the simulations with the theoret-ical estimate (86), assuming (weakly) nonlinear equilib-rium filaments. As before, the PIC simulation data arelight colored where they cannot be measured accurately.In both simulation cases, the theoretical estimates appearin reasonable agreement with the PIC results, especiallyin the near precursor x (cid:46) c/ω p , where the velocity ofthe “Weibel frame” is well defined. V. CONCLUSIONS
This paper belongs to a series of articles in which webuild a theoretical model of unmagnetized, relativisticcollisionless pair shocks, and compare it with dedicatedPIC simulations. More specifically, we have discussedin the present work the physics of the purely transverseCFI that results from the interpenetration of the back-ground plasma and the beam of suprathermal particlesin the shock precursor. We have argued that there ex-ists a frame R w , referred to as the “Weibel frame” inwhich the instability is mostly magnetic in nature. Wehave derived the velocity of this frame at each point ofthe shock precursor, through a kinetic description of thelinear stage of the instability, as well as through a qua-sistatic model of the nonlinear phase of the filamentationinstability. This “Weibel frame” is of particular impor-tance to the physics of the shock and of the accelerationprocess, because it represents the frame of the scatteringcenters.We emphasize the following properties of this “Weibelframe”: (i) it is found to move at subrelativistic velocitiesrelative to the background plasma, ( i.e. , with 4-velocity u w | p < β w | p is of oppo-site sign to β p , the velocity of the background plasma inthe shock rest frame, implying that the “Weibel frame”moves slightly less fast than the background plasma rel-ative to the shock front. Furthermore, (iii), the relativevelocity β w | p typically scales as ξ b , which represents thebeam pressure normalized to the incoming momentumflux at infinity. Finally, it is to be emphasized that the“Weibel frame” is not globally inertial, because its veloc-ity β w depends on the distance to the shock. As discussedin [51], and more particularly [52], which addresses indetail the consequences of this noninertial nature, the“Weibel frame” decelerates from the far to the near pre-cursor. In the present article, we have determined the ve-locity of this frame at each point in the precursor, giventhe local physical conditions, in a WKB-like approxima-tion.Our PIC simulations confirm these various features. Inparticular, they reveal that δE ⊥ < δB ⊥ close to the shockfront and over most of the precursor, and hence that the“Weibel frame” is well defined; relatively to the shockfront, the “Weibel frame” velocity, measured through theratio of δE ⊥ /δB ⊥ , is also found to be slightly below thebackground plasma velocity; finally, the relative velocity β w | p generally decreases away from the shock front, like ξ b .Our kinetic model relies on the use of Maxwell-J¨uttnerdistribution functions for the beam and the backgroundplasma, and it provides (at the expense of rather complexcalculations) simple approximations to the velocity of the“Weibel frame” and to the growth rate of the purely6transverse CFI. Our quasistatic model of the nonlinearstage of the instability describes the current filaments asperiodic magnetostatic structures in the “Weibel frame”,in pressure equilibrium with the plasma. Our PIC sim-ulations indicate that, over most of the precursor, thosefilaments are in a mildly nonlinear stage, with a nonlin-earity parameter below unity. We expect that both mod-els should capture the salient features of the instabilityand indeed, the resulting formulae turn out to bracketrather well the value β w | p seen in PIC simulations.In subsequent papers, it will be shown that the “Weibelframe” plays an essential role in shaping the microphysicsof the shock transition, in particular the physics of heat-ing and deceleration of the background plasma. ACKNOWLEDGMENTS
We acknowledge financial support from the Pro-gramme National Hautes ´Energies (PNHE) of the C.N.R.S., the ANR-14-CE33-0019 MACH project andthe ILP LABEX (under reference ANR-10-LABX-63)as part of the Idex SUPER, and which is financed byFrench state funds managed by the ANR within the ”In-vestissements d’Avenir” program under reference ANR-11-IDEX-0004-02. This work was granted access to theHPC resources of TGCC/CCRT under the allocation2018-A0030407666 made by GENCI. We also acknowl-edge PRACE for awarding us access to resource JoliotCurie-SKL at TGCC Center.
Appendix A: Reduced expressions of the kinetic dielectric tensor
The triple integrals involved into Eq. (29) after substitution of the model distribution (36) can be recast in the formof one-dimensional quadratures [63]. To do so, we mostly follow the lines sketched in Ref. [65], where it is shown thattwo of the momentum integrations can be carried out in closed form.Let us first consider the dielectric tensor element ε xx , which is involved in the purely electromagnetic dispersionrelation of the CFI ( i.e. , in a frame close to the “Weibel frame”). After straightforward algebra, it can be rewrittenas ε xx = 1 + (cid:88) α ω α ω C α γ α µ α β α (cid:90) d u v x e − γ α µ α ( γ − β α u x ) − ω α ζk C α γ α µ α (cid:90) d u v x e − γ α µ α ( γ − β α u x ) ζ − v y , (A1)where we have defined C α = µ α πγ α K ( µ α ) . (A2)The first integral in the right-hand side of (A1) can be exactly solved by noting that β α = (cid:90) d u v x e − γ α µ α ( γ − β α u x ) . (A3)We then proceed by changing to velocity variables in cylindrical coordinates along the wave vector: v =( v (cid:107) , v ⊥ cos θ, v ⊥ sin θ ). Making use of d u = γ d v , Eq. (A1) then becomes ε xx = 1 + (cid:88) α ω α ω γ α β α µ α − ω α ζk C α γ α µ α (cid:90) d v (cid:107) ζ − v (cid:107) (cid:90) /γ (cid:107) d v ⊥ v ⊥ γ e − γ α µ α γ (cid:90) π d θ cos θ e γ α β α µ α γv ⊥ cos θ , (A4)with γ (cid:107) = (1 − v (cid:107) ) − / . Given the integral representation of the modified Bessel functions of the first kind [70], I n ( z ) = 1 π (cid:90) π d t e z cos t cos( nt ) , (A5)we obtain ε xx = 1 + (cid:88) α ω α ω C α µ α γ α β α − πω α k ζ C α γ α µ α (cid:90) − d v (cid:107) S α + S α ζ − v (cid:107) . (A6)7Here we have introduced S αmn ( v (cid:107) ) = (cid:90) /γ (cid:107) d v ⊥ γ v m ⊥ I n ( µ α γ α β α v ⊥ ) e − γ α µ α γ . (A7)From the definition t = γ/γ (cid:107) follow the relations γv ⊥ = ( t − / and v ⊥ d v ⊥ = d t/γ (cid:107) t , allowing us to express S α and S α as S α ( v (cid:107) ) = γ (cid:107) (cid:90) ∞ d t ( t − I (cid:16) µ α γ α β α (cid:112) t − (cid:17) e − γ α µ α γ (cid:107) t , (A8) S α ( v (cid:107) ) = γ (cid:107) (cid:90) ∞ d t ( t − I (cid:16) µ α γ α β α (cid:112) t − (cid:17) e − γ α µ α γ (cid:107) t , (A9)which can be further recast as S α ( v (cid:107) ) = γ (cid:107) (cid:32) γ α µ α d dγ (cid:107) − (cid:33) (cid:90) ∞ d t I (cid:16) µ α γ α β α (cid:112) t − (cid:17) e − γ α µ α γ (cid:107) t , (A10) S α ( v (cid:107) ) = γ (cid:107) (cid:32) γ α µ α d dγ (cid:107) − µ α ddγ (cid:107) + 1 (cid:33) (cid:90) ∞ d t I (cid:16) µ α γ α β α (cid:112) t − (cid:17) e − γ α µ α γ (cid:107) t , (A11)We now take advantage of the following formula [71] (cid:90) ∞ d x (cid:18) x − x + 1 (cid:19) ν/ e − δx I ν (cid:16) ζ (cid:112) x − (cid:17) = e − √ δ − ζ (cid:112) δ − ζ (cid:32) ζδ + (cid:112) δ − ζ (cid:33) ν , (A12)valid for (cid:60) ν > − δ > ζ . There follows S α = γ (cid:107) γ α µ α (cid:32) γ α µ α d dγ (cid:107) − (cid:33) (cid:20) e − γ α µ α Π α Π α (cid:21) , (A13) S α = β α γ (cid:107) γ α µ α (cid:32) γ α µ α d dγ (cid:107) − γ α µ α ddγ (cid:107) + 1 (cid:33) (cid:20) e − γ α µ α Π α Π α ( γ (cid:107) + Π α ) (cid:21) . (A14)where Π α = ( γ (cid:107) − β α ) / . After evaluation of the first and second-order derivatives, one finds S α + S α = 2 γ (cid:107) γ α µ α (cid:34) γ (cid:107) + 2 β α Π α + γ α µ α ( γ (cid:107) + 2 β α )Π α + γ α µ α β α Π α (cid:35) e − γ α µ α Π α . (A15)Substituting Eq. (A15) into Eq. (A6) finally gives the following simplified expression for ε xx : ε xx = 1 + (cid:88) α ω α ω µ α γ α β α − πω α k ζ C α γ α (cid:90) − d v (cid:107) γ (cid:107) ζ − v (cid:107) (cid:34) γ (cid:107) + 2 β α γ α µ α Π α + γ (cid:107) + 2 β α γ α µ α Π α + β α Π α (cid:35) e − γ α µ α Π α . (A16)Let us now consider ε yy , which writes ε yy = 1 − ω α ζk (cid:88) α γ α µ α (cid:90) d u v y e − γ α µ α ( γ − β α u x ) ζ − v y . (A17)After changing to cylindrical velocity coordinates, one obtains ε yy = 1 − (cid:88) α ω α ζk C α γ α µ α (cid:90) − d v (cid:107) v (cid:107) ζ − v (cid:107) (cid:90) /γ (cid:107) d v ⊥ γ v ⊥ e − γ α µ α γ (cid:90) π d θ e γ α µ α β α γv ⊥ cos θ , (A18)which, making use of Eq. (A5), reduces to ε yy = 1 − (cid:88) α πω α ζk C α γ α µ α (cid:90) − d v (cid:107) v (cid:107) S α ζ − v (cid:107) . (A19)8Making use of Eq. (A12), one obtains S α ( v (cid:107) ) = γ (cid:107) γ α µ α d d γ (cid:107) (cid:90) ∞ d t I (cid:16) µ α γ α β α (cid:112) t − (cid:17) e − γ α µ α γ (cid:107) t = (cid:18) γ (cid:107) γ α µ α (cid:19) d d γ (cid:107) (cid:20) e − γ α µ α Π α Π α (cid:21) , (A20)and therefore S α ( v (cid:107) ) = (cid:18) γ (cid:107) γ α µ α (cid:19) (cid:34) γ (cid:107) + β α Π α + γ α µ α (2 γ (cid:107) + β α )Π α + γ α µ α γ (cid:107) Π α (cid:35) e − γ α µ α Π α . (A21)Combining Eqs. (A19) and (A21) then readily yields ε yy = 1 − (cid:88) α πω α ζk C α γ α (cid:90) − d v (cid:107) v (cid:107) γ (cid:107) ζ − v (cid:107) (cid:34) γ (cid:107) + β α γ α µ α Π α + 2 γ (cid:107) + β α γ α µ α Π α + γ (cid:107) Π α (cid:35) e − γ α µ α Π α . (A22)Finally, let us rewrite the non-diagonal dielectric tensor element ε xy as ε xy = (cid:88) α ω α k ζ µ α γ α β α − ω α k C α γ α µ α (cid:90) d u v x e − γ α µ α ( γ − β α u x ) ζ − v y , (A23)Changing to cylindrical velocity variables leads to ε xy = (cid:88) α ω α k ζ µ α γ α β α − πω α k C α γ α µ α (cid:90) − d v (cid:107) ζ − v (cid:107) (cid:90) /γ (cid:107) d v ⊥ γ v ⊥ e − γ α µ α γ (cid:90) π d θ e γ α µ α β α γ α v ⊥ cos θ . (A24)Using Eq. (A5) gives ε xy = (cid:88) α ω α k ζ µ α γ α β α − πω α k C α γ α µ α (cid:90) − d v (cid:107) S α ζ − v (cid:107) , (A25)where S α can be recast in the form S α ( v (cid:107) ) = γ (cid:107) (cid:90) ∞ d t (cid:18) t − t + 1 (cid:19) / t ( t + 1) I (cid:16) γ α µ α β α (cid:112) t − (cid:17) e − γ α µ α γ (cid:107) t . (A26)Exploiting Eq. (A12) yields S α ( v (cid:107) ) = γ (cid:107) β α γ α µ α (cid:32) d dγ (cid:107) − γ α µ α dd γ (cid:107) (cid:33) (cid:20)(cid:18) γ (cid:107) Π α − (cid:19) e − γ α µ α Π α (cid:21) , (A27)After evaluation of the derivatives, one finds S α ( v (cid:107) ) = β α (cid:18) γ (cid:107) γ α µ α (cid:19) (cid:20) γ α µ α Π α + 3 γ α µ α Π α + 3Π α (cid:21) e − γ α µ α Π α , (A28)Combining Eqs. (A25) and (A28), one obtains the simplified expression ε xy = (cid:88) α ω α k ζ µ α γ α β α − πω α k C α γ α β α (cid:90) − d v (cid:107) γ (cid:107) ζ − v (cid:107) (cid:20) α + 3 γ α µ α Π α + 3 γ α µ α Π α (cid:21) e − γ α µ α Π α . (A29)It should be stressed that the compact expressions (A16), (A22) and (A29) are strictly valid for (cid:61) ζ > ζ -plane of the integrals, owing to the presence of branchpoints at ζ = ± v (cid:107) → s , such that γ (cid:107) = (cid:112) s /γ α . This gives the alternative expressions ε xx = 1 + (cid:88) α ω α k ζ µ α γ α β α − ω α k ζ µ α γ α K ( µ α ) A αxx , (A30) ε yy = 1 − (cid:88) α ω α k ζ µ α γ α K ( µ α ) A αyy , (A31) ε xy = (cid:88) α ω α k ζ µ α β α γ α − ω α k µ α γ α K ( µ α ) A αxy , (A32)where the A αlm terms are defined as A αxx = 2 µ α (cid:112) − ζ (cid:90) ∞−∞ d s χ α − s × (cid:26) γ α β α ( s + 1) / + 1 µ α (cid:20) s + 1 + 3 β α γ α ( s + 1) (cid:21) + 1 µ α (cid:20) s + 1) / + 3 β α γ α ( s + 1) / (cid:21)(cid:27) e − µ α √ s +1 , (A33) A αyy = 1 µ α γ α (cid:112) − ζ (cid:90) ∞−∞ d s s χ α − s × (cid:26) s + 1) / + γ α β α ( s + 1) / + 1 µ α (cid:20) s + 1 + 3 γ α β α ( s + 1) (cid:21) + 1 µ α (cid:20) s + 1) / + 3 γ α β α ( s + 1) / (cid:21)(cid:27) e − µ α √ s +1 , (A34) A αxy = β α µ α (cid:112) − ζ (cid:90) ∞−∞ d s χ α − s × (cid:26) s + 1) / + γ α β α ( s + 1) / + 3 µ α (cid:20) s + 1 + γ α β α ( s + 1) (cid:21) + 3 µ α (cid:20) s + 1) / + γ α β α ( s + 1) / (cid:21)(cid:27) e − µ α √ s +1 , (A35)and we have introduced χ α = γ α ζ (cid:112) − ζ . (A36) Appendix B: Low-temperature expression of the kinetic dielectric tensor
In the following we expand the dielectric tensor elements in the nonrelativistic thermal limit µ α (cid:29)
1, of particularrelevance to the background plasma particles. Our starting point is Eqs. (A30)-(A35).Let us first address ε xx , by rewriting (A33) as A αxx = β α µ α (1 − ζ ) B αxx ( µ α , χ α ) (B1)where the integral function B αxx ( µ α , χ α ) is of the form B αxx ( µ α , χ α ) = (cid:90) ∞−∞ d s f α ( s ) χ α − s e − µ α ϕ ( s ) . (B2)Applying Laplace’s method to B αxx ( µ α , χ α ) [74], one obtains to first order in 1 /µ α : B αxx ( µ α , χ α ) (cid:39) e − µ α ϕ (0) (cid:90) ∞−∞ d s e − s (cid:113) µ α χ α ϕ (cid:48)(cid:48) (0)2 − s × (cid:26) f α (0) + 1 µ α (cid:20) f (cid:48)(cid:48) α (0) ϕ (cid:48)(cid:48) (0) s − (cid:18) ϕ (4) f α (0)6 ϕ (cid:48)(cid:48) (0) + 2 ϕ (3) f (cid:48) α (0) ϕ (cid:48)(cid:48) (0) (cid:19) s + ϕ (3) (0) f α (0) ϕ (cid:48)(cid:48) (0) s (cid:21)(cid:27) . (B3)0Using ϕ (0) = 1, ϕ (cid:48)(cid:48) (0) = 1, ϕ (3) (0) = 0 and ϕ (4) (0) = −
3, as well as the expansions f α (0) (cid:39) γ α β α + (1 + 3 γ α β α ) /µ α and f (cid:48)(cid:48) α (0) (cid:39) − γ α β α , one gets B αxx ( µ α , χ α ) (cid:39) e − µ α (cid:90) ∞−∞ d s e − s ˜ χ α − s (cid:26) γ α β α + 1 µ α (cid:20) γ α β α − γ α β α s + γ α β α s (cid:21)(cid:27) , (B4)where we have defined ˜ χ α = (cid:114) µ α χ α = (cid:114) µ α γ α ζ (cid:112) − ζ . (B5)Introducing the well-known plasma dispersion function [66] Z ( η ) = 1 √ π (cid:90) ∞−∞ d s e − s x − η , (B6)and noting that Z (cid:48) ( η ) = − η Z ( η )], Eq. (B4) can be conveniently expressed as B αxx ( µ α , χ α ) (cid:39) −√ πe − µ α (cid:26)(cid:20) γ α β α + 1 µ α (1 + 3 γ α β α ) (cid:21) Z ( ˜ χ α ) + ˜ χ α γ α β α µ α (cid:20) Z (cid:48) ( ˜ χ α ) + 14 (cid:0) − ˜ χ α Z (cid:48) ( ˜ χ α ) (cid:1)(cid:21)(cid:27) . (B7)By combining this equation with Eqs. (A30) and (B1), and using the small-argument expansion K ( x ) (cid:39) (cid:112) π x e − x (1 + x ) [70], one obtains the following low-temperature approximation of ε xx : ε xx (cid:39) − (cid:88) α ω α k ζ µ α γ α (cid:26) β α Z (cid:48) ( ˜ χ α ) + 18 µ α (cid:18) − γ α (cid:19) (cid:18) Z (cid:48) ( ˜ χ α )2 (cid:19) − β α µ α ˜ χ α (cid:20) Z (cid:48) ( ˜ χ α ) + 14 (cid:0) − ˜ χ α Z (cid:48) ( ˜ χ α ) (cid:1)(cid:21)(cid:27) . (B8)Note that ˜ χ α is the correct relativistic equivalent of the standard argument ( η = (cid:112) µ α ζ ) of the Z ( η ) function involvedin the nonrelativistic CFI dispersion relation [3]. A similar result was obtained in Ref. [75] in the case of electrostaticplasma waves.Likewise, the integral involved in ε yy (A31) can be expanded to first order in 1 /µ α A αyy (cid:39) e − µ α µ α (cid:112) − ζ (cid:90) ∞−∞ d s s ˜ χ α − s e − s , (B9)so that ε yy (cid:39) − (cid:88) α ω α k (1 − ζ ) µ α γ α Z (cid:48) ( ˜ χ α ) . (B10)Finally, we expand the integral involved in ε xy (A32) as A αxy ( µ α , χ α ) (cid:39) β α e − µ α µ α (1 − ζ ) (cid:90) ∞−∞ d s e − s ˜ χ α − s (cid:26) γ α + 1 µ α (cid:20) γ α + (2 − γ α ) s + γ α s (cid:21)(cid:27) . (B11)Combining this expression with Eq. (A32) and identifying the Z and Z (cid:48) functions yields ε xy (cid:39) (cid:88) α ω α k ζ µ α β α γ α (cid:26) χ α µ α Z ( ˜ χ α ) − Z (cid:48) ( ˜ χ α ) + ˜ χ α µ α (cid:20)(cid:18) − γ α (cid:19) Z (cid:48) ( ˜ χ α ) + 12 (cid:0) − ˜ χ α Z (cid:48) ( ˜ χ α ) (cid:1)(cid:21)(cid:27) . (B12)Expressions (B8), (B10) and (B12) provide the sought-for low-temperature expansions of the CFI dielectric tensor.The involved Z and Z (cid:48) functions can be readily evaluated in the entire complex plane using, e.g. , the fast solverdeveloped in Ref. [76].1 Appendix C: Series expansion of the dielectric tensor
In similar fashion to Ref. [75], the dielectric tensor can be expanded in the form of an infinite series, which provesconvenient for deriving approximations in the kinetic regime.Let us first address ε xx , remarking that when (cid:61) χ α > A αxx = − iµ α (cid:112) − ζ (cid:90) ∞ d t e itχ α (cid:90) ∞−∞ d s e − µ α √ s +1 − its × (cid:26) γ α β α ( s + 1) / + 1 µ α (cid:20) s + 1 + 3 β α γ α ( s + 1) (cid:21) + 1 µ α (cid:20) s + 1) / + 3 β α γ α ( s + 1) / (cid:21)(cid:27) e − µ α √ s +1 , (C1)where we have exploited the identity 1 χ α − s = − i (cid:90) ∞(cid:61) χ α dt e it ( χ α − s ) , (C2)Equation (C1) can be put in the form A αxx = 4 iµ α (cid:112) − ζ (cid:90) ∞ d t e itχ α (cid:20) ∂∂b J ( t, µ α , − β α γ α ∂ ∂b J ( t, µ α , (cid:21) . (C3)where we have introduced [77] J ( t, λ, b ) = (cid:90) ∞−∞ d s e − λ √ s + b − its √ s + b = 2 K (cid:104)(cid:112) b ( λ + t ) (cid:105) . (C4)Substituting the above expression into Eq. (C1) yields A αxx = − iµ α (cid:112) − ζ (cid:90) ∞ d t e itχ α (cid:104)(cid:112) µ α + t K (cid:16)(cid:112) t + µ α (cid:17) + β α γ α (cid:0) t + µ α (cid:1) K (cid:16)(cid:112) t + µ α (cid:17)(cid:105) . (C5)Expanding e itχ α = (cid:80) ∞ n =0 ( itχ α ) n /n !, and making use of the identities [71] (cid:90) ∞ d t t n K (cid:16) x (cid:112) t + µ (cid:17) = 2 n − Γ (cid:0) n +12 (cid:1) x n +12 K n +12 ( xµ ) , (C6) d dx K ( xz ) = z K ( xz ) + K ( xz )] , (C7)we obtain after some algebra, A αxx = − √ πγ α ζ (1 − ζ ) (cid:40) i ∞ (cid:88) m =0 ( − m Γ( m + 1) (cid:18) µ α χ α (cid:19) m − (cid:20) µ α K m + ( µ α ) + β α γ α K m + ( µ α ) (cid:21) − ∞ (cid:88) m =0 ( − m Γ( m + ) (cid:18) µ α χ α (cid:19) m (cid:20) µ α K m +2 ( µ α ) + β α γ α K m +3 ( µ α ) (cid:21)(cid:41) . (C8)The above series can be further simplified using the multiplication theorem of Bessel functions [70], giving finally ε xx = 1 + (cid:88) α ω α ω µ α β α γ α + √ πω α k (1 − ζ ) µ α γ α K ( µ α ) (cid:26) i (cid:115) µ α χ α K / (cid:16) µ α (cid:112) χ α (cid:17) µ α (1 + χ α ) / + β α γ α K / (cid:16) µ α (cid:112) χ α (cid:17) (1 + χ α ) / − ∞ (cid:88) m =0 ( − m Γ( m + ) (cid:18) µ α χ α (cid:19) m (cid:20) µ α K m +2 ( µ α ) + β α γ α K m +3 ( µ α ) (cid:21)(cid:27) . (C9)Likewise, the integral involved in ε yy (A31) can be recast as A αyy = − iµ α γ α (cid:112) − ζ (cid:90) ∞ d t e itχ α (cid:20) I ( t, µ α , − µ α ∂∂b I ( t, µ α ,
1) + 4 β α γ α µ α ∂ ∂b I ( t, µ α , (cid:21) . (C10)2Here, we have defined I ( t, λ, b ) = (cid:90) ∞−∞ d s s e − λ √ s + b − its √ s + b , (C11)which can be solved in closed form as [77] I ( t, λ, b ) = − bt t + λ K (cid:104)(cid:112) b ( λ + t ) (cid:105) + 2 b / ( λ − t )( t + λ ) / K (cid:104)(cid:112) b ( λ + t ) (cid:105) . (C12)As in the previous derivation, we substitute Eqs. (C12) into Eq. (C10), expand the exponential factor and exploit [71] (cid:90) ∞ d t t λ +1 K ν (cid:2) √ t + z (cid:3) ( t + z ) ν/ = 2 λ Γ( λ + 1) z ν − λ − K ν − λ − ( z ) . (C13)There follows A αyy = − √ πγ α ζµ α (1 − ζ ) (cid:40) i ∞ (cid:88) m =0 ( − m Γ( m + 1) (cid:18) µ α χ α (cid:19) m + (cid:20) K m + ( µ α ) + 2 µ α (cid:18) β α (cid:18) m + 12 (cid:19)(cid:19) K m + ( µ α ) (cid:21) + ∞ (cid:88) m =0 ( − m Γ( m + ) (cid:18) µ α χ α (cid:19) m (cid:20) K m ( µ α ) + 2 µ α (1 + β α m ) K m +1 ( µ α ) (cid:21)(cid:41) . (C14)Applying again the multiplication theorem of Bessel functions, we find the following expression for ε yy ε yy = 1 + (cid:88) α √ πω α k (1 − ζ ) µ α γ α K ( µ α ) (cid:40) i (cid:114) µ α χ α (cid:34) K / (cid:16) µ α (cid:112) χ α (cid:17) (1 + χ α ) / + 2 µ α (cid:18) β α (cid:19) K / (cid:16) µ α (cid:112) χ α (cid:17) (1 + χ α ) / − β α χ α K / (cid:16) µ α (cid:112) χ α (cid:17) (1 + χ α ) / (cid:35) + ∞ (cid:88) m =0 ( − m Γ( m + ) (cid:18) µ α χ α (cid:19) m (cid:20) K m ( µ α ) + 2 µ α (1 + β α m ) K m +1 ( µ α ) (cid:21)(cid:41) . (C15)Reiterating the previous procedure, we first express the integral term involved in ε xy (A32) as A αxy = − iβ α µ α (cid:112) − ζ (cid:90) ∞ d t e itχ α (cid:20) µ α J ( t, µ α , − ∂∂b J ( t, µ α ,
1) +4 β α γ α ∂ ∂b J ( t, µ α , (cid:21) . (C16)Evaluating the derivatives of the function J defined by Eq. (C4), we obtain A αxy = √ πγ α β α ζ (1 − ζ ) (cid:40) i ∞ (cid:88) m =0 ( − m Γ( m + 1) (cid:18) µ α χ α (cid:19) m − (cid:20) mµ α K m + ( µ α ) − γ α K m + ( µ α ) (cid:21) − ∞ (cid:88) m =0 ( − m Γ( m + ) (cid:18) µ α χ α (cid:19) m (cid:20) (2 m + 1) µ α K m +2 ( µ α ) − γ α K m +3 ( µ α ) (cid:21)(cid:27) . (C17)Finally, Bessel-function identities allow us to obtain the following alternative expression for ε xy ε xy = (cid:88) α ω α k ζ µ α γ α β α + √ πω α ζk (1 − ζ ) µ α γ α β α K ( µ α ) (cid:40) i (cid:114) µ α χ α (cid:18) γ α χ α (cid:19) K / (cid:16) µ α (cid:112) χ α (cid:17) (1 + χ α ) / + ∞ (cid:88) m =0 ( − m Γ( m + ) (cid:18) µ α χ α (cid:19) m (cid:2) (2 m + 1) K m +2 ( µ α ) − µ α γ α K m +3 ( µ α ) (cid:3)(cid:41) . (C18) [1] E. S. Weibel, Phys. Rev. Lett. , 83 (1959). [2] B. D. Fried, Phys. Fluids , 337 (1959). [3] R. C. Davidson, D. A. Hammer, I. Haber, and C. E.Wagner, Phys. Fluids , 317 (1972).[4] F. Califano, F. Pegoraro, and S. V. Bulanov, Phys. Rev.E , 963 (1997).[5] A. Achterberg and J. Wiersma, Astron. Astrophys. ,1 (2007).[6] S. S. Moiseev and R. Z. Sagdeev, Journal of Nuclear En-ergy. Part C, Plasma Physics, Accelerators, Thermonu-clear Research , 43 (1963).[7] M. V. Medvedev and A. Loeb, Astrophys. J. , 697(1999).[8] A. Gruzinov and E. Waxman, Astrophys. J. , 852(1999), astro-ph/9807111.[9] Y. Sentoku, K. Mima, P. Kaw, and K. Nishikawa, Phys.Rev. Lett. , 155001 (2003).[10] J. C. Adam, A. H´eron, and G. Laval, Phys. Rev. Lett. , 205006 (2006).[11] B. Allen, V. Yakimenko, M. Babzien, M. Fedurin,K. Kusche, and P. Muggli, Phys. Rev. Lett. , 185007(2012).[12] A. Debayle, J. J. Honrubia, E. d’Humi`eres, and V. T.Tikhonchuk, Phys. Rev. E , 036405 (2010).[13] P. E. Masson-Laborde, W. Rozmus, Z. Peng, D. Pesme,S. H¨uller, M. Casanova, V. Y. Bychenkov, T. Chapman,and P. Loiseau, Phys. Plasmas , 092704 (2010).[14] S. Mondal, V. Narayanan, W. J. Ding, A. D. Lad, B. Hao,S. Ahmad, W. M. Wang, Z. M. Sheng, S. Sengupta,P. Kaw, A. Das, and G. R. Kumar, Proc. Natl. Acad.Sci. USA , 8011 (2012).[15] K. Quinn, L. Romagnani, B. Ramakrishna, G. Sarri,M. E. Dieckmann, P. A. Wilson, J. Fuchs, L. Lancia,A. Pipahl, T. Toncian, O. Willi, R. J. Clarke, M. Not-ley, A. Macchi, and M. Borghesi, Phys. Rev. Lett. ,135001 (2012).[16] F. Fiuza, R. A. Fonseca, J. Tonge, W. B. Mori, and L. O.Silva, Phys. Rev. Lett. , 235004 (2012).[17] C. Ruyer, L. Gremillet, and G. Bonnaud, Phys. Plasmas , 082107 (2015).[18] W. Fox, G. Fiksel, A. Bhattacharjee, P.-Y. Chang,K. Germaschewski, S. X. Hu, and P. M. Nilson, Phys.Rev. Lett. , 225002 (2013).[19] C. M. Huntington, F. Fiuza, J. S. Ross, A. B. Zyl-stra, R. P. Drake, D. H. Froula, G. Gregori, N. L. Kug-land, C. C. Kuranz, M. C. Levy, C. K. Li, J. Meinecke,T. Morita, R. Petrasso, C. Plechaty, B. A. Remington,D. D. Ryutov, Y. Sakawa, A. Spitkovsky, H. Takabe, andH.-S. Park, Nat. Phys. , 173 (2015).[20] R. P. Drake and G. Gregori, Astrophys. J. , 171(2012).[21] H. Chen, F. Fiuza, A. Link, A. Hazi, M. Hill, D. Hoarty,S. James, S. Kerr, D. D. Meyerhofer, J. Myatt, J. Park,Y. Sentoku, and G. J. Williams, Phys. Rev. Lett. ,215001 (2015).[22] M. Lobet, C. Ruyer, A. Debayle, E. d’Humi`eres,M. Grech, M. Lemoine, and L. Gremillet, Phys. Rev.Lett. , 215003 (2015).[23] C. Ruyer, L. Gremillet, G. Bonnaud, and C. Riconda,Phys. Rev. Lett. , 065001 (2016); Phys. Plasmas ,041409 (2017).[24] J. Wiersma and A. Achterberg, Astron. Astrophys. ,365 (2004).[25] Y. Lyubarsky and D. Eichler, Astrophy. J. , 1250(2006).[26] M. V. Medvedev, M. Fiore, R. A. Fonseca, L. O. Silva, and W. B. Mori, Astrophys. J. , L75 (2004).[27] M. Milosavljevic and E. Nakar, Astrophys. J. , 978(2006).[28] A. Achterberg, J. Wiersma, and C. A. Norman, Astron.Astrophys. , 19 (2007).[29] A. Bret, A. Stockem, D. Fiuza, C. Ruyer, L. Gremillet,R. Narayan, and L. O. Silva, Phys. Plasmas , 042102(2013).[30] C. Ruyer, L. Gremillet, A. Debayle, and G. Bonnaud,Phys. Plasmas , 032102 (2015).[31] A. Vanthieghem, M. Lemoine, and L. Gremillet, Phys.Plasmas , 072115 (2018).[32] R. Lee and M. Lampe, Phys. Rev. Lett. , 1390 (1973).[33] L. O. Silva, R. A. Fonseca, J. W. Tonge, J. M. Dawson,W. B. Mori, and M. V. Medvedev, Astrophys. J. Lett. , L121 (2003).[34] J. T. Frederiksen, C. B. Hededal, T. Haugbølle, and˚A. Nordlund, Astrophys. J. Lett. , L13 (2004).[35] C. H. Jaroschek, H. Lesch, and R. A. Treumann, Astro-phys. J. , 822 (2005).[36] K.-I. Nishikawa, J. Niemiec, P. E. Hardee, M. Medvedev,H. Sol, Y. Mizuno, B. Zhang, M. Pohl, M. Oka, andD. H. Hartmann, Astrophys. J. Lett. , L10 (2009).[37] G. Shvets, O. Polomarov, V. Khudik, C. Siemon, andI. Kaganovich, Phys. Plasmas , 056303 (2009).[38] A. Bret, L. Gremillet, and M. E. Dieckman, Phys. Plas-mas , 120501 (2010).[39] T. N. Kato, Astrophys. J. Lett. , 974 (2007).[40] A. Spitkovsky, Astrophys. J. Lett. , L5 (2008).[41] S. F. Martins, R. A. Fonseca, L. O. Silva, and W. B.Mori, Astrophys. J. Lett. , L189 (2009).[42] U. Keshet, B. Katz, A. Spitkovsky, and E. Waxman,Astrophys. J. , L127 (2009).[43] L. Sironi and A. Spitkovsky, Astrophys. J. , 1523(2009).[44] L. Sironi, A. Spitkovsky, and J. Arons, Astrophys. J. , 54 (2013).[45] T. Haugbølle, Astrophys. J. Lett. , 42 (2011).[46] R. Kumar, D. Eichler, and M. Gedalin, Astrophys. J. , 165 (2015).[47] M. Lemoine and G. Pelletier, Mon. Not. Roy. Astron.Soc. , 321 (2010).[48] I. Rabinak, B. Katz, and E. Waxman, Astrophys. J. ,157 (2011).[49] M. Lemoine and G. Pelletier, Mon. Not. Roy. Astron.Soc. , 1148 (2011).[50] R. Shaisultanov, Y. Lyubarsky, and D. Eichler, Astro-phys. J. , 182 (2012).[51] M. Lemoine, L. Gremillet, G. Pelletier, and A. Van-thieghem, Phys. Rev. Lett. , 035101 (2019),arXiv:1907.07595.[52] M. Lemoine, A. Vanthieghem, G. Pelletier, andL. Gremillet, Phys. Rev. E, submitted (Pap. II) (2019),arXiv:1907.08219 [astro-ph.HE].[53] M. Lemoine, G. Pelletier, A. Vanthieghem, andL. Gremillet, Phys. Rev. E, submitted (Pap. III) (2019),arXiv:1907.10294 [astro-ph.HE].[54] A. Vanthieghem, M. Lemoine, L. Gremillet, and G. Pel-letier, Phys. Rev. E, in prep. (Pap. IV) (2019).[55] R. D. Blandford and C. F. McKee, Phys. Fluids , 1130(1976).[56] L. O. Silva, R. A. Fonseca, J. W. Tonge, W. B. Mori,and J. M. Dawson, Phys. Plasmas , 2458 (2002).[57] S. Ichimaru, Basic Principles of Plasma Physics: A Sta- tistical Approach , A lecture note and reprint series (W.A. Benjamin, Reading, MA (USA), 1973).[58] K. Molvig, Phys. Rev. Lett. , 1504 (1975).[59] J. R. Cary, L. E. Thode, D. S. Lemons, M. E. Jones, andM. A. Mostrom, Phys. Fluids , 1818 (1981).[60] T. Okada and K. Niu, J. Plasma Phys. , 423 (1980).[61] J. M. Hill, M. H. Key, S. P. Hatchett, and R. R. Freeman,Phys. Plasmas , 082304 (2005).[62] A. Bret, M.-C. Firpo, and C. Deutsch, Phys. Rev. E ,046401 (2004).[63] A. Bret, L. Gremillet, and J. C. Bellido, Phys. Plasmas , 032103 (2007).[64] F. J¨uttner, Annalen der Physik , 856 (1911).[65] T. P. Wright and G. R. Hadley, Phys. Rev. A , 686(1975).[66] J. D. Huba, Plasma Physics (Naval Research Laboratory,Washington, DC, 2013) pp. 1–71.[67] E. Lefebvre, N. Cochet, S. Fritzler, V. Malka, M.-M. Al´eonard, J.-F. Chemin, S. Darbon, L. Disdier,J. Faure, A. Fedotoff, O. Landoas, G. Malka, V. M´eot,P. Morel, M. Rabec LeGloahec, A. Rouyer, C. Rubbe-lynck, V. Tikhonchuk, R. Wrobel, P. Audebert, andC. Rousseaux, Nucl. Fus. , 629 (2003). [68] B. B. Godfrey and J.-L. Vay, J. Comput. Phys. , 1(2014).[69] A. Bret, L. Gremillet, and D. B´enisti, Phys. Rev. E ,036402 (2010).[70] M. Abramowitz and I. A. Stegun, Handbook of Mathe-matical Functions (Dover Publications, New York, 1972).[71] I. S. Gradshteyn and I. M. Rizhik,