Simulation study of the formation of a non-relativistic pair shock

TL;DR

This study uses particle-in-cell simulations to explore how non-relativistic pair shocks form through plasma instabilities, revealing the roles of electrostatic waves, phase space vortices, and the Weibel instability in shock development.

Contribution

It provides a detailed simulation-based analysis of shock formation mechanisms in pair plasma, highlighting the interplay of electrostatic and magnetowaves in non-relativistic regimes.

Findings

Electrostatic two-stream and oblique modes dominate initial wave growth.

Phase space vortices facilitate wave saturation and plasma heating.

Strong magnetowaves develop downstream despite electrostatic shock mediation.

Abstract

We examine with a particle-in-cell (PIC) simulation the collision of two equally dense clouds of cold pair plasma. The clouds interpenetrate until instabilities set in, which heat up the plasma and trigger the formation of a pair of shocks. The fastest-growing waves at the collision speed c/5 and low temperature are the electrostatic two-stream mode and the quasi-electrostatic oblique mode. Both waves grow and saturate via the formation of phase space vortices. The strong electric fields of these nonlinear plasma structures provide an efficient means of heating up and compressing the inflowing upstream leptons. The interaction of the hot leptons, which leak back into the upstream region, with the inflowing cool upstream leptons continuously drives electrostatic waves that mediate the shock. These waves heat up the inflowing upstream leptons primarily along the shock normal, which results in an anisotropic velocity distribution in the post-shock region. This distribution gives rise to the Weibel instability. Our simulation shows that even if the shock is mediated by quasi-electrostatic waves, strong magnetowaves will still develop in its downstream region.

Equations

$${\mu }_0{ϵ}_0\frac{\partial \mathbf{E}}{\partial t}=\nabla \times \mathbf{B}-{\mu }_0\mathbf{J},$$

$$\frac{\partial \mathbf{B}}{\partial t}=-\nabla \times \mathbf{E}.$$

$$\left({\omega }^2{ϵ}_{xx}-k_y^2{c}^2\right)\left({\omega }^2{ϵ}_{yy}-k_x^2{c}^2\right)-{\left({\omega }^2{ϵ}_{yx}+k_x{k}_y{c}^2\right)}^2=0.$$

$${ϵ}_{\alpha \beta }(\mathbf{k},\omega )={\delta }_{\alpha \beta }\left(1-\frac{{\omega }_p^2}{{\omega }^2}\right)+\frac{{\omega }_p^2}{{\omega }^2}\sum _j\int d^{3}p\frac{p_{\alpha }{p}_{\beta }\mathbf{k}\cdot \left(\frac{\partial f_j^0}{\partial \mathbf{p}}\right)}{m\omega -\mathbf{k}\cdot \mathbf{p}}.$$

$$\frac{{\delta }_{TS}}{{\omega }_p}=\frac{\sqrt{2}}{2}.$$

$$\frac{{\delta }_W}{{\omega }_p}=2{\beta }_0^{\prime }.$$

$$A=\frac{\int f_t(v_x,v_y){(v_x-p_0/2m_e)}^{2}d{v}_{x}d{v}_y}{\int f_t(v_x,v_y)v_y^{2}d{v}_{x}d{v}_y}-1.$$

Full text

1 Introduction

Compact objects like neutron stars or black holes that accrete material can emit relativistic jets. These jets are composed of electrons, positrons and ions. The emission of relativistic jets by microquasars (Fabian & Rees, 1979; Margon, 1984), which are stellar-size black holes that gather material from a companion star, and by some of the supermassive black holes in the centers of galaxies has been observed directly (Bridle & Perley, 1984). The fireball model attributes gamma-ray bursts (GRBs) to ultrarelativistic jets, which are emitted during strong supernovae. A direct observation of the ultrarelativistic jets that trigger the GRBs and occur at cosmological distances is not possible. Their existence can thus not be established unambiguously (Woosley & Bloom, 2006). However, observations of a mildly relativistic plasma outflow during the supernova 1998bw by Kulkarni et al. (1998) lend some support to the fireball model.

The efficiency, with which the accreting object can accelerate the jet plasma, is not constant in time. A variable plasma acceleration efficiency results in a spatially varying velocity profile of the jet plasma. Internal shocks can form at locations with a large velocity change and these shocks can constitute strong sources of electromagnetic radiation (Rees, 1978). The prompt emissions of gamma-ray bursts, which are associated with internal shocks in ultrarelativistic jets, are visible across cosmological distances and internal shocks should thus be sources of intense electromagnetic radiation.

The relativistic factors of the internal shocks in GRB jets are probably of the order of a few. A wide range of theoretical and numerical studies have addressed the collision of lepton clouds at relativistic speeds and the instabilities that sustain the shock and thermalize the plasma that crosses it (Kazimura et al., 1998; Medvedev & Loeb, 1999; Brainerd, 2000; Sakai et al., 2000; Silva et al., 2003; Haruki & Sakai, 2003; Jaroschek et al., 2004; Medvedev et al., 2005; Milosavljevic & Nakar, 2006; Chang et al., 2008; Stockem et al., 2008; Bret et al., 2008, 2013; Sironi & Giannios, 2014; Marcowith et al., 2016). Such shocks thermalize plasma via the magnetic fields that are driven by the filamentation instability of counter-streaming beams of charged particles, which is also known as the beam-Weibel instability.

