Interaction between electrostatic collisionless shocks generates strong magnetic fields

The head-on collision between electrostatic shocks is studied via multi-dimensional particle-in-cell simulations. A strong magnetic field develops after the interaction, which causes the shock velocities to drop significantly. This transverse magnetic field is generated by the Weibel instability, which is driven by pressure anisotropies due to longitudinal electron heating while the shocks approach each other. The possibility to explore the physics underpinning the shock collision in the laboratory with current laser facilities is discussed.


Introduction
Collisionless shocks, i.e. shocks mediated solely by electromagnetic fields, are pervasive in space, astrophysics and laboratory plasmas (see discussion in [1][2][3]). In a collisionless unmagnetised electron-ion plasma with electrons much hotter than ions, shocks may be sustained only by an ambipolar electric field, which develops because of the different inertia of the particles [3,4]. These shocks are electrostatic in nature and sometimes are referred to as ion acoustic shocks (exceptions include the so called Weibel mediated shocks [5][6][7], electromagnetic shocks, where the shock formation is connected to the magnetic field turbulence generated by the Weibel instability [8]).
Electrostatic shocks derive from processes of wave steepening and breaking of ion-acoustic modes and present no magnetic field jump [3,6]. They have been observed in space [9,10] and laboratory plasmas [11][12][13] and are thought to play a role during the late expansion phase of supernova remnant blast shells [14]. It was shown that their formation and evolution exhibit remarkable similarities with perpendicular magnetic shocks more common in astrophysics [15]. The theoretical framework to describe these shocks was first developed by Sagdeev [16], who demonstrated that the electrostatic potential of an ion acoustic soliton and its properties could be computed using the equation of a particle moving in a pseudopotential, which became later known as Sagdeev potential. A shock solution is then obtained by considering mechanisms that break the soliton symmetry [17], such as Landau damping [17], upstream reflection of ions [4,18,19] or electron trapping in the downstream [20][21][22]. Pioneering numerical simulations based on the particle-in-cell (PIC) technique granted insight into the shock formation mechanism and particle dynamics [15,23]. In these simulations, electrostatic shocks were generated through the collision of plasma flows with sub-relativistic bulk speeds and large electron thermal velocities [6]. Ion acoustic shocks have also been triggered in simulations by letting a dense plasma expanding into a rarefied medium [24,25]. While these works focus mainly on the initial stages of the shock development, more recently the long time evolution of ion acoustic shocks has been explored. Kato and Takabe showed that the presence of counterstreaming ions upstream of the shock may result into the oblique electrostatic ion-ion instability [26,27] and the transverse ion current filamentation instability [28,29], which eventually lead to the shock being destroyed [30]. Stockem et al found that the electron Weibel instability [8,[31][32][33][34] can develop downstream of the shock due to electrons being heated longitudinally and thus becoming anisotropic [35]. Sarri et al did not detect any Weibel-generated magnetic field in their simulations of electrostatic shocks [36]. However, they observed filamentary structures across the shock, which were attributed to the ion-ion instability [36]. Differently from Stockem et al, where the shock forms from counterstreaming plasma clouds [35], Sarri et al drove ion acoustic shocks by letting denser plasmas expanding into more rarefied ones [36]. In this case, electron mixing at the discontinuity interface acts to keep the anisotropy under the threshold for the development of the Weibel instability [36].
In recent years, collisionless electrostatic shocks have attracted considerable attention in the quest for a compact particle accelerator based on laser and plasma technologies [37][38][39][40]. Ions at rest are reflected by the electrostatic potential of the shock to twice the shock velocity, making the electrostatic shock an efficient mechanism for particle acceleration. The generation of electrostatic shocks in laser-driven plasmas and the production of high-quality proton beams have now been demonstrated in several experiments [41][42][43][44][45][46][47][48][49]. Computer simulations performed under realistic experimental conditions revealed that a shock is launched when the radiation pressure exerted by the laser steepens the plasma density and creates a density spike and, at the same time, target electrons are heated to high temperatures, such that T e T i , with T e and T i electron and ion temperatures, respectively [50][51][52]. To ensure that the shock propagates with a constant speed, the electron temperature should be as uniform as possible over the plasma [50,51]. This places a constraint on the target size and inhibits the development of temperature anisotropy instabilities prior to shock formation.
In this work, we will explore for the first time (at least to the authors' knowledge) the microphysics underpinning the head-on collision between two electrostatic shocks. The interaction between collisionless shocks, independently on their electrostatic or electromagnetic nature, is a fundamental physical process (see for instance [53]) and understanding how these collisionless non-linear waves interact represents a major advance in plasma physics. Performing ab initio PIC simulations, we investigate the microphysics of the interaction in detail and identify a regime for shock collisions, where kinetic effects determine the interaction dynamics. The results indicate that the collision does not change the nature of the waves; the shocks can pass through each other and proceed in opposite directions, similarly to what is predicted for electrostatic solitons [54]. However, the collision is highly inelastic causing the shocks to slow down considerably after the interaction. A strong perpendicular magnetic field develops right after the collision. The field is generated by the Weibel instability [8,[31][32][33][34], which is driven by an electron pressure anisotropy in the interaction region, caused by strong longitudinal heating while the shocks approach. Given that electrostatic shocks have now been produced in several experiments [41][42][43][44][45][46][47][48][49], we suggest that the binary collision of shock waves could be explored in the laboratory exploiting near critical density target(s) heated by laser beams. Thus, we provide estimates of the laser parameters and target density which will allow for the generation of counterpropagating shocks such that the magnetic field can be observed after their collision. The proposed setup could then be used not only to study how collisionless shocks interact, but also to explore magnetic field generation at electron scales via Weibel instability in the laboratory [55][56][57][58][59][60]. Experiments based on this configuration could provide insight into astrophysical environments, such as gamma-ray bursts, supernovae remnants, active galactic nuclei, where the instability is believed to produce very intense fields [61][62][63][64].

The microphysics of the interaction between collisionless shocks
To study the interaction of two collisionless electrostatic shocks, we performed one-and two-dimensional (1D and 2D, respectively) simulations using the PIC code OSIRIS [65]. A plasma composed by uniformly hot Maxwellian electrons and cold hydrogen ions with m i /m e = 1836, where m i and m e are the proton and electron mass, respectively, akin to a laser-generated system, is considered. The two species do not have any initial drift velocity. This is a valid assumption for laser pulses interacting with near critical density targets, where most of the laser energy is transferred to the plasma in the form of electron heating [50]. The plasma density profile presents two perfectly symmetric sharp discontinuities (inset in figure 1(a)). They constitute the jump in density that triggers the shock formation, similarly to that explained in [25,51]. We carried out simulations considering electron temperatures in the range T e = 0.02 − 1.5 MeV and initial density jumps in the interval σ ≡ n 0,M /n 0,m = 3.33-10, where M and m indicate the highest and the lowest density respectively. Similar electron temperatures and density spikes (with respect to the background plasma density) have been measured in PIC simulations of laser-driven electrostatic shocks [42,45]. The largest simulation window was 4800c/ω p long and 240c/ω p wide, with c/ω p the electron skin depth, c the speed of light in vacuum, ω p ≡ 4πe 2 n 0,M /m e the electron plasma frequency and e the elementary charge. The spatial discretisation was selected to resolve the smallest electron scale; thus the smaller between the skin depth or the electron Debye length λ D ≡ k B T e /4πe 2 n 0,M , where k B is the Boltzmann constant, is resolved with four cells. To model the plasma dynamics correctly 36 (or 1000 in 1D) particles per cell per species and a fourth-order interpolation scheme were employed. Periodic boundary conditions have been used, but the simulation window is large enough that they do not interfere with the plasma dynamics. Our PIC simulations do not include Coulomb collisions. For typical laser-driven plasmas of relevance for electrostatic shock generation (electron temperatures T e ∼ MeV, ion temperatures T i ∼ hundreds eV and densities 10 19 -10 21 cm −3 [50]), the mean free path is much larger than the system. Therefore collisional effects can be neglected.
In figure 1, the main stages of a typical simulation are reported. In this case, T e = 1.5 MeV and σ = 10 were used. When the more dense plasma expands into the less dense component, it drives several nonlinear waves that eventually develop into two shock waves streaming towards each other in opposite directions. The shock structures are clearly recognizable at x 1 ±35c/ω p in figures 1(a) and (d), which depict the longitudinal ion phase space and the longitudinal electric field E 1 at t = 335ω −1 p . In particular, E 1 presents a double layer structure, which is a typical signature of electrostatic shocks [4,16,21]. The shock waves collide in the middle of the simulation window (x 1 = 0) at t = 913.5ω −1 p (figures 1(b) and (e)). The collision results in an enhancement of the electric field (figure 1(e)). Remarkably, until t = 913.5ω −1 p , the magnetic field in the direction out of the box B 3 is negligible (figures 1(g) and (h)). After the collision, the shocks pass through each other and continue streaming in opposite directions (figures 1(c) and (f)). At that stage, B 3 in the middle of the simulation window presents an intense filamentary structure, which has developed only after the interaction (figure 1(i)). The field is confined in a small region between x 1 = [−30, 30]c/ω p downstream of the shocks. The filaments are characterised by a wavelength λ 50.7c/ω p . The field reaches a maximum at around t 1309ω −1 p and then starts to decay. The measured shock velocities before and after the interaction are shown in figure 2(a), along with the evolution of the ion density averaged along the x 2 direction. The shocks move at a constant speed v s = ±0.056c corresponding to a Mach number M = v s /c s = 1.4, where c s = k B T e /m i is the sound speed for ions of mass m i , in excellent agreement with the theoretical model in [4]. After the collision, the speed of both shocks drops to v s ±0.032c, which corresponds to a reduction Δv s /v s of about 42%. Despite the different simulation setup, the slowdown agrees moderately well with the estimates in [66], which predicts a decrease in velocity ≈1/ √ 2 and ≈1/ √ 3 in 2D and 3D simulations of interpenetrating plasmas following the development of a Weibel-like instability. The reasons behind this slowdown will be clear after the analysis of the particle distribution functions, presented in the next paragraphs. For comparison, we report the same results as obtained from a 1D simulation (figure 2(b)), where, due to the geometry, the growth of modes perpendicular to x 1 is inhibited and therefore no B 3 is observed. In this case, the shocks only slightly slow down: their velocity passes from v s = ±0.058c to v s ±0.054c after the collision. We verified that despite the slowdown, the nature of the shock waves remain unchanged: the electric field maintains the double layer structure (figures 1(d) and (f)), the density jump conditions for low Mach number shocks [15] are satisfied (figures 2(c) and (d)) and the waves move with M > 1. These results indicate that the multidimensional interaction of collisionless shocks is an inelastic process, and that the waves lose energy during the collision much more substantially in 2D than 1D. Given this and the prediction in [66], we expect that in 3D simulations the shock velocity after the collision would be even lower. Indeed, a preliminary 3D simulation performed using the same parameters seems to confirm this trend.
To gain insight on what happens immediately after the collision, we have checked the electron distribution functions f(p 1 ) and f (p 2 ), with p 1 and p 2 longitudinal and traverse momentum respectively, in the interaction region (figure 3) and far from it (not presented here) much before, during, and after the encounter. At the beginning of the simulation, f (p 1 ) and f (p 2 ) are perfectly Maxwellian everywhere in the simulation window ( figure 3(a)). At t = 913.5ω −1 p and around x 1 = 0, f (p 1 ) is much broader than f (p 2 ) ( figure 3(b)), while for x 1 0 or x 1 0, f(p 1 ) and f (p 2 ) remain equivalent. Electrons near x 1 = 0 have been heated longitudinally. This preferential heating is the result of the non-linear physics characterising the interaction between the two shocks. The combination of the two shock waves produces a central region with electrostatic potential higher than in the far downstream. Some of the downstream electrons have sufficient energy to overcome the potential barriers that mark the transition between downstream and upstream. Upon entering the region with higher electrostatic potential, electrons from either side experience an increase of their longitudinal absolute velocities. As a consequence, the electron distribution in this central region has a higher thermal spread in the longitudinal direction than in the perpendicular direction. We observe that both our simulation setup and our results are fundamentally different from [35]. In [35], wave breaking and electron trapping are at the origin of strong longitudinal electron heating and temperature anisotropy in the downstream of an electrostatic shock driven by the interpenetration of plasma shells moving towards each other with an equal and opposite drift velocity. As a result, Weibel generated fields are detected downstream of the shock. We stress that simulations of single electrostatic shocks triggered by the expansion of a dense plasma into a more rarefied one show no magnetic field evidence in the downstream region and lead to constant shock velocities [36,51,67]. Further in time, the distribution functions in what is now the downstream of the shocks become very similar, but their spread is bigger than in the far upstream ( figure 3(c)). The electron isotropization could then explain the observed slowdown of the shock waves. In the reference frame of the shock, the electrons are moving towards the shock with longitudinal velocity −v s and the magnetic field B 3 bends them in the x 2 direction, thus partitioning their kinetic energy among the two in-plane directions. As a result, the longitudinal component of their velocity decreases by 1/ √ 2, which translates in the same shock slowdown when moving back to the simulation reference frame (e.g. downstream reference frame). We computed the anisotropy Δ = α /α ⊥ − 1, where α ⊥ and α are connected to the effective electron kinetic temperatures T e⊥ and T e according to [68]: Here, the subscripts ⊥ and indicate the direction perpendicular and parallel to the instability wave vector and f (p 2 ⊥ , p ) is the relativistic generalisation of the bi-Maxwellian electron distribution function [68]: where is the normalisation constant, K n the modified Bessel function of the second kind of order n, γ = (1 + p 2 /(m 2 e c 2 )) 1/2 the electron Lorentz factor and p the momentum. We calculated T e⊥ and T e by performing a matrix diagonalisation of the temperature tensor T as described in [69]. The result of the latter operation is such that T 1,1 = T e⊥ and T 2,2 = T e identify the direction with the highest and lowest temperature, respectively. Figure 4(a) reports Δ in the simulation window at t = 974.4ω −1 p , when the anisotropy reaches its maximum Δ max = 0.7. Apart from the central region, the anisotropy in the rest of the domain is almost zero. Given α ⊥ 0.19, Δ = 0.7 and the wavenumber k = 2π/λ = 0.12ω p /c as provided by the simulation and solving the dispersion relation in [68] for these values of Δ and k, we have computed the growth rate of the electron Weibel instability Γ theory = 0.019ω p . We have compared the theoretical prediction with the numerical results. Figure 4(b) shows the evolution of the energy of B 3 in the region between x 1 = −30c/ω p and x 1 = 30c/ω p . The magnetic field energy has an exponential growth until t 1309ω −1 p , where it reaches a maximum and then starts decaying. The slope of the linear phase allows for inferring the growth rate of the instability, which is measured to be 0.02ω p in excellent agreement with the theory. This constitutes a strong evidence that the electron Weibel instability is at play and causes the magnetic field to grow, leading to the isotropisation of the electrons. We also considered the possibility that the magnetic field is caused by the current filamentation instability due to the counterstreaming flow of ions reflected by the shocks or by these ions drifting through the upstream plasma at rest. However, for the range of ion velocities observed in the simulation, the fastest growing mode would have a growth rate and wavelength one order of magnitude smaller and larger, respectively, than the measured growth rate and wavelength [28,29]. Furthermore, simplified simulations of proton beams similar to those reflected by one shock streaming in a plasma at rest did not show any growth in the magnetic field on time scales comparable with our counterpropagating shocks. We notice that the simulation wavenumber and the relative growth rate do not correspond to the fastest growing mode. We conjecture this is because the fastest mode quickly grows and saturates; hence, the only modes surviving and observable are those at lower ks. The evolution of the Fourier transform of B 2 3 (figure 4(c)) displays indeed low energetic modes at k 0.25ω p /c (wave number corresponding to the fastest growing mode according to [68]), which saturate  in a time interval of ≈100ω −1 p after the collision. In this short time interval, the magnetic field in the transverse out-of-plane direction shows a larger number of filaments at the instability onset, which rapidly merge to produce the topology shown in figure 1(i). We verified that this is not the result of a numerical artefact via a thorough convergence study and we remark that the transverse size of our simulation window and resolution are such that they support 0.013ω p /c k 125.63ω p /c.
It is thus the Weibel instability that causes the shocks to slow down. The intense longitudinal electron heating occurring while the shocks are approaching produces a temperature anisotropy, which in turn drives the electron Weibel instability in the interaction region. The instability saturates when the distribution functions are isotropic again. In the 1D case, the Weibel instability, whose modes are perpendicular to x 1 , cannot occur and therefore the shock velocity is only slightly decreased by the longitudinal heating (see figure 2).

Three-dimensional simulations of shock collision
In order to confirm our findings in a more realistic geometry, we perform an additional 3D simulation employing the same setup detailed in section 2. For this simulation we considered a plasma composed of hot electrons with temperature T e = 1.5 MeV and cold ions. The initial density jump that drove the shock was σ = 10. A simulation box of size 1200 × 240 × 240(c/ω p ) 3 discretised with 4800 × 960 × 960 cells was used. The time step was chosen to satisfy the Courant condition. Each species was modelled using eight particles per cell and quartic interpolation. Also in this case, periodic boundary conditions were selected along all the directions, but the simulation box is large enough that they do not interfere with the plasma dynamics.
The expansion of the more dense lateral plasmas into the less dense central plasma leads to the generation of two counterpropagating electrostatic shocks. The shocks collide at the centre of the box at t 974.69ω −1 p . Their velocity is observed to decrease after the collision, as shown in figure 5. Specifically, before the interaction the shocks move towards each other at v s = 0.059c, while after the collision, their velocity is reduced to v s = 0.018c. As in the 2D case, the slowdown is connected with the development of the Weibel instability, which causes the growth of an intense magnetic field in the region downstream of the shocks (see figure 6(a), which shows a cut of B 3 at x 3 = 120c/ωp or the supplementary material for a 3D visualisation). We observe that the magnitude of the generated magnetic field in 3D is comparable with the 2D results ( figure 6(b)).  Maximum amplitude of the magnetic field (a), instability growth rate (b), shock velocity before the collision (triangles) and shock velocity after the collision (asterisks) as a fraction of the initial shock velocity (c) versus T e for σ = 10 (black) and σ = 5 (red). The magnetic field amplitude and shock speed scales as √ T e , while the data for the growth rate are fitted with quadratic curves.

Experimental perspectives
The ability to drive two counterpropagating shocks experimentally, whose collision will lead to the growth of the Weibel instability, depends on the laser intensity and the target maximum density and thickness. These parameters affect the initial density jump σ and electron temperature T e , which trigger the shock formation and determine the shock initial speed, thus influencing the development of the instability. To understand the dependence of the maximum amplitude of the magnetic field and the growth rate of the instability, both measurable in experiments, on T e and σ, we have performed a parameter scan varying these quantities (figure 7). The field amplitude at saturation varies as √ T e (figure 7(a)): saturation of the field is expected when the magnetic pressure ∝B 2 equals the plasma pressure ∝T e . The growth rate values are fitted by quadratic curves (figure 7(b)); the quadratic trend can be explained by the larger electron inertia of more energetic electrons. Given the higher inertia, they are more difficult to isotropise and the time required for the Lorentz force to accomplish this, therefore, longer. It is then expected that for even higher temperature, not of interest here, Γ will decrease with increasing T e . Figure 7(c) reports the shock velocity before the interaction and the ratio between values of the shock velocity before and after the collision (v s,1 /v s,0 ) as measured in our simulations. The shock speed scales as √ T e [51]. In all cases analysed, we observe the shocks slowing down. This slowdown seems to be larger when the shock is triggered by higher values of σ.
To compare our findings with current laser experiments, we compute the laser normalised intensity using the ponderomotive scaling formula [70]. To reach the considered electron temperatures, laser pulses with normalised vector potential a 0 = k B T e /m e c 2 + 1 2 − 1 1-14 are necessary. Such intensities are already or will be soon available, creating an opportunity to probe the setup proposed here. The other crucial factor to drive an electrostatic shock is the formation of a sharp density variation in the plasma [50]. Sharp gradients comparable to those employed in our simulations can be easily obtained using thin targets whose density profile reaches a maximum close to the critical density [43,46,51]. The laser will be stopped around the critical density and its radiation pressure will contribute to further steepen the density, which will lead to the formation of the shock. Overdense targets will prevent an optimal heating of the plasma electrons, while underdense targets are transparent to the laser and the laser will penetrate through the target without creating any density pile-up [71]. Therefore in both of these cases, conditions for electrostatic shock formation will not be met. Based on these considerations, mid-infrared pulses (λ 0 = 10 μm, where λ 0 is the central wavelength) with intensities 1.3 × 10 16 -2.7 × 10 18 W cm −2 illuminating two symmetric gas jets with peak densities 10 19 -10 20 cm −3 will drive shocks adequate to explore their collision. If an intensity of 2.7 × 10 18 W cm −2 is probably out of reach for current CO 2 laser systems, experiments at laser intensities of 1-8 × 10 16 W cm −2 have already shown the generation of electrostatic shocks using gas jets [41][42][43]. For near infra-red lasers (λ 0 1 μm), the derived values of a 0 correspond to intensities between 1.3 × 10 18 -2.7 × 10 20 W cm −2 , which could be achieved at current laser facilities. Indeed electrostatic shocks have been driven with these lasers in highly pressurised gas jets [45,[47][48][49] or tailored solid targets [46]. In order to launch two counterpropagating shocks, two lasers or one laser split into two pulses using a technique similar to that described in [72] could be used to illuminate two frontal targets or the opposite sides of the same target. Our considerations were validated through a 2D simulation where two laser-driven counterpropagating shocks were modelled (not shown here). For this simulation, we used two laser pulses with intensity I = 9 × 10 16 W cm −2 , wavelength λ 0 = 10 μm, normalised vector potential a 0 = 8.55 −10 λ (μm)(I (W cm −2 )) 1/2 = 2.5 and duration at full width half maximum (FWHM) τ FWHM = 14 ps. The lasers impinge upon two near critical density targets, whose longitudinal density profile raises from 0 to 2.5n c , where n c = 10 19 cm −3 is the critical density for the described lasers, over a distance of 100c/ω p and then exponentially decays, similarly to that described in [50]. In this simulation we observed the shocks forming as a result of a combination of laser heating and steepening of the targets and then colliding. The velocity of the shocks decreased from ±0.12c before the collision to ±0.09c after the collision and a strong magnetic field with a filamentary structure appeared in the interaction region, as in our simulations without lasers.

Summary
We have performed fully kinetic simulations of the interaction between two electrostatic shocks. The resulting head-on collision was highly inelastic; the shocks slow down to up to 50% of their initial velocity. The decrease in shock velocity is due to a strong magnetic field in the direction perpendicular to the simulation plane produced by the Weibel instability. The temperature anisotropy, which drives the electron Weibel instability after the collision, is caused by intense longitudinal electron heating occurring upstream of the shocks, as they approach each other. The magnetic field starts decaying when the electron distributions are again isotropic. We have confirmed the magnetic field growth due to the Weibel instability for different initial plasma conditions attainable with available or near future lasers, thus suggesting that the physics of shock collision can be investigated in the laboratory. For example, the 10 μm CO 2 laser systems available at University of California at Los Angeles [73] and at the Brookhaven National Laboratory [74] with a 0 in the range 1.5-2.5 would easily allow for testing this setup with standard hydrogen gas jets. In principle, the scheme could also be tested with near-infrared lasers, as long as higher densities are used. Finally, we note that the proposed setup allows for the electron Weibel instability to be probed in a fully collisionless regime in the laboratory, similarly to what has previously been done for the ion Weibel instability [7,75,76] and would complement the experiments of Zhang et al, which explore the growth and saturation of the electron Weibel instability in a collisional non-relativistic regime [59,60].