Physics of collisionless shocks - theory and simulation

Collisionless shocks occur in various fields of physics. In the context of space and astrophysics they have been investigated for many decades. However, a thorough understanding of shock formation and particle acceleration is still missing. Collisionless shocks can be distinguished into electromagnetic and electrostatic shocks. Electromagnetic shocks are of importance mainly in astrophysical environments and they are mediated by the Weibel or filamentation instability. In such shocks, charged particles gain energy by diffusive shock acceleration. Electrostatic shocks are characterized by a strong electrostatic field, which leads to electron trapping. Ions are accelerated by reflection from the electrostatic potential. Shock formation and particle acceleration will be discussed in theory and simulations.


Introduction
Collisionless shocks show interesting features for particle acceleration, which is why they are of importance in many fields of physics [1]- [7]. They are generated by the interaction of charged particles with the surrounding self-generated fields -in contrast to collisiondominated shocks, where two-particle interactions determine the physical behaviour. In a simple setup of two interpenetrating plasma slabs, either a strong electrostatic field [8] is generated due to the two-stream instability [9], or an electromagnetic field due to Weibel-like instabilities [10,11]. The feedback of such fields mediates the shock formation process.
The Weibel or filamentation instability generate a turbulent field in which the charged particles are deflected. Parallel momentum is transferred into perpendicular momentum, and thus the particles are accumulated and the plasma slabs are compressed. In relativistic initially unmagnetised plasmas, the density compression is three times the initial density [12]. The time scales of the shock formation are determined by the time scales of the electromagnetic instability [13]. The details are given in Sec. 3.
Electrostatic (ES) shocks form in plasmas with different components of electrons and (heavy) ions. A high mass and temperature difference is favourable for the generation of these shocks. An electrostatic potential is built-up which traps electrons in the downstream region and the electron distribution function is widened in parallel direction. The perpendicular components are less affected. This temperature anisotropy gives rise to the electromagnetic Weibel instability. The time scales of the formation of such electromagnetic modes are calculated in Sec. 4 and compared against the time scales of shock formation in order to determine their relevance during the shock formation process. We perform the analysis for plasmas with relativistic temperatures and fluid velocities. Laser-generated plasmas in the laboratory are now entering the relativistic regime, which is why it is important to explore the physics also in this parameter range.

Initial setup for shock formation
We study shock formation in a simple symmetric setup (see Fig. 1). Let us consider two charge-neutral counterstreaming beams of electrons and positrons or ions with bulk velocities ±v 0 , temperatures T e , T i and thermal parameter µ = mc 2 /k B T . The overlapping region turns unstable due to collisionless plasma instabilities (Phase 1), which will mediate two shocks propagating into the upstream regions (Phase 2). The collisionless shock is of electromagnetic or electrostatic nature, depending on the initial parameters for the velocities and temperatures. We will discuss these details in the following sections.

Shock formation in electromagnetic shocks
Electromagnetic shocks develop if electromagnetic instabilities dominate in the overlapping region. In such a symmetric setup, the electromagnetic filamentation (Weibel) instability is the fastest mode for cold beams T e ≈ T i ≈ 0 and relativistic velocities v 0 /c ≈ 1. The growth rate of the cold filamentation instability is given by δ a = 2/γ 0 β 0 ω pa with β 0 = v 0 /c, Lorentz factor γ 0 = (1−β 2 0 ) −1/2 and plasma frequency of species a being ω pa = 4πn 0 q 2 a /m a . The saturation time of this instability is a function of the initial and final magnetic field strengths. The initial magnetic field strength B i is obtained from the evaluation of the spectra of spontaneous magnetic fluctuations, while the final field B f ≃ 8πγ 0 n 0 mc 2 results from a trapping condition in the magnetic field structure [14]. The saturation time in electron-positron shocks can then be expressed as For electron-positron pair shocks, the shock formation time is simply twice the saturation time of the instability in two dimensions, τ f,e = 2τ s,e , and a factor 3 in three dimensions, τ f,e = 3τ s,e [15]. For electron-ion shocks the theory has to be extended by an extra term for the merging time of the filaments. At the saturation time of the filamentation modes in electron-positron plasmas, the transverse size of the filaments is already large enough in order to deflect particles strongly enough for efficient accumulation of particles. On the other hand, in electron-ion plasmas, the filaments are still on the electron scale. An additional merging time is required in order to bring the filaments to the required size to significantly deflect the ions [13]. The total shock formation time in electron-ion shocks is thus given by with τ s,i the saturation time of the cold ion filamentation instability and d the number of dimensions. The predicted scaling is broadly consistent with particle-in-cell simulation results.

