Change of a Weibel-type to an Alfvénic shock in pair plasma by upstream waves

We examine with particle-in-cell simulations how a parallel shock in pair plasma reacts to upstream waves, which are driven by escaping downstream particles. Initially, the shock is sustained in the two-dimensional simulation by a magnetic ﬁlamentation (beam-Weibel) instability. Escaping particles drive an electrostatic beam instability upstream. Modiﬁcations of the upstream plasma by these waves hardly affect the shock. In time, a decreasing density and an increasing temperature of the escaping particles quench the beam instability. A larger thermal energy along than perpendicular to the magnetic ﬁeld destabilizes the pair-Alfv (cid:2) en mode. In the rest frame of the upstream plasma, the group velocity of the growing pair-Alfv (cid:2) en waves is below that of the shock and the latter catches up with the waves. Accumulating pair-Alfv (cid:2) en waves gradually change the shock in the two-dimensional simulation from a Weibel-type shock into an Alfv (cid:2) enic shock with a Mach number that is about 6 for our initial conditions.


I. INTRODUCTION
Black hole microquasars are binary systems, in which a stellar-mass black hole accretes material from a nearby companion star and ejects a relativistic jet [1][2][3] . The interior of the jet is composed of electrons and positrons and an unknown fraction of ions. It is thought that the low number density and the high temperature of the jet material let the effects of binary collisions between particles be small compared to those of the electromagnetic fields, which are driven by the ensemble of all particles. We call such a plasma collisionless.
A nonuniform flow speed of the jet can result in collisionless shocks at those locations where a faster flow is catching up with a slower one. A shock can also form between the relativistically moving jet material and the inner cocoon of the jet 4 , which is separated by a collisionless magnetic discontinuity from the ambient plasma into which the jet expands 5 . If the jet is leptonic then we would expect that its internal shocks slow down and heat a flow of electrons and positrons.
Shocks in collisionless plasma have a finite thickness. The shock transition layer is defined as the spatial interval where the inflowing upstream plasma is slowed down, heated and compressed by the electromagnetic fields that mediate the shock. The upstream material moves at a supersonic speed in the reference frame of the shock. Supersonic means in this context that the flow speed exceeds the speed of the wave that mediates the shock. Plasma, which crossed the transition layer, enters the downstream region. The downstream material is a hot and dense plasma in a thermal equilibrium that moves at a subsonic speed relative to the shock. a) mark.e.dieckmann@liu.se Collisionless leptonic shocks have been studied widely in the past with particle-in-cell (PIC) simulations. Relativistic collision speeds and a low temperature of the upstream plasma yield shock transition layers that are mediated by the filamentation instability, which is also known as the beam-Weibel instability [6][7][8][9][10][11] . The wavevectors k of these magnetowaves are oriented primarily orthogonally to the flow direction of the upstream plasma. The magnetowaves heat the upstream plasma, which crosses the transition layer, to a relativistic temperature before it enters the downstream region.
The exponential growth rate of the filamentation instability decreases as the collision speed is decreased. Oblique modes can outgrow the filamentation modes in particular if the interacting plasma beams have different densities. A nonrelativistic leptonic shock was examined in the simulation in Ref. 12 . Escaping energetic downstream particles interacted with the inflowing upstream plasma via an electrostatic oblique mode instability. The temperature of the preheated upstream plasma was higher along the collision direction than orthogonal to it after it crossed this layer. This thermal anisotropy resulted in a magneto-instability similar to the one found by Weibel 13 , which thermalized the pair plasma and established its thermal equilibrium. Increasing the shock speed to a mildly relativistic value triggered a filamentation instability between the intervals where electrostatic waves and the Weibel modes grew 14 .
Collisionless shocks in pair plasma, which is not permeated by a background magnetic field B 0 , are mediated by waves driven by the two-stream instability, the oblique mode instability and the filamentation instability 15 . Aligning B 0 with the shock normal modifies the spectrum of unstable waves 16 and only the twostream instability operates if B 0 = |B 0 | is sufficiently strong 17 . Such a magnetic field does not only modify the dispersion relation of the aforementioned modes. It can also introduce new unstable wave modes.
Consider for example a pair shock with a normal that is aligned with B 0 . This field can maintain different temperatures perpendicularly and parallel to B 0 . If the electromagnetic fields in the original shock transition layer cannot establish a thermally isotropic distribution of the plasma then instabilities like the mirror-or firehose instabilities [18][19][20] can grow behind the original transition layer and broaden it. Another example is provided by shocks in magnetized electron-ion plasma. Such shocks can emanate beams of particles with a super-Alfvénic speed into their upstream regions. Such beams can trigger the growth of Alfvénic waves upstream of the shock, thereby broadening its transition layer. The cosmic-ray driven Bell instability 21-23 falls into this category.
We study here with one-(1D) and two-dimensional (2D) PIC simulations the instabilities that grow in the transition layer of a shock in a pure pair plasma. A magnetic field is aligned with the shock normal. Its amplitude is not large enough to suppress the filamentation instability. The shock is created when the pair plasma, which is reflected at one boundary of the simulation box, interacts with the inflowing pair plasma. The two-stream instability grows first in the 1D simulation, which excludes geometrically the filamentation instability. It creates a plasma close to the reflecting boundary that is thermally anisotropic. Eventually a magneto-instability is triggered, which results in the growth of pair-Alfvén waves. In what follows, we refer to this instability as the pair-Alfvén wave instability. Pair-Alfvén waves mediate the shock in the 1D simulation. The filamentation instability outgrows two-stream instability in the 2D simulation and its magneto-waves sustain initially the shock. Some electrons and positrons outrun the shocks in both simulations. Initially these particles drive electrostatic instabilities upstream of the shock 14 . In time pair-Alfvén waves 24 grow in the upstream region of the shocks. The shock catches up with these slow waves and piles them up. Pair-Alfvén waves partially replace the filamentation modes in the 2D simulation as the means to sustain the shock. A similar replacement of the filamentation mode by Alfvén waves driven by Bell's instability has been observed at quasi-parallel electron-ion shocks 25 .
This work is the first, to our knowledge, to identify in simulation results, and to extensively study these pair-Alfvén-mediated parallel shocks and their associated upstream turbulence in plasmas with significant background magnetic fields. This work is therefore significant in bettering our understanding of trans-relativistic, non-collisional, pair-plasma environments, such as those predicted to occur at the base of relativistic jets in black hole microquasars.
Our paper is structured as follows. Section 2 summarizes the PIC method and it lists our initial conditions. It also gives a brief summary of the instabilities that develop in our simulation. Section 3 presents the simulation results, which are discussed in Section 4.