This magnetic instability outgrows the competing electrostatic instabilities if two equally dense lepton clouds collide or interpenetrate at a relativistic speed. Electrostatic waves and instabilities can, however, not always be neglected and they exist even in lepton plasmas where positrons and electrons are equally dense. An external electric field accelerates electrons and positrons into opposite directions, which creates a current. Even if the initial electric field perturbation is removed, the current leads to a self-generated electric field in the plasma. The electric field eventually becomes strong enough to reverse the flow directions of electrons and positrons. The motion overshoots though and electrons and positrons oscillate around the equilibrium position. This oscillation is sustained by electrostatic fields, which can couple resonantly to other beams.

The jets of microquasars contain a significant fraction of positrons (Trigo et al., 2013) alike the jets that trigger GRBs. The high variability of a microquasar jet, which is emitted by a stellar-size black hole that accretes material from a companion star, suggests that internal shocks are present in such jets (Kaiser et al., 2000; Miller-Jones et al., 2005). Microquasar jets expand at a lower speed than GRB jets. Typical velocity changes are thus likely to be smaller and the internal shocks of microquasars may not always be relativistic, in which case electrostatic processes may become more important. Nonrelativistic pair shocks have so far not received much attention and the structure of their transition layers remains unknown.

We explore with a particle-in-cell (PIC) simulation the initial evolution phase of a leptonic shock that forms when two equally dense pair clouds collide at a speed that is one fifth of the speed of light. The electrons and positrons of each cloud have an equal number density and mean speed and the net charge and current of each cloud vanishes. The thermal spread of the particle velocities is small compared to the cloud collision speed and the instabilities can initially be described in the cold plasma limit. Both, the pair temperature and the collision speed are probably too low to be realistic for the plasma distribution close to an internal shock of the jet of a microquasar. We chose these initial conditions because electrostatic instabilities tend to be more important than magnetic ones for nonrelativistic collisions of cold plasma. We can thus study a shock, for which electrostatic effects are maximized in its transition layer and that thus constitutes a lower bound for pair shocks with respect to the degree of magnetization.

Our simulation shows that the transition layer of the pair shock is indeed mediated by nonlinear and predominantly electrostatic waves. However, the incomplete thermalization of the inflowing upstream plasma by the shock results in a downstream plasma with a thermally anisotropic distribution. This anisotropy is strong enough to trigger the growth of the Weibel instability in its original form (Weibel, 1959). Even a practically electrostatic pair shock thus generates a magnetic field in its downstream region. The thermal anisotropy and, hence, the growth rate of the Weibel instability can be increased further by the presence of ions and by the ambipolar electric field they drive at the shock front (Stockem et al., 2014).

Our paper is structured as follows. The shock formation mechanism, the PIC simulation method, our initial conditions and the expected spectrum of growing waves are discussed in Section 2. Section 3 presents the simulation results and section 4 is the summary.

2 Shock formation, the simulation code and the initial conditions

2.1 The formation mechanism of a collisionless leptonic shock

We examine the formation of shocks out of the collision of two charge- and current-neutral clouds of electrons and positrons. The plasma we consider is initially unmagnetized, no ions are present and all lepton species have the same temperature. The absent binary collisions imply that both lepton clouds will move through each other until plasma instabilities start to grow. Only three wave modes can develop for our initial conditions. The two-stream modes are purely electrostatic and their wave vector is aligned with the collision direction. The quasi-electrostatic oblique modes have a wave vector that is oriented obliquely to the collision direction and they belong to the same wave branch as the two-stream modes. The third mode is the filamentation mode, which is also known as the beam-Weibel mode (Califano et al., 1998). The wave vectors of these waves form the angle $\pi/2$ with the collision direction. These three modes grow simultaneously during the shock formation stage. Their growth is eventually halted by nonlinear processes, which heat up the plasma in the overlap layer and bring it closer to a thermal equilibrium.

A leptonic shock can be created in a PIC simulation by the collision of one lepton cloud with a reflecting wall. The reflected leptons move against the inflowing leptons that have not yet reached the wall and an overlap layer develops. The instabilities in this overlap layer let waves grow that heat up the plasma when they saturate. The expansion of the heated plasma is limited on one side by the wall and a shock forms on the other side. The shock evolution is resolved correctly once a downstream region has formed that is thick enough to decouple the shock from the wall. The formation phase of the shock may, however, not be resolved correctly by this computationally efficient method. The mechanism that triggers the filamentation or beam-Weibel instability is that particles with oppositely directed current vectors repel each other and particles with parallel current vectors attract each other. The instability saturates by forming current channels that collect particles with the same direction of the current vector. Current channels that contain particles with oppositely directed current vectors are separated by magnetic fields. A reflection of a particle by the wall changes its velocity component along the wall’s normal direction and, thus, the direction of its current vector. Spatially separated current channels can, however, not form at the wall because the particle is not spatially displaced by the reflection. The suppression of the filamentation instability at the reflecting wall will affect the spectrum of the unstable waves.