Shock formation in electrostatic shocks
Electrostatic shocks require a mass and temperature difference between ions and electrons. The electrostatic two-stream instability is dominant in the overlapping region, which is fast for low streaming velocities. A strong electrostatic potential is generated in the downstream region which can trap electrons (see Fig. 2), which can be approximated as a flat-top distribution. In the non-relativistic case, we follow Ref. [16] for a steadystate solution of the shock and describe the electron distribution by a population of free streaming and a population of trapped electrons. If we introduce normalisations for the velocity, β = v/c, and the electrostatic potential, ϕ = eφ/m e c 2 , electrons are free streaming if the condition |β x | > √ ϕ is fulfilled and trapped otherwise. The nonrelativistic distribution of the electrons is then given by [16] f e = C 0 µ 2π

Relativistic generalisation of the distribution function
For the relativistic generalisation of the electron distribution eq. (5), we start with a Maxwell-Jüttner distribution in the mean rest frame (subscript R) with the Lorentz factor always defined as γ = (1 − β 2 ) −1/2 , and perform a Lorentz transformation with γ 0 into the moving (= laboratory) frame (subscript L) in ±x directions. Thus, with γ R = γ 0 γ L (1 ± β 0 β L,x ) we obtain the relativistic distribution in the laboratory frame with µ L = µ R γ 0 . C R and C L are normalisation constants in the mean rest frame and in the laboratory frame, respectively. We now introduce the dependency on the electrostatic potential ϕ. From energy conservation, we obtain a balance between the kinetic energy in the upstream, where the bulk is moving with γ L , and the total energy in the shock transition region, where kinetic energy is transformed into potential energy. The bulk Lorentz factor is reduced to γ ′ L . This is expressed by Electrostatic shocks are mainly one-dimensional, which is why we assume that the longitudinal and transverse processes can be separated. Introducing We can then write the relativistic generalisation of the electron distribution in eq. (5), if we replace γ L and β L,x in eq. (7) by their ϕ dependent terms and with γ := γ ′ L and with the definitions γ = γ(u) = √ 1 + u 2 and γ ⊥ = 1 + u 2 ⊥ . The definition of the trapping velocity u c is derived from the balance of electrostatic and kinetic energy in eq. (8) for γ L = 1 which gives γ = γ c := 1 + ϕ and u c = γ 2 c − 1. The non-relativistic approximation of this trapping velocity agrees with the non-relativistic definition in eq.
The normalisation constant C r0 is obtained from Eqs. (5) and (9) describe the electron distribution in the quasi-steady ES shock. In fig. 3 we show the change of the distribution function from the initial Maxwell distribution for ϕ = 0 to a flat-top distribution for ϕ > 0 for u 0 = β 0 γ 0 = 0.01 and µ = 50, corresponding to k B T e = 10 keV. The relativistic expression eq. (9) is also compared against the non-relativistic expression eq. (5). In the far upstream, the distribution is Maxwellian, see fig. 3a. In the shock transition region, where 0 < ϕ ≤ ϕ max , the distribution f (u x ) broadens and becomes flat-top. ϕ max denotes the saturation value of the electrostatic potential in the far downstream which can be obtained by solving numerically the differential equation 1 2 ∂ϕ ∂x 2 + Ψ(ϕ) = 0 with the Sagdeev potential Ψ(ϕ) [17,18]. Fig. 3b shows the distribution function at the ion reflection condition, where the kinetic energy of the ions equals the electrostatic energy. In normalised quantities, this is given by ϕ = ϕ ref l := (γ 0 − 1)m i /m e . In the far downstream, for ϕ = ϕ max the non-relativistic description breaks since ϕ max > 1. In this configuration, the electrons are usually treated kinetically, while the higher inertia ions are described with a fluid model. In the following section, we summarise the initial conditions for which an ES shock forms.