II. THE PIC CODE AND THE INITIAL CONDITIONS
Ampère's law µ 0 0Ė = ∇ × B − µ 0 J and Faraday's lawḂ = −∇ × E are approximated on a numerical grid, where E, B and J are the electric field, the magnetic field and the current density. The vacuum permittivity and permeability are 0 and µ 0 . Each plasma species j is approximated by computational particles (CP's). The i th CP has a charge-to-mass ratio q i /m i that must match that of the species j it represents. The electromagnetic fields are coupled to the CPs and the CPs are coupled to J via suitable numerical schemes as implemented in the EPOCH code we use 26 .
Our two-dimensional simulation resolves x by 2 × 10 4 grid cells and y by 2000 grid cells. Boundary conditions are reflective along x and periodic along y. We model one electron species and one positron species, which are uniformly distributed in space. Each species has the density n 0 /2 and is resolved by a total of 8 × 10 8 CPs, which corresponds to 20 CPs per cell for each species. We perform also a one-dimensional simulation with the same plasma parameters, where we resolve only x and use 10 7 CPs for each species. This number amounts to 500 CPs per cell for electrons and the same number for the positrons. We do not inject CPs while the simulation is running.
The plasma frequency is ω p = (e 2 n 0 /m e 0 ) 1/2 with e, m e being the elementary charge and the electron mass. We normalize the time t to ω −1 p , space to the plasma skin depth λ s = c/ω p (c : light speed) and densities to n 0 . Frequencies ω are normalized to ω p and wavenumbers k x in 1D or (k x , k y ) in 2D to λ −1 s . The simulation box spans the intervals 0 ≤ x ≤ 2650 and 0 ≤ y ≤ 265. Both species have a Maxwellian velocity distribution with the temperature T 0 = 10 keV, which gives the thermal speed v t = 4.2 × 10 7 m/s with v t = (k B T 0 /m e ) 1/2 (k B : Boltzmann constant), and the mean speed v b = −3v t (0.42c) along x. The pair cloud moves with the mean speed v b to the boundary at x = 0. Particles are reflected by this boundary and flow back to increasing x. A shock is triggered by the interaction of the reflected and inflowing pair plasma. The Debye length is λ D = 0.14. We set the grid cell size to λ D and hence we resolve one skin depth λ s by just over 7 grid cells. Electric fields E and magnetic fields B are normalized to m e cω p /e and m e ω p /e respectively. A magnetic field B 0 = (B 0 , 0, 0) with B 0 = 0.1 is present at t = 0. This value equals the normalized electron gyro-frequency ω c = eB 0 /m e ω p . Our pair plasma has a value β = n 0 k B T 0 /(B 2 0 /2µ 0 ) = 4. The pair Alfvén speed v A = B 0 /(µ 0 n 0 m e ) 1/2 equals 0.7v t (0.1c). Both simulations are stopped at t sim = 2500.
Pair-Alfvén waves play a pivotal role in sustaining our shocks. Thermal noise provides us with insight into their dispersion relation. Most of its power is concentrated on the undamped modes in the simulation plasma 27 . Figure 1 shows the frequency-wavenumber spectrum for a pair plasma with the aforementioned plasma parameters in its rest frame (the bar denotes the complex conju- gate). The thermal noise at low wave numbers peaks on the dispersion relation of the pair-Alfvén wave. Its phase The solution of the linear dispersion relation in cold pair plasma 24 follows closely ω = v A k x in this wavenumber interval and its phase speed hardly decreases for larger k x in the displayed interval. Its frequency converges to the cyclotron resonance ω c = 0.1 with increasing k x . The phase speed of the mode in the simulation decreases faster than that of the cold plasma mode for larger values of k x and the sharp noise peak disappears for k x > 0.25, which usually implies that the wave is damped (See also Ref. 20 ). This damping is likely to be caused by wave- Due to the presence of the background magnetic field, we have several wave modes in our simulation that can be rendered unstable by interacting pair beams and by a thermal anisotropy. These instabilities have been analysed under the assumption that the wave amplitudes are small and that the plasma can be approximated either by a bi-Maxwellian particle velocity distribution or by counter-streaming beams.
We consider first the case of counterstreaming beams. The large initial speed modulus |v b | = 3v t implies that the inflowing and reflected pair plasmas form two beams close to the reflecting boundary that are separated in velocity space. Both beams flow along B 0 . It is of interest to determine whether or not our value B 0 = 0.1 is large enough to suppress the magnetic filamentation instability of counter-streaming beams. Bret et al. 17 determine the value of B 0 that is needed to suppress the filamentation instability of counter-streaming cold electron beams. For our non-relativistic flow speed and assuming that the inflowing and reflected pair clouds consist only of electrons and are equally dense, we obtain the critical magnetic field value B c = √ 2v b /c giving B c = 0.6; we expect that our magnetic field is not strong enough to suppress the filamentation instability. This instability is, however, suppressed geometrically in our 1D simulation. We expect that the inflowing and reflected pair cloud trigger an electrostatic two-stream instability in the 1D simulation, which saturates by forming phase space holes 14,28 . Individual phase space holes are stable in 1D but unstable otherwise 29 . Their collapse heats the plasma. As long as the simulation resolves more than one spatial dimension, two-stream instabilities can mediate a narrow shock transition layer. This shock transition layer widens in a 1D simulation 14 because planar phase space holes can only thermalize by their slow coalescence 30 .
Let us assume that the initial instabilities have heated up the plasma along the collision direction. Weibel considered in his work 13 a single electron species with a bi-Maxwellian non-relativistic velocity distribution. Electrons had a lower temperature along one direction than along the other two. He found aperiodically growing waves with a wave vector along the cool direction. He also showed that aligning B 0 with this direction cannot stabilize the plasma. Weibel's work was extended to pair plasma 31 . Aligning B 0 with the cool direction of both species gives rise to two modes; Alfvén-like waves are positronic modes while the electrons give rise to magnetosonic-like waves. Both modes have an equal dispersion if electrons and positrons have equal distributions giving the combined mode a linear polarization 32 . We refer here to this mode as the pair-Alfvén mode. Gary et al. 31 consider first the case where the plasma temperature is the same in all directions. The pair-Alfvén mode is undamped in the wavenumber interval where the wave has the dispersion relation ω = v A k x , which is k x < 0.1 in Fig. 1. Increasing only the plasma temperature perpendicular to B 0 gives rise to a mirror-like instability.
Schlickeiser 18 examines also the case where the pair plasma is hotter along B 0 than perpendicular to it, which is relevant for our simulation. The initial value β = 4 in our simulation is boosted along B 0 because particles get reflected by the boundary at x = 0 and mix with the inflowing plasma. Typical particle speeds along this direction are thus increased from the thermal speed v t to the beam speed modulus |v b | = 3v t in the simulation frame. If we assume that the thermal pressure perpendicular to the magnetic field remains unchanged, we obtain a ther- : thermal speeds perpendicular and parallel to B 0 ). According to Fig. 8 in Ref. 18 , only the firehose instability can grow at low wavenumbers k x 1 if only wave vectors parallel to B 0 are taken into account. It yields aperiodically growing fluctuations. The Fourier spectrum of a wave with an amplitude that grows exponentially and non-oscillatory involves waves with a wide frequency spread. Aperiodically growing waves can thus couple their energy into the low-frequency pair-Alfvén mode 20 .