This spectrum is resolved correctly if we let two separate lepton clouds collide. If both clouds differ only in their mean speed, then we have to resolve in the simulation two identical shocks that enclose the expanding downstream region. It is computationally expensive and unnecessary to track both shocks for a long time. Here we let a long and a short lepton cloud collide. We increase the time interval during which we can observe the shock between the downstream region and the long lepton cloud. The second shock moves into the opposite direction and it eventually reaches the simulation boundary. By that time, it does no longer affect the evolution of the other shock.

2.2 The particle-in-cell (PIC) simulation method

We model the collision of the lepton clouds with a particle-in-cell (PIC) simulation. The PIC simulation code is based on the kinetic plasma model, which approximates each plasma species $i$ by a phase space density distribution $f_{i}(\bm{x},\bm{v},t)$. The position vector $\bm{x}$ and the velocity vector $\bm{v}$ are treated as independent coordinates, which allows for arbitrary velocity distributions at any given position. The number density of this species is the zero’th moment of the distribution $n_{i}(\bm{x},t)=\int f_{i}(\bm{x},\bm{v},t)\,d\bm{v}$ and the mean speed $\bar{\bm{v}}_{i}(\bm{x},t)=\int\bm{v}f_{i}(\bm{x},\bm{v},t)\,d\bm{v}$ corresponds to its first moment. The number density and the mean speed yield the charge density $\rho_{i}(\bm{x},t)=q_{i}n_{i}(\bm{x},t)$ and the current density $\bm{J}_{i}(\bm{x},t)=q_{i}\bar{\bm{v}}_{i}(\bm{x},t)n_{i}(\bm{x},t)$ of the species $i$. The total charge density $\rho(\bm{x},t)=\sum_{i}\rho_{i}(\bm{x},t)$ and current density $\bm{J}(\bm{x},t)=\sum_{i}\bm{J}_{i}(\bm{x},t)$ update the electromagnetic fields via a finite difference approximation of Ampere’s and Faraday’s laws on a numerical grid.

[TABLE]

[TABLE]

The EPOCH code (Arber et al., 2015) we use fulfills $\nabla\cdot\bm{B}=0$ and $\nabla\cdot\bm{E}=\rho/\epsilon_{0}$ to round-off precision.

An ensemble of computational particles (CPs) with the charge $q_{i}$ and mass $m_{i}$ approximates the phase space density distribution $f_{i}(\bm{x},\bm{v},t)$. The relativistic momentum $\bm{p}_{j}$ of the $j^{th}$ CP of species $i$ is updated via a discretized form of the relativistic Lorentz force equation $\frac{\partial}{\partial t}\bm{p}_{j}=q_{i}\left(\bm{E}(\bm{x}_{j})+\bm{v}_{j}\times\bm{B}(\bm{x}_{j})\right)$ and its position is updated via $\frac{\partial}{\partial t}\bm{x}_{j}=\bm{v}_{j}$. The electric field and the magnetic field are interpolated from the numerical grid to the particle’s position $\bm{x}_{j}$ to update its momentum. The current density on the grid, which is used to update the electromagnetic fields, is the sum over all particle currents after they have been interpolated from the particle positions to the grid nodes.

2.3 The simulation setup

Our two-dimensional simulation box has the length $L_{x}$ along $x$ and $L_{y}$ along $y$. The simulation box is subdivided into the two intervals $-0.65\,L_{x}<x<0$ and $0\leq x<0.35\,L_{x}$. The boundary at $0.35\,L_{x}$ is reflecting and that at $-0.65\,L_{x}$ is open. The boundary conditions at $y=0$ and $y=L_{y}$ are periodic. We place electrons and positrons with equal densities $n_{0}$ and temperatures $T_{0}$ = 10 eV everywhere in the box at $t=0$. The electrons and positrons in the interval with $x\geq 0$ have a vanishing mean speed. The electrons and positrons in the interval $x<0$ have the mean speed $v_{0}=0.2c$ along $x$. No new particles are introduced while the simulation is running and the simulation is stopped well before the end of the inflowing lepton cloud encounters the shock or before the leptons that are reflected by the boundary at $x=0.35\,L_{x}$ return to the shock.

We normalize the position to the electron skin depth $\lambda_{s}=c/\omega_{p}$, where $\omega_{p}={(n_{0}e^{2}/\epsilon_{0}m_{e})}^{1/2}$ is the electron plasma frequency of one cloud. Velocities are normalized to $c$. Momenta are normalized to $cm_{e}$ and we define $p_{0}=v_{0}m_{e}$ as the mean momentum of a lepton of the plasma cloud in the half-space $x<0$. The box size $L_{x}\times L_{y}=60\times 2.4$ is resolved by $1.9\times 10^{4}$ grid cells along $x$ and by 760 grid cells along $y$. Electrons and positrons are represented by 25 CPs per cell, respectively. The time is normalized to $\omega_{p}^{-1}$. The simulation time $t_{sim}=120$, which is subdivided into $57200$ equal time steps. We normalize the electric field to $\omega_{p}cm_{e}/e$ and the magnetic field to $\omega_{p}m_{e}/e$.