Conditions for ES shock formation
Here, we summarise the results on ES shock formation conditions. In Refs. [19,20] a condition for the maximum Mach number was found as 1 < M max 3.1, with the Mach number M defined as the ratio of the upstream velocity in the shock rest frame to the ion sound speed, M = v ′ 0 /c s = (v 0 + v sh )/c s . The shock velocity in the laboratory frame can be approximated from the shock jump conditions as v sh /c = √ γ ad − 1 (γ 0 − 1)/(γ 0 + 1) with the ideal gas adiabatic constant γ ad [12]. The ion sound speed is given by c s = k B T e /m i . This can be generalised for relativistic plasmas with the relativistic Mach number M = u ′ 0 /u s , where u ′ 0 = β ′ 0 γ ′ 0 is the dimensionless momentum in the shock rest frame and u s = β s γ s with β s = k B T e /m i c 2 [19,21]. This imposes a condition for the upstream fluid velocity and electron temperature

Calculation of unstable modes
During the early stage of ES shock formation the broadening of the electron distribution, which we observed in the previous section, occurs mainly in the longitudinal direction. The transverse directions stay almost unaffected. This has also been observed by particle-in-cell simulations [22]. The generated temperature anisotropy in the electron distribution gives rise to electromagnetic modes.
We develop now a model to describe the growth rate of the electromagnetic modes in such a setup. Later, we will compare them to the shock formation time scales and the growth rate of the cold ion-ion instability in order to determine their relevance.

Dispersion relation of EM waves in relativistic plasmas
We start from the electron distribution in eq. (9) in order to evaluate the dispersion relation of electromagnetic waves in plasma and We consider only fluctuations k = ke z perpendicular to the fluid velocity u 0 = u 0 e x in order to simplify the geometry. An evaluation of the integrals leads to and where the last integrations have to be done numerically for the general case. We use eqs. (16) and (17) to solve the dispersion relation (13) numerically to obtain the growth rate σ(k) := ℑ(ω(k)). In Figs. 4a and b we plot the maximum growth rate σ max := max(σ(k)) for different values of µ and u 0 . The maximum growth rate was evaluated at the ion reflection condition ϕ ref l = (γ 0 −1)m i /m e . The role of the potential will be discussed in sec. 5.3.
σ max shows a maximum as a function of the velocity. The higher the electron temperature, i.e. the lower µ ∝ T −1 e , the higher the velocity at which the maximum of σ max appears ( fig. 4a). On the other hand, σ max increases with µ. For an easier interpretation, we plotted the dependence of σ max against the electron temperature in Fig. 4b. The higher the initial electron temperature, the lower is the growth rate σ max .

Dispersion relation of EM waves in non-relativistic plasmas
In order to further discuss the properties of this model and to check the consistency with previous models, dispersion relation (13) is approximated for non-relativistic velocities, i.e. β 0 ≪ 1 and µ ≫ 1. We obtain U n ≈ −1 and For zero beam velocity and ϕ = 0, this is exactly the well-known solution of the Weibel instability [11]. For small frequencies ω ≪ kc, the plasma dispersion function Z can be further approximated and we get The dispersion relation then reads with V (ϕ) = C 0 e µϕ erfc √ µϕ + 2 µϕ π + 4 3 µ 3 ϕ 3 π e −µβ 2 0 /2 . The solution of the dispersion relation can now be derived analytically, which is given by with k 2 0 c 2 = ω 2 pe (V (ϕ) − 1)/3 the location of the maximum and the maximum growth rate For comparison, fig. 5 shows the relativistically correct solution of the dispersion relation, obtained from eqs. (16) and (17), and the non-relativistic approximations from eqs. (20) and (21). The parameters used are β 0 = 0.01, µ = 100 and ϕ = 0.05 (γ 0 −1) m i /m e . The maximum growth rate matches with the approximation in eq. (22), σ max = 2.0×10 −3 ω pe and the location of the maximum at k 0 = 0.25 ω pe /c.