A. 1D simulation
Geometrical constraints suppress the filamentation instability in the 1D simulation because it grows by letting the plasma currents rearrange themselves in the direction orthogonal to the plasma collision direction. Pair-Alfvén waves and electrostatic Langmuir waves, which propagate along B 0 = (B 0 , 0, 0), can still grow. Figure 2 shows the time evolution of the plasma density, that of the electric field E x and that of the magnetic field components B y and B z . The B x component cannot change since ∇ · B = 0. The plasma density close to x = 0 increases to 2 after the simulation begins, due to the superposition of the inflowing and reflected plasma. We observe at early times density modulations, which are tied to oscillations of E x in Fig. 2(b). These waves are the result of a two-stream instability between the incoming and the reflected particles and they are present until t = 500 in an interval that expands with time.
Magnetowaves grow in Figs. 2(c, d) after the electrostatic waves collapse in the interval 0 ≤ x ≤ 150. Both magnetic field components show wave activity in the interval x ≤ 100 during 500 ≤ t ≤ 1000. Figure 2(a) shows that the plasma density close to the boundary increases from 2 to 3.5 during this time. The front of the dense plasma expands at the speed v f = 0.12c and it is correlated with strong waves. Some electrostatic waves propagate at a much larger speed.
Figures 2(c, d) reveal waves that are propagating from the upstream direction towards the front of the dense plasma after t = 500. These waves are transported with the upstream flow towards the dense plasma, then enter it and finally are damped out. The wave amplitude is almost stationary to the left of the white line, where the plasma has a low mean speed.
The frequency of the waves to the right of the white line is low in the rest frame of the upstream plasma. However, pair-Alfvén waves are the only low-frequency waves that can propagate along B 0 in the thermally isotropic upstream plasma 24 . Their wavelength 2π/k x ≈ 50 falls into the wavenumber interval where we find undamped pair-Alfvén waves in Fig. 1. The phase speed of the magnetowaves is −0.33c in the box frame and 0.09c in the rest frame of the upstream plasma, which has a mean speed −0.42c. The waves we observe to the right of the white line in Figs. 2(c, d) are thus pair-Alfvén waves that propagate in the upstream direction. They are connected to the waves to the left of the white line, which suggests that they belong to the same wave branch.
The phase space density distribution sheds light on why the collapse of the electrostatic waves in Fig. 2(b) at t ≈ 500 coincides with the growth of magnetowaves in Figs. 2(c, d). We select the electron distribution at the time t = 380 and show its projections onto the three momentum axes in Fig. 3. Figure 3(a) reveals phase space vortices in the interval 0 ≤ x ≤ 130. They are the product of a two-stream instability between the in-flowing and reflected electron beams. Dilute phase space vortices are present for 130 ≤ x ≤ 250. They are responsible for the fast structures in Fig. 2(b). The other two projections show an increase of the phase space density for 0 ≤ x ≤ 150 and some oscillations but the momentum range covered by the electrons is well below that in Fig.  3(a). The plasma in the overlap layer x ≤ 130 is hotter along x than along the other directions.
A pair-Alfvén wave instability lets the magnetic B y and B z components grow. These fields can deflect electrons and positrons from the parallel into the perpendicular direction and bring the plasma closer to thermal equilibrium; its density increases as observed in Figs 2(a). Pair-Alfvén waves are low-frequency waves, which explains why they do not grow immediately after the thermal anisotropy has developed. The pair plasma does not have a bi-Maxwellian distribution in the interval where the waves grow and the velocity distribution varies along x. The observed instability is thus not the firehose instability that was analyzed by Schlickeiser 18 . Both instabilities may, however, be related. Figure 4 displays the electron distribution close to the front of the dense plasma at the time t sim . It reveals that this front is a shock, which is located at x ≈ 370. It thermalizes rapidly the inflowing upstream electrons as they cross a transition layer in the interval 250 ≤ x ≤ 370. Our simulation box is at rest in the downstream frame of reference and v f = 0.12c thus corresponds to the shock speed in this frame. The phase speed of the pair-Alfvén wave in the downstream plasma is reduced to 3.5 due to the higher plasma density ≈ 3.5. The shock moves at the speed v f ≈ 2.3v * A . Pair-Alfvén waves downstream of the shock can not keep up with this and so the waves we observe must be generated directly behind the shock front by the thermal anisotropy. The growth of these waves is accelerated by the seed waves, which are transported with the upstream plasma to the shock. The lower Alfvén speed and the higher plasma temperature downstream of the shock imply that the waves can interact with the dense thermal bulk population of the particles, explaining the damping of the pair-Alfvén waves in this region. The oscillations of the mean velocities along y and z in the region x > 370 can be attributed to the pair-Alfvén waves that arrive at the shock from the upstream direction.
The large velocity oscillations of the upstream plasma in Fig. 4 may couple the pair-Alfvén waves, which are polarized along y, with those that are polarized along z. Such a coupling would break the linear polarization these waves have when their amplitude is low.