2.4 The solution of the linear dispersion relation

We have to verify that our box is large enough to resolve the competing unstable modes and we want to determine the wave mode, which grows fastest for the selected initial conditions. We solve for this purpose the linear dispersion relation in order to determine the spectrum of the growing waves. The solution is computed under the assumption that the overlap layer has an infinite size. This condition is approximately fulfilled if the colliding clouds can interpenetrate for some time before the instabilities grow.

The initial velocity spread for $T_{0}=10$ eV is about $v_{th0}=4.5\times 10^{-3}c$ and both clouds drift toward each other at $v_{0}=0.2c$. Thermal effects can be neglected for the ratio $v_{0}/v_{th0}=44$ and the lepton beams are cold. We solve the linear dispersion relation in the frame of reference in which the total momentum vanishes. The pair clouds move in this reference frame into opposite x-directions at the speed modulus $\beta_{0}^{\prime}\equiv v_{0}^{\prime}/c=1/10$.

The non-relativistic dispersion equation for a perturbation of the form $\exp(i\bm{k}\cdot\bm{r}-i\omega t)$ and a wave vector $\bm{k}$ with an arbitrary orientation is (Bret et al., 2010)

[TABLE]

where $\delta_{\alpha\beta}$ is the Kronecker symbol and

[TABLE]

The problem of finding the fastest growing mode has been solved (Bret & Deutsch, 2005; Bret et al., 2013) for cold distributions of the form $f^{0}_{j}(\bm{p})=\delta(p_{y})\delta(p_{x}-P_{j})$.

Figure 1 shows the solution of the linear dispersion relation for our plasma parameters.

The growth rate peaks at the wave number $k_{x}\lambda_{s}\approx 12$ and its value does not depend on $k_{y}\lambda_{s}$ for the considered wave number interval. The fastest-growing modes are thus the two-stream/oblique modes. Their peak exponential growth rate is

[TABLE]

The filamentation modes are characterized by a flow aligned component $k_{x}=0$. Califano et al. (1998) estimated their growth rate as

[TABLE]

Figure 1 demonstrates that the growth rate of the filamentation modes with $k_{x}=0$ is smaller than that of the two-stream/oblique modes, which confirms the aforementoned approximations since $\delta_{W}<\delta_{TS}$ for $\beta_{0}^{\prime}=1/10$.

We can estimate with the help of Fig. 1 if and how our limited box size will affect the spectrum of growing waves. The simulation employs periodic boundary conditions along $y$ and the box length is $L_{y}$ in this direction. The smallest resolved wave number is thus $k_{c}=2\pi/L_{y}$ or $k_{c}\lambda_{s}=2.6$ and waves with $k_{y}<k_{c}$ can not grow. Figure 1 shows that the growth rate of the filamentation modes decreases below $\delta_{W}$ for $k_{y}<k_{c}$ while that of the two-stream/oblique modes remains unchanged. The main effect of the limited box size along $L_{y}$ is thus to suppress the wave numbers where the growth of the filamentation instability is negligible. If our simulation shows that the plasma dynamics is governed by the two-stream/oblique modes, then we would obtain the same result also for larger $L_{y}$.

3 Simulation results

We discuss the simulation results at selected times and focus on the shock that forms at lower values of $x$. The first part addresses the wave modes that trigger the formation of shocks. The second part discusses the structure of the shock and the electric fields that mediate it and the final part examines the growth of magnetic fields.

3.1 Instability and nonlinear saturation

The two clouds of initially unmagnetized collisionless lepton plasma will move through each other for some time before plasma instabilities set in. Figure 2(a-c) displays the electric $E_{x}$ and $E_{y}$ components as well as the magnetic $B_{z}$ component at the time $t_{1}=7.6$.

The lepton cloud, which was initially located in the half-space $x<0$, has moved by $v_{0}t_{1}=1.5$ towards increasing values of x. Waves have grown in the cloud overlap layer, which spans the interval $0<x<1.5$ at this time. The distribution of $E_{x}$ reveals waves with a wavelength $\lambda\approx 0.4$. The $E_{y}$ and $B_{z}$ components are closely correlated and both oscillate rapidly along y.

The in-plane electric field components and the out-of-plane magnetic field at the time $t_{2}=14.7$ are displayed in Fig. 2(d-f). Figures 2(a) and (d) show the same distribution of $E_{x}$ except for the larger amplitude. Their spatial confinement demonstrates that these waves do not propagate along $y$. The wave structures belong to electrostatic two-stream modes. The patterns in $E_{y}$ resemble those in $B_{z}$ and their amplitude ratio is comparable to that at the earlier time. The spatial correlation of the field structures in the distribution of $E_{y}$ and $B_{z}$ suggests that they belong to the same waves.