The role of the electrostatic potential
The electrostatic potential ϕ is quickly built up during the electrostatic shock formation. Nevertheless, the electromagnetic growth rate σ(k) depends on ϕ and thus, the dependence should be investigated. We calculate the growth rate numerically for values 0 ≤ ϕ ≤ ϕ ref l , which are plotted in fig. 6. We find a power-law dependence of the growth rate on the electrostatic potential, which gives with α > 0. This ϕ dependence can be interpreted as the changing growth rate across the steady-state shock. If we assume, as a simple approach, that the potential grows linear in time, we can connect the growth rate σ max with time, which will be further discussed in the next section.

Comparison of different time scales
In order to determine the relevance of the electromagnetic electron Weibel modes during the ES shock formation process, we compare the growth rate given by eq. (13), which we label from now on σ EM,ee , with other processes.

ES Shock formation time
The time scales of ES shock formation are determined by the growth rate of the ion-ion electrostatic instability which is given by [9] σ ES,ii = 1 Estimating the shock formation time with 5 times the inverse maximum growth rate, we obtain t sf = 10γ −3/2 0 m i /m e ω −1 pe , which was confirmed in simulations (see [23]).

Growth rate of the electromagnetic cold ion-ion instability
Another competing process in ES shocks is the cold ion-ion filamentation instability, which has a growth rate [24] σ EM,ii = β 0 2 γ 0 ω pi .
The EM mode of the ion-ion instability (25) grows faster than the ES mode (24) for fluid velocities v 0 /c > 1/3. In the parameter range v 0 /c ≤ 0.1, which we looked on in this paper, the ion-ion EM modes can be neglected.

Dominant regimes
A comparison of the ES shock formation time scales with the EM modes makes it possible to determine parameter regimes, for which a shock stays electrostatic. For this, we compare the growth rate from eq. (13) with eq. (24) for different electron temperatures, expressed by k B T e /m e c 2 , and fluid velocities, given by u 0 = β 0 γ 0 . For simplicity, we choose an electrostatic potential ϕ = ϕ ref l . In the dark green region in fig. 7 the growth rate of the EM electron instability σ EM,ee is smaller than the ES ion-ion growth rate σ ES,ii . We label this domain the purely ES domain. In contrast, in the light green region, an ES shock will turn EM since σ ES,ii < σ EM,ee . In the white region (EM), no electrostatic shock develops because condition (12) is not fulfilled. This regime is electromagnetically dominated from the beginning, with shock formation according to section 3.
A more detailed analysis of the scenario in fig. 7 has been done in fig. 8 for different values of the electrostatic potential ϕ. For different fluid velocities u 0 and temperature parameters µ = m e c 2 /k B T e , the maximum growth rate σ max is calculated as in fig. 6. The characteristic time is then calculated as t char ≈ 5/σ max and plotted against the time of shock formation, where we assumed that the potential grows linearly as ϕ ∝ t.  The grey area in fig. 8 shows the region, where EM modes develop faster than the ES shock, which is the case e.g. for u 0 = 0.01 for µ 10 ( 51 keV) or for u 0 = 0.1 for µ 0.1 ( 5 MeV). Fig. 8 also shows that a detailed analysis of the growth rate is not necessary, since the characteristic time t char quickly saturates. The rough estimate in fig. 6 is thus sufficient for a qualitative determination of the regimes.

Summary and conclusions
Electromagnetic and electrostatic shocks can both develop from the same simple symmetric setup of counterstreaming beams. The choice of the initial plasma parameters determines the final shock character.
Electromagnetic shocks develop in plasmas with large beam velocities and low temperatures where the electromagnetic filamentation instability is the fastest mode. In electron-positron pair plasmas, the shock formation time is simply twice the saturation time of the instability in 2D. In electron-ion plasmas an extra merging time is necessary in order to obtain the right scale of the filaments to efficiently scatter the ions. Thus, the shock formation in electron-ion plasmas is delayed by almost a factor 3 compared with the pair shock, τ f,i ω pi ≈ 3τ f,e ω pe . Electrostatic shocks form in electron-ion plasmas with small beam velocities and large electron temperatures. Here, the electrostatic two-stream instability dominates. An electrostatic potential is generated in which electrons are trapped. The subsequent deformation of the electron distribution gives rise to electromagnetic Weibel modes which can destroy the electrostatic shock features.