B. 2D simulation
We consider first the time evolution of the total plasma density, that of the field energy density P E = 0 (E 2 x + E 2 y ) y /(2n 0 k B T 0 ) of the in-plane electric field and that of the field energy densities P Bx = B 2 x y /(2µ 0 n 0 k B T 0 ), P By = B 2 y y /(2µ 0 n 0 k B T 0 ) and P Bz = B 2 z y /(2µ 0 n 0 k B T 0 ). We integrate these quantities over y as indicated by the subscript of the brackets. Figure 5 shows the box-averaged plasma density and the square roots of the field energy densities. Furthermore we give in Figs. 6-8 snapshots of the plasma at three representative times; early time before pair-Alfvn growth (t = 380), when the pair-Alfvn mode has saturated (t = 1100), and late time when transient effects have ended (t = 2500). Figure 5(a) shows that a shock front rapidly forms close to the wall in the 2D simulation. We find already at t ≈ 100 a boundary that separates the downstream plasma with mean density ≈ 3.5 from the upstream plasma. This shock front propagates at the speed v f 2 = 0.16c in the positive x direction, approximately 4/3 the speed of the 1D shock. This difference can be accounted for by the increased efficiency of particle-wave scattering arising from the filamentation-driven turbulence in the 2D case. This efficiency lets the shock establish itself quickly and its propagation speed is set by the pressure of a thermal downstream plasma. The shock forms later in the 1D simulation and its downstream region is initially not in a thermal equilibrium. We measure the shock speed v f 2 in the box frame. Its speed in the rest frame of the upstream plasma is v f 2 + |v b | = 0.58c, which exceeds the pair-Alfvén speed v A by the factor 6.
A dilute plasma beam outruns the shock in Fig. 5(a) and reaches an x-position of 700 at t ≈ 1700, which amounts to the speed 0.4c. This speed matches the mean speed modulus of the upstream plasma. The front of the dilute beam is initially trailed by a second beam with the density ≈ 1.3, which reaches x ≈ 300 at t ≈ 1300 and is absorbed by the shock at later times. Figure 5(b) shows a broad electrostatic pulse that outruns the shock that is centered on the front of the dilute plasma beam. We can understand the nature of these structures with the help of the spatial distributions of the electromagnetic fields and the phase space density distributions of the electrons and positrons. Figure 6 shows these at time t = 380. The density jump caused by the shock is located at x = 50. A downstream region with the density 3.5 has formed in Fig. 6(a), which is separated by a sharp shock boundary from the inflowing upstream region. Figure 6(b) reveals strong electrostatic waves just ahead of the shock. On average, their wave vector is aligned with the x-direction, which means that these structures can be resolved by a 1D simulation. We thus identify them with the phase space vortices we found in Fig. 3(a) in the interval 150 ≤ x ≤ 250. The width of these phase space vortices was on the order of few plasma skin depths, which matches the wavelength of the waves in Fig. 6(b). The mean speed of the vortices is larger than the shock speed, explaining why they outrun the shock in Fig. 5(a). The fastest phase space vortices involve only a very dilute plasma, which is not resolved by the color scale in Fig. 5(a). Hence the electrostatic waves seem to outrun the front of the dilute plasma. Figures 6(c-e) show the spatial distributions of the magnetic field components. The strongest component is B z . Initally the filamentation sets up a periodic oscillation of plasma density along y. This then drives turbulence with B z. This asymmetry between the y and z directions is, however, an artefact of the dimensionality of the simulation, in full 3D the filamentation instability would similarly drive magnetic turbulence along y. The filamentation-driven turbulence initially dominates since its growth rate is proportional to ω p whereas the growth rate of the pair-Alfvn mode is proportional to the electron cyclotron frequency ω c which is smaller by a factor of 10. Hence the fastest growth is seen in B z which saturates near the shock front almost instantly. Figures 5(c,  d) show that the peak of B z energy density is colocated with that of B x but not with that of B y . The modulation of B x is thus a consequence of the modulation of B z and presumably needed to fulfill ∇ · B = 0. The magnetic B y component shows oscillations in the downstream region of the shock, which suggests that pair-Alfvén waves have also grown in the 2D simulation. Their amplitude is barely a third of that of the filamentation modes. At first glance, one may thus conclude that the filamentation The electron phase space density distributions fe(x, px), fe(x, py) and fe(x, pz), which have been integrated over y, are shown in (f-h), respectively. They are normalized to the peak density in the upstream region and the color scale is 10-logarithmic. modes mediate the shock at this time. However, resonant interactions between particles and the pair-Alfvén wave could enhance their scattering. It has been reported that such resonant interactions can play an important role in the transition layer of Alfvénic shocks in electron-ion plasma 33,34 and the same may hold in our simulation. We also have to take into account that the magnetic amplitude of the pair-Alfvén wave is comparable to B 0 = 0.1. The wave is thus already in a non-linear regime. Figure 6(f) reveals a bulk distribution of the plasma centered on p x ≈ −0.5m e c for x > 50; this is the upstream plasma. The mean velocity jumps to p x = 0 at the location x = 50 of the shock. A larger momentum spread of the particles downstream of their shock implies that the plasma has been heated by crossing the shock. A dilute hot plasma beam with p x > 0 outruns the shock. The phase space vortices in this beam cannot be seen due to the phase space density distribution being integrated over y. The beam also drives the electrostatic waves in Fig. 6(b), a two-stream instability. Such waves typically have speed comparable to that of the beam with the lower plasma frequency (in this case the dilute beam reflected from the wall). Here the waves have velocity ≈ 0.27c in the simulation frame, and so quickly outrun the shock. Figures 6(g, h) demonstrate that the plasma has been heated by the shock passage also along the other directions. Both momentum distributions are evidence of energetic particles ahead of the shock. Their peak momentum decreases with distance from the shock and the distribution of the energetic particles joins with the distribution of the upstream plasma at x = 200. This implies that the energetic beam at x > 50 and p x > 0 in Fig.  6(f) consists of particles with relatively low values of p y and p z . The beam is thus fed by energetic downstream particles with a momentum vector that is almost aligned with B 0 . Figure 7 shows the same distributions at time t = 1100. The plasma density distribution in Fig. 7(a) reveals a shock at x ≈ 160. Its transition layer has not broadened along x and it still compresses the plasma to a downstream density ≈ 3.5. The electrostatic waves in Fig. 7(b) have propagated ahead of the shock as already demonstrated by Fig. 5(b). These waves apparently form only at early times. The two-stream instability grows quickly if two dense beams interact while remaining wellseparated along the velocity direction. This is the case here prior to the formation of the shock because the incoming plasma and the plasma reflected by the wall can interact. We find two well separated beams in Fig. 6(f) for x > 180. Once the shock has formed, the plasma that makes it upstream is hotter and less dense, and we no longer have a bimodal distribution in p x , hence the two-stream instability is quenched.
The amplitudes of B y and B z are now comparable. The B z component shows small patches, which oscillate rapidly along y, and are correlated with the oscillations of B x . The wave vector of the oscillations of B y is practically aligned with x; the magnetic field oscillations are thus tied to pair-Alfvén modes. We find oscillations of B y with the same orientation to the left and right of the shock at x = 160. We observe waves to both sides of the shock also in Fig. 7(e). They appear more turbulent. The oscillations of B y and B z appear identical to those in the 1D simulation. We thus infer from their different distributions in Fig. 7 that the oscillations of B y are caused by pair-Alfvén waves while those of B z are tied to filamentation modes and pair-Alfvén modes. The growth rate of the filamentation instability is greatest for particles streaming at relativistic speeds. These modes can then be driven even by dilute relativistic beams. This explains why we find the growth of waves in Fig. 5(e) in the B z distribution well ahead of the x−interval in which the pair-Alfvén modes grow in Fig. 5(d).
Figure 7(f) shows that the particles escaping upstream of the shock no longer form a beam that could drive electrostatic instabilities. The existing waves propagate ahead of the shock. The momentum distributions along p y and p z show oscillations ahead of the shock. Upstream particles are deflected by the magneto-waves close to and ahead of the shock. The large oscillation amplitude implies that the particle speed is large relative to the wave speed. This suggests in turn that the speed of the pair-Alfvén waves changes as they approach the shock. The piled-up waves form a shock precursor. Figure 8 shows the equivalent spatial distributions at the time t = 2500. The shock is now located at x = 400 in Fig. 8(a) and the uniform density of the downstream plasma suggests that it has been fully thermalized. The electrostatic waves have moved far upstream at this time and are no longer captured by the resolved x-interval in Fig. 8(b). Electric field oscillations are seen close to the shock boundary; they are not necessarily electrostatic since we find time-varying magnetic fields that can induce electric fields via Faraday's law. The magnetic B y and B z components show oscillations with a wave vector that is aligned with x. This suggests that pair-Alfvén waves are now mediating the shock. Figures 8(f-h) show a thermalized downstream region that is separated by the shock from the upstream plasma. The mean velocity along x of the plasma changes drastically across the shock and so does the plasma temperature. The inflowing upstream particles still gyrate in the magnetic field of the pair-Alfvén modes. Figure 9 zooms in on the shock front. The shock front in Fig. 9(a) is sharp even on this spatial scale. The B x component shows rapid oscillations along y, which suggests that the filamentation instability has not completely disappeared. Figures 9(c, d) demonstrate that waves at the shock front are quasi-planar. The B z component shows oscillations along y ahead of the shock and also at the shock.
How can we explain that initially the filamentation instability dominated while we observe predominantly pair-Alfvén modes at later times? The fact that we observed the filamentation instability at early times indicates that its waves have a larger linear growth rate meaning that they grow faster from noise levels to their nonlinear saturation 35 . The filamentation modes thermalized the pair plasma before the pair-Alfvén wave instability could grow. The growth of waves upstream of the shock, which are convected to the shock, implies that instabilities at the shock no longer grow from noise levels. Oscillations in the B y component can only be caused by pair-Alfvén waves. The upstream waves provide a seed for the instability at the shock, which causes the growth of strong pair-Alfvén modes at the shock. The upstream waves in the B z distribution in Fig. 9(d) are composed of filamentation modes and pair-Alfvén modes, which implies that both instabilities at the shock are boosted by these seed perturbations.
Finally we want to test how well the plasma has been thermalized by the shock crossing and what temperature it reached. We subdivide the downstream regions in Fig.  8(f-g) into the intervals 0 ≤ x ≤ 130, 130 ≤ x ≤ 260 and 260 ≤ x ≤ 390 and integrate them separately over x. We integrate the three momentum distributions separately and obtain a total of 9 momentum distributions. We plot in Fig. 10 all distributions into the same panel. All distributions agree well and can be approximated by a Maxwellian with the temperature 43 keV. We find no thermal anisotropies and no variations of the temperature with x; the downstream plasma has thus been well thermalized by the shock crossing. The electron phase space density distributions fe(x, px), fe(x, py) and fe(x, pz), which have been integrated over y, are shown in (f-h), respectively. They are normalized to the peak density in the upstream region and the color scale is 10-logarithmic.
FIG. 9. As in 7 this figure shows the plasma and field distributions at the simulation end time t = tsim = 2500, but zoomed in around x = 400 to show detail at the shock front. Panel (a) shows the plasma density. The magnetic Bx, By and Bz components are depicted by panels (b-d), respectively.