We can extract some properties of the waves from a comparison of the amplitude of $E_{y}$ and $B_{z}$ at the times $t_{1}$ or $t_{2}$. The ratio of the field energy densities $\epsilon_{0}(E_{x}^{2}+E_{y}^{2})/2$ and $B_{z}^{2}/2\mu_{0}$ is in the given normalization $(E_{x}^{2}+E_{y}^{2})/B_{z}^{2}\approx 100$. The particles of both clouds move at a speed $\approx v_{0}/2$ relative to the waves, which are slow-moving in the reference frame of the overlap layer. The electric force imposed on a charged particle, which moves with $v_{0}/2$ = 0.1, is 50 times larger than the magnetic force. We conclude that the wave’s magnetic field neither has a significant energy density nor does it affect the lepton dynamics. The waves are thus quasi-electrostatic and their wavelength along the collision direction is $\approx 0.4$. The amplitude of the waves has increased by a factor $\approx 50$ during the time interval $t_{2}-t_{1}=7.1$. If we assume that the waves grow exponentially, then their growth rate is $\delta\approx 0.5$ in units of $\omega_{p}$, which matches that in Fig. 1.

The wave modes that yield the observed electric field can be identified with its spatial power spectrum. We Fourier-transform the in-plane electric field distribution $E_{p}(x,y)=E_{x}(x,y)+iE_{y}(x,y)$ over the spatial interval $0.2<x<2.7$ and over all $y$ and multiply it with its complex conjugate. Figures 3(a, b) show the power spectra at the times $t_{1}$ and $t_{2}$ in the quadrant $k_{x}>0$ and $k_{y}>0$.

The power spectrum at $t_{1}=7.6$ shows wave power at $k_{x}\lambda_{s}\approx 14$, which extends up to maximum perpendicular wave number $|k_{y}\lambda_{s}|\approx 200$. The wavenumber $k_{x}\lambda_{s}=14$ corresponds to a wavelength along $x$ of about $0.45$.

The flow-aligned wave number $k_{x}\lambda_{s}\approx 14$ of the fastest-growing waves and the extension of wave power to large values of $k_{y}$ agree with the numerical solution of the linear dispersion relation in Fig. 1. The solution of the linear dispersion relation predicts a peak growth rate that does not depend on the value of $k_{y}$ for the considered wave numbers. The wave spectrum on Fig. 3(a) does however suggest that waves with a low value of $k_{y}$ grow faster. The growth rate is proportional to the amplitude the wave would reach after a given time if its growth would not be limited by nonlinear effects. The electric field amplitude, which is necessary to form phase space vortices, decreases with increasing values of $k={|k_{x}^{2}+k_{y}^{2}|}^{1/2}$ (O’Neil, 1965) and the discrepancy between the spectral distribution in Fig. 3(a) and the solution of the linear dispersion relation in Fig. 3(c) can thus be explained with a stronger nonlinear damping that is imposed on the modes with a large value of $k_{y}$.

The power spectrum in Fig. 3(b) is still concentrated on the two-stream / oblique mode branch. Its width along $k_{y}$ has diminished, which suggests that thermal damping is at work; the range of wave numbers $k_{y}$ that are unstable to the oblique mode instability is large in a cold plasma, while the wave growth is concentrated at low values of $k_{y}$ if the plasma is hot (Silva et al., 2002). A first and second harmonic along $k_{x}$ have emerged. The wave amplitudes have thus reached a non-linear regime (Umeda et al., 2003).

Non-linear effects in the wave distribution should be tied to changes in the lepton distribution. Figure 4 shows the phase space density distributions $f(x,p_{x})$ of the electrons and positrons.

The overlap layer of both clouds spans the interval $-0.2<x<3$. The counterstreaming clouds have not yet merged along $p_{x}$. However, the substantial particle acceleration demonstrates that the instability is about to saturate. The density in the overlap layer is twice that of a single cloud and the density flucutations caused by the waves are of the order of $5\%-10\%$ (not shown).

The filamentation instability starts to grow immediately but two-stream outgrows it.

3.2 Shock formation

The two-stream instability saturates by forming stable phase space vortices in the electron and positron distributions (Berk & Roberts, 1967) and the same holds during the initial saturation stage of the oblique mode instability (Dieckmann et al., 2006b). Electron phase space vortices are characterized by strong bipolar pulses in the electric field distribution, which correspond to a localized positive excess charge. Positron phase space vortices correspond to a localized negative excess charge.

The in-plane electric field components at the time $t_{3}=58.8$ are displayed in Fig. 5(a,b).

The electric $E_{x}$ component shows such bipolar field structures. A large quasi-planar field pulse is located at $x\approx 2.5$ in the interval starting from $y\approx 1$ that goes through the periodic boundary at $y=2.4$ until $y\approx 0.5$. The polarity of $E_{x}$ indicates the presence of a positive excess charge in between both electric field bands. If this quasi-planar bipolar pulse is associated with an electron phase space vortex, then the latter should be detectable in the electron phase space density distribution even if it has been integrated over all values of $y$.