IV. DISCUSSION
In summary, we have demonstrated the existence of a novel type of collisionless shock which may occur in mildly-relativistic pair plasmas, such as those thought to exist in the inner regions of jets that are emitted by microquasars. These shocks are initially fed by turbulence generated by the filamentation instability as it was found in previous PIC simulations of unmagnetized leptonic shocks [6][7][8][9][10] . We aligned a magnetic field with the average shock normal. The magnetic field amplitude was not large enough to suppress the filamentation instability as in the simulation by Hededal et al. 16 . Instead it provided a new unstable wave branch, namely the pair-Alfvén wave. We compared the results of a 1D PIC simulation study, where the filamentation instability was suppressed due to the alignment of the simulation box with the magnetic field, with that of a 2D study.
Initially, both simulations provided different results. The instability that mixed the inflowing pair plasma with the one reflected by the simulation box boundary was an electrostatic two-stream instability in the 1D case and the filamentation instability in the 2D simulation. Electrostatic waves with an electric field pointing along the magnetic field can mix particles only along the magnetic field, while particles can be deflected in all directions by the magnetic filamentation modes together with the background magnetic field. The plasma downstream of the shock was thus far from a thermal equilibrium in the 1D simulation while it was almost thermal in the 2D one. Eventually, a pair-Alfvén wave instability thermalized the downstream plasma also in the 1D simulation.
Energetic downstream particles were able to outrun the forming shock before it settled into its final state. The particles formed a beam that drove electrostatic instabilities upstream of the forming shock. These fast waves outran the shock and they damped out eventually. Once the shock was established, the escaping particles formed a diffuse energetic upstream population that was no longer able to drive electrostatic waves. The leaking downstream particles increased the mean thermal energy of the upstream plasma in the direction of B 0 without introducing separate particle beams. Such a distribution is close to a bi-Maxwellian. The pair-Alfvén wave instability that developed ahead of the shock may thus be similar to the firehose instability discussed by Schlickeiser 18 . The pair-Alfvén waves it seeded were slower than the electrostatic ones and the shock could catch up with them. Their pile-up changed the shock into an Alfvénic one in the 1D simulation. Pair-Alfvén modes coexisted with filamentation modes in the transition layer of the two-dimensional shock. Both modes had comparable amplitudes of the out-of-plane magnetic field. Filamentation modes with in-plane magnetic fields are excluded in a 2D geometry and the in-plane magnetic field was tied exclusively to pair-Alfvén modes. We note in this context that the pair-Alfvén wave in a pair plasma with identical distributions of electrons and positrons has a linear polarization. Pair-Alfvén waves with an in-plane magnetic field may thus be decoupled from waves with a magnetic field that points out of the simulation plane. Future 3D PIC simulation studies have to test if the pair-Alfvén wave can replace the filamentation mode in a realistic geom-etry. Future work will also have to determine how the shock structure changes with the amplitude and direction of the background magnetic field. Finally one also has to test if pair-Alfvén waves keep their linear polarization in the non-linear regime.
Since potential sites for these shocks are ubiquitous in the high-energy sky we expect that the radiative signal produced by the accelerated leptons should be observable, at least in regimes with low optical depth to the base of the jet, e.g. radio. To achieve this, further work should examine more closely the expected radiation signature and provide a prediction for observations. Additionally we may expect to see significant acceleration for ions at these shocks 36 , and so it should be investigated if this can be shown to occur by adding a low number density "cosmic-ray" plasma species to future simulations. Furthermore the cosmic-ray instability and back-reaction of the accelerated ions on the turbulence may materially affect the shock structure. This was recently demonstrated for electron-ion shocks 25 but not for electron-positron-ion shocks. We defer this to a future publication.