Figures 5(c, d) show the corresponding electron and positron distributions. Figure 5(c) confirms the existence of a phase space vortex in the electron distribution at this location. The vortex in Fig. 5(c) spans the spatial interval $2<x<3$ and the momentum interval $-1<p_{x}/p_{0}<1$. The mean momenta of the upstream electrons and positrons are modulated by the electrostatic potential of the vortex when they pass it, but they are not trapped by it. The upstream leptons continue to move to increasing values of $x$ until they are thermalized upon entering the downstream region $x>4$, which is characterized by a dense phase space density distribution between $0<p_{x}/p_{0}<1$. This thermalization can only be accomplished by the field structures seen in the in-plane electric field between $x\approx 3$ and $x\approx 4$ in Fig. 5(a, b).

The distribution of the positrons shows two smaller vortices that surround the large electron phase space vortex. The positron vertices are centered at $x\approx 1.8$ and $x\approx 3$. The zero-crossing of the electric $E_{x}$ component and, thus, the extremal point of the electrostatic potential at $x\approx 2.5$ in Fig. 5(a) corresponds to a stable equilibrium point for the trapped electrons. Hence it is an unstable equilibrium point for the positrons, explaining why the vortices of positrons and electrons are staggered along $x$.

A small localized cloud of electrons and positrons is centred at $x\approx 0$ and $p_{x}\approx 0$. The cloud is an artifact from our initial conditions. The finite growth time of the electrostatic instabilities implies that the waves start to grow well behind the front of the plasma cloud that was initially located in the half-space $x>0$. This charge- and current neutral cloud is stable against electrostatic instabilities, because its extent along $x$ is not sufficiently large to allow it to interact with the inflowing upstream leptons via a two-stream instability.

Figure 6 shows the in-plane electric field distribution and the associated lepton phase space density distributions at the time $t_{4}=120$.

We observe strong quasi-planar electric field structures in the $E_{x}$-distribution in the interval $-3<x<0$. Their amplitude is comparable to the one that gave rise to phase space vortices in the electron- and positron distributions at the earlier time. These electrostatic structures in $E_{x}$ have propagated well beyond the initial collision boundary $x=0$ reaching a position $x\approx-3$. We find relatively strong electric field oscillations in $E_{x}$ and $E_{y}$ between $0<x<8$. The transition layer of this shock thus spans at this time an interval with the width $\Delta x\approx 10$.

The strong planar waves in the interval $-3<x<0$ in Fig. 6 are correlated with phase space vortices in the hot lepton population at low speeds. The vortices of electrons and positrons are staggered along $x$. The electrons and positrons that gyrate in these vortices originate from the hot plasma component and they are well-separated along $p_{x}$ from the inflowing upstream leptons. The mean speed of these phase space vortices is less than $p_{x}=0$, which implies that they move towards decreasing values of $x$. The mean speed of the vortices decreases with an increasing distance from the shock transition layer and they are thus accelerated away from the downstream plasma. The leptons, which gyrate in the vortices, reach a peak momentum $\approx-p_{0}$.

The simultaneous presence in the interval $-6<x<5$ of the hot leptons that have leaked from the downstream region and the cooler drifting upstream leptons implies that the overall plasma distribution is non-thermal and thus unstable. The electric field of the phase space vortices seeds the instability and we observe momentum oscillations along $p_{x}$ in the cool inflowing electrons and positrons that increase with $x$ in the interval $-3<x<-1$. The oppositely directed oscillations of electrons and positrons result in a strong current, which induces an electric field. The electric field oscillates in space and its oscillation amplitude decreases in unison with the net current in the direction of decreasing values of $x$. We can describe this oscillation in terms of a product between a sinusoidally oscillating electric field and an envelope function.

A spatially varying envelope function gives rise to a ponderomotive force (Kono et al., 1980) that does not depend on the sign of the particle charge and accelerates electrons and positrons in the direction of decreasing values of the envelope function. This ponderomotive force is excerted by the modulated upstream plasma onto the hot leptons that form the phase space vortices and it accelerates them.

The leptons in the interval $1<x<6$ in Fig. 6(c,d) are composed of a hot dilute component and the cool dense upstream leptons. Both populations gradually mix and they merge to a single one at $x\approx 6$. We observe electric fields in this interval in Fig. 6. These fields show some piecewise planar structures, which correspond to phase space vortices with a limited extent along $y$. The two strongest localized structures at $x\approx 1$ are separated by a perpendicular $E_{y}$ field at $y\approx 0.8$. These localized structures are likely to be the result of an instability of initially planar phase space vortices or phase space tubes. Indeed, two-dimensional PIC simulations (Oppenheim et al., 1999) of phase space tubes in a stabilizing magnetic field show that the phase space tubes gradually break up along their axes. The collapse of a phase space vortex is an effective way to scatter the leptons in phase space, which results in a mixing of the downstream and upstream leptons in the $x,p_{x}$ plane. Strong small-scale electric fields are observed up to $x\approx 8$. The absence of phase space vortices with $x>8$ demonstrates that the lepton distribution in this interval is no longer unstable to electrostatic instabilities.

The phase space density distribution of the leptons was uniform in the interval $4<x<6$ at the time $t_{3}=144$ and in the interval $7<x<13$ at the time $t_{4}=297$. We conclude that these intervals correspond to a downstream region that is close to being in a thermal equilibrium, at least with respect to electrostatic waves and instabilities. The density distribution along $x$ sheds light onto how much the plasma is compressed by the shock crossing. Figure 7 compares the density distributions at the times $t_{3}$ and $t_{4}$.

The density converges at low $x$ to the initial density. A density peak is observed close to $x\approx 0$ in Fig. 7(a), which corresponds to the dense lepton cloud at this position shown in Fig. 5(a). The density peak has disappeared in Fig. 7(b). The electric fields that grew in this spatial interval in response to the instability between the inflowing upstream leptons and the leaked hot downstream leptons teared the cloud apart by accelerating the electrons and positrons into opposite directions. The density rises from about $N(x)=1$ to $N(x)\approx 3$ over a few electron skin depths. The plasma compression factor of about 3 is the one expected for a strong nonrelativistic shock (Zel’Dovich et al., 1967).

3.3 Secondary instabilities and magnetic field generation

The electric fields associated with the phase space vortices heat up the leptons via Landau damping (Landau, 1946; O’Neil, 1965) and their collapse scatters them in phase space. The effects of this heating on the lepton distribution is visualized by Fig. 8, which shows the phase space density as a function of $p_{x}$ and $p_{y}$ at three positions along $x$. The distribution has been integrated over $y$ and over an interval along $x$ of width 0.3.

The distribution in Fig. 8(a) has been sampled far upstream of the shock. The upstream leptons constitute the cold dense beam that is located at $p_{x}\approx p_{0}$. The leaked leptons form a hot and dilute beam that moves at $p_{x}\approx-p_{0}$. The mean speed of the hot lepton beam exceeds that expected from a specular reflection, since the shock is moving to increasing values of $x$. Figure 8(b) reveals that the inflowing upstream leptons have been heated up by the time they reach the position $x\approx 2$. They are distributed over a wider velocity range and their peak value of the phase space density has thus decreased. The temperature is of the order of 100 eV. These leptons are immersed in a hot dilute lepton component. Its thermal momentum spread is of the order of $p_{0}$ and the temperature is thus about one keV. The inflowing upstream leptons form a hot beam at $p_{x}\approx p_{0}$ in 8(c) that is only about twice as dense as the leptons in the hot population.

The waves observed close to $x=2$ in Fig. 6(a,b) suggest that the velocity distribution in Fig. 8(b) is still unstable to an electrostatic instability. It can not be the two-stream instability because that one requires two beams that are well-separated along $p_{x}$. This distribution can, however, still be unstable to the electron acoustic instability. Alike the well-known ion acoustic instability, which is driven by a drift between cold ions and hot electrons, the electron acoustic instability can develop if cold electrons drift relative to a hot electron species. Waves grow if the drift speed between the hot and the cold electron species exceeds several times the thermal spread of the cold electron species (Gary, 1987). This condition is fullfilled in Fig. 8(b). We note in this context that although the phase space density of the hot leptons is two orders of magnitude less than that of the inflowing upstream leptons their number density, which we obtain by integrating the phase space density along $p_{x}$, is of the same order. The interaction of counterstreaming lepton beams with a similar density results in rapidly growing instabilities. To the best of our knowledge the acoustic instability in pair plasma has not yet been explored. Here we can not unambigously show that it exists in pair plasma, because the electric field may also be the residual field of a phase space hole that formed previously.

The distribution in Fig. 8(c) appears to be stable against electrostatic instabilities since we do not observe significant electric field oscillations in the region $x>8$. The velocity distribution of the leptons in this region is, however, not a Maxwellian. Therefore the plasma contains free energy that can be released by a collisionless instability. The thermal anisotropy contained in the total lepton velocity distribution $f_{t}(v_{x},v_{y})$ and measured in the rest frame of the downstream plasma can be estimated as

[TABLE]

A value $A=0$ would say that the thermal energy in the $x$ direction equals that in the $y$ direction, which would imply that there is no thermal anisotropy. We obtain the value $A\approx 6$ from the data shown in Fig. 8(c). Such a large anisotropy value results in the Weibel instability in its original form (Weibel, 1959; Morse & Nielson, 1969).

Figure 9 confirms that a magnetic field has grown in the downstream region.

Strong magnetic fields with approximately the same peak amplitude are present at both times. The amplitude of the magnetic field exceeds that observed in Fig. 2(c, f) by a factor 3 and it equals that of the electric field in the given normalization in Fig. 6(a,b). Nevertheless, the magnetic force, which is acting on a lepton that is moving at the speed of $0.1$ relative to these field patches, will still be an order of magnitude weaker than the electric force in the shock transition layer.

The strong downstream magnetic fields are not correlated with electric field structures (See Figs. 5(a, b) and 6(a,b)) and they are thus not driven by an oblique mode instability. The Weibel instability drives magnetowaves with a negligible electric field and with the same magnetic field direction as the one observed here and this instability is thus compatible with the simulation data. The Weibel instability yields magnetic fields with an energy density that can reach up to 10% of the thermal energy for strong temperature anisotropies (Morse & Nielson, 1971; Stockem et al., 2009). The velocity spread in Fig. 8(c) corresponds to a lepton temperature of the order of 1 keV or $10^{7}$ K. The magnetic pressure of a field $B_{z}=0.03$ is thus a few percent of the cumulative thermal pressure of electrons and positrons and within the range that is accessible to the Weibel instability.

We can compare the maximum size of the magnetic patches in Fig. 9 to the gyroradius of an electron in that field. The magnetic field $\tilde{B}_{z}$ and the collision speed $\tilde{v}_{0}$ in physical units can be calculated from the normalized ones via $B_{z}=e\tilde{B}_{z}/m_{e}\omega_{p}$ and $\tilde{v}_{0}=v_{0}c$. The gyroradius of an electron that moves at the speed $v_{0}/2$ relative to the stationary perpendicular magnetic field, which is normalized to $\lambda_{s}=c/\omega_{p}$, is then $r_{g}/\lambda_{s}=v_{0}/2B_{z}$. Taking $v_{0}=0.2$ and $B_{z}=0.03$ gives $r_{g}\approx 3\lambda_{s}$, which is about three times the coherence scale of the largest magnetic field patches. The magnetic field patches are also not stationary on the time scale needed to perform a gyro-orbit. This time equals for a maximum amplitude $B_{max}=0.03$ in our normalization $2\pi/B_{max}\approx 200$, which exceeds the simulation time. The leptons can thus not complete a full gyro-orbit. The magnetic field will instead deflect leptons by a small angle that depends on where the lepton entered the patch and on how long it stayed inside the patch. The magnetic field will thus scatter the leptons of the directed beam in Fig. 8(c). The repeated scattering of the leptons will eventually thermalize their distribution.

4 Summary

We have examined the formation and the initial evolution of a non-relativistic leptonic shock. The shock was created by letting two spatially uniform clouds of equally dense electrons and positrons collide at a relative speed of 0.2c. The absence of binary collisions implied that both clouds initially interpenetrated and formed an overlap layer. The two-stream instability grew in this overlap layer. The nonlinear saturation of the growing waves heated up the plasma in the overlap layer and transformed it into a downstream region that was enclosed by two shocks. We followed the evolution of one of them.

Some of the hot downstream leptons escaped upstream and interacted with the inflowing upstream plasma. Nonlinear and predominantly electrostatic phase space vortices formed, which mediated the shock transition layer. These structures are unstable (Wharton et al., 1968; Berk & Roberts, 1967; Morse & Nielson, 1969) and their collapse scattered and heated the leptons and gave rise to electric field flucutations. The interaction of electric field fluctuations and charged particles has a similar effect as binary collisions between particles (Dum, 1978; Dieckmann et al., 2006a; Bale et al., 2002; Baalrud et al., 2009; Bret, 2015) with respect to the thermalization of the inflowing plasma and this interaction contributed to the particle heating by the shock.

A spatially confined region formed, in which the plasma density exceeded the cumulative density of the colliding clouds. We found the compression factor 3, which is expected for a strong shock. However, the electrostatic shock transition layer could not fully thermalize the upstream plasma that crossed it. The residual thermal anisotropy of the downstream plasma drove magnetowaves via the instability proposed by Weibel (1959) in its original form. The magnetic fields were weak and spatially inhomogeneous. Nevertheless, the leptons that would enter these field patches would experience a small-angle deflection by the magnetic field.

The magnitude of the deflection angle depends on the time the particle needs to cross the magnetic patch. Repeated deflections will thus randomize the particle paths. This randomization will result in a thermalization of the downstream plasma on a time scale that exceeds by far the one accessible to our simulation. We note that the Weibel instability grows already at low values of the thermal anisotropy (Morse & Nielson, 1969). Magnetic field patches will thus grow and scatter particles until the lepton population has reached a thermal equilibrium.

The motivation of our work has been to better understand the properties of the internal shocks of microquasar jets. The initial temperature, which we gave to the leptons, and the shock speed are, however, lower than the values we may find close to the internal shocks of microquasars jets. We chose these low values to make the waves and plasma structures in the shock transition layer quasi-electrostatic. We will study in future work larger non-relativistic collision speeds and consider the effects of a larger initial temperature and determine up to which values the shocks resemble the one we have examined here.

We will also study in more detail the spectrum of the unstable waves that can grow in a pair plasma that consists of a cool beam immersed in a hot background population. Our simulation showed the presence of electrostatic waves in such a plasma. We proposed that such waves could be driven by an instability, which is similar to the electron acoustic instability Gary (1987), and we have to determine if a similar instability exists in the pair plasma we considered here.

Acknowledgements: Grant ENE2013-45661-C2-1-P from the Ministerio de Educación y Ciencia, Spain and Grant PEII-2014-008-P from the Junta de Comunidades de Castilla-La Mancha. Computer time and support was provided by the HPC2N and by the Swedish National Infrastructure for Computing (SNIC) through the grant SNIC2015-1-305.