Cocoon formation by a mildly relativistic pair jet in unmagnetized collisionless electron-proton plasma

Cocoon formation by a mildly relativistic pair jet

in unmagnetized collisionless electron-proton

plasma

Mark E Dieckmann, Gianluca Sarri, Doris Folini, Rolf Walder and Marco Borghesi

The self-archived postprint version of this journal article is available at Linköping

University Institutional Repository (DiVA):

http://urn.kb.se/resolve?urn=urn:nbn:se:liu:diva-153199

N.B.: When citing this work, cite the original publication.

Dieckmann, M. E, Sarri, G., Folini, D., Walder, R., Borghesi, M., (2018), Cocoon formation by a mildly relativistic pair jet in unmagnetized collisionless electron-proton plasma, Physics of Plasmas, 25(11), 112903. https://doi.org/10.1063/1.5050599

Original publication available at:

https://doi.org/10.1063/1.5050599

Copyright: AIP Publishing

electron-proton plasma

M. E. Dieckmann,1 _{G. Sarri,}2 _{D. Folini,}3 _{R. Walder,}3 _{and M. Borghesi}2

1_{Department of Science and Technology (ITN), Link¨opings University,}

Campus Norrk¨oping, SE-60174 Norrk¨oping, Sweden

2

Centre for Plasma Physics (CPP), Queen’s University Belfast, BT7 1NN, UK

3_{Ecole Normale Sup´´ erieure, Lyon, CRAL, UMR CNRS 5574, Universit´e de Lyon, France}

(Dated: October 27, 2018)

By modelling the expansion of a cloud of electrons and positrons with the temperature 400 keV that propagates at the mean speed 0.9c (c : speed of light) through an initially unmagnetized electron-proton plasma with a particle-in-cell (PIC) simulation, we find a mechanism that collimates the pair cloud into a jet. A filamentation (beam-Weibel) instability develops. Its magnetic field collimates the positrons and drives an electrostatic shock into the electron-proton plasma. The magnetic field acts as a discontinuity that separates the protons of the shocked ambient plasma, known as the outer cocoon, from the jet’s interior region. The outer cocoon expands at the speed 0.15c along the jet axis and at 0.03c perpendicularly to it. The filamentation instability converts the jet’s directed flow energy into magnetic energy in the inner cocoon. The magnetic discontinuity cannot separate the ambient electrons from the jet electrons. Both species rapidly mix and become indistinguishable. The spatial distribution of the positive charge carriers is in agreement with the distributions of the ambient material and the jet material predicted by a hydrodynamic model apart from a dilute positronic outflow that is accelerated by the electromagnetic field at the jet’s head.

PACS numbers: 52.65.Rr,52.72.+v,52.27.Ep

INTRODUCTION

Accreting black holes emit jets, which are composed of pairs of electrons and positrons and an unknown fraction of ions [1]. Their velocity can be moderately relativistic in the case of the microquasars [2–4]. Supermassive black holes in the centers of active galactic nuclei [5] can ac-celerate the jet plasma to ultrarelativistic speeds. Some of these jets are collimated, which allows them to cross astronomical distances at a relativistic speed. The fast-moving inner part of the jet remains separated from the surrounding material, which implies that the jet has an internal structure. Their internal structure was discussed for example by [6–10].

A hydrodynamic model was proposed in Ref. [6], which described the structure of a relativistic pair jet that prop-agated into ambient material with a larger mass density.

In the case of microquasars the latter can be the

inter-stellar medium (ISM) [11] or inter-stellar wind. Reference [6]

assumed that the jet is cylindrically symmetric. The jet expelled the ambient material in its way and carved out a path along which the jet material could move freely. Their model considered collimated and uncollimated jets. We focus on the collimated jets, which form if the interac-tion between the jet and the ambient material is strong. Reference [6] considered a jet with a planar head that propagated along the axial jet direction and had a nor-mal that is aligned with the cylinder axis. The front of the head was the forward shock between the pristine am-bient material and the shocked amam-bient material. The reverse shock on the rear side of the head slowed down

the material of the pair jet that flowed towards the head. The slowed-down material entered the inner cocoon. A contact discontinuity separated the inner cocoon from the outer cocoon, which consisted of the shocked ambi-ent material. Its high thermal pressure led to a lateral expansion of the shocked material on both sides of the contact discontinuity, which then flowed around the jet. A contact discontinuity separated the outer cocoon from the inner cocoon also on the sides of the jet.

Hydrodynamic jet models assume that the collisional-ity in the plasma of the jet and in its surroundings is large enough to establish and sustain thin shocks and discontinuities on the spatio-temporal scales of interest.

A positron with the energy 1 MeV is slowed down to a nonrelativistic speed on a distance of the order of kilo-parsecs by its collisional interaction with the interstellar medium if the latter has a number density of the order of

one particle per cm3[12]. Collisions may thus not be able

to sustain a contact discontinuity between the jet plasma and the ambient plasma on the scale of a microquasar jet, which is considerably smaller than that.

It is important to determine to what degree the model proposed in Ref. [6] is valid also for pair jets, for which the average time between particle collisions is large com-pared to the growth time of plasma instabilities. Al-though the hydrodynamic model will still be valid on the global scale of the jet, electromagnetic instabilities and structures will shape the shocks and discontinuities.

Indeed we know that electromagnetic fields are present in jets. Spectral properties of the synchrotron emissions of electrons and positrons suggest that the jets are per-meated by a magnetic field that is coherent on a large

scale with superimposed fluctuations [14]. Three sources have been proposed for them.

Firstly, coherent magnetic fields exist at the base of the jet close to the black hole, where they and the radiation are strong enough to generate huge clouds of electrons and positrons [15]. Simulations [16] show how this field extracts energy from the accretion disk and accelerates and collimates the jet outflow with this energy. The jet carries the magnetic field with it.

Secondly, the magnetic field of the interstellar medium (ISM) [11]or of the stellar wind of the companion staris compressed as it crosses the jet’s external shock, which results in an incoherent downstream magnetic field [17]. Thirdly, incoherent magnetic fields are generated by nonthermal plasma distributions in the pair plasma of the jet. Separate plasma populations can move through each other during a time that is short compared to the characteristic time between Coulomb collisions of parti-cles. The interpenetration of the jet plasma with the ISM or with plasma clouds within the jet that move at a different speed gives rise to anisotropic particle velocity distributions and to charged particle beams. They relax via the Weibel [18, 19] or the filamentation instability

also known as the beam-Weibel instability [20–29] that

drive strong magnetic fields that are coherent on small scales. Particle-in-cell simulations of cylindrical plasma clouds, which consist of cool electrons and positrons, that propagate through an ambient medium show that beam-Weibel instabilities, mushroom instabilities and kinetic Kelvin-Helmholtz instabilities (See Ref. [13] and refer-ences cited therein) develop in and close to the jet.

The coherent large-scale magnetic fields and probably also the small-scale magnetic fields, which result from plasma instabilities, play an important role especially on kinetic scales and they need to be studied further.

We examine here how a cloud of electrons and positrons interacts with an ambient plasma, which con-sists of electrons and protons. We select values for the cloud’s density, temperature and mean speed for which the ensuing instabilities accelerate protons [30] and gen-erate strong magnetic fields [31]. In these works, the ambient plasma and the pair cloud were uniform in the direction orthogonal to the expansion direction. Here we consider a pair cloud with a density that decreases quadratically with the distance from its central axis un-til it reaches zero. The cloud is truncated at its front.

We show with a simulation using the relativistic and electromagnetic particle-in-cell (PIC) code EPOCH [32] how the electrons and protons in the ambient plasma compela part of the pair cloudto form a magnetized jet that expels the protons in its path. The magnetic field, which is driven by a filamentation instability, acts as a

pistonthat piles up the protons in the lateral direction;

an outer cocoon forms, which is separated by an electro-static shock from the pristine ambient plasma. The shock propagates at about 3% of the light speed c. A flat head

of the jet deflects the ambient protons around it into the outer cocoon. The lateral width of the head is compara-ble to the wavelength of the filamentation instability and it propagates at the speed 0.15c.

The magnetic field of the piston is sustained by a fil-amentation instability, which develops in a small spatial interval around the piston that contains both ambient

plasma and particles of the pair cloud. The instability

converts the directed flow energy of thepair cloud into magnetic energy and heat and this interval is thus the in-ner cocoon. Our pair plasma has a thermal momentum spread that is comparable to its drift speed and a shock between the jet plasma and the inner cocoon would be weak. We could not detect one in the simulation.

An electric field is induced around the magnetic pis-ton. It accelerates the positrons that flow out of the jet. It is also responsible for the deflection of ambient protons around the jet’s head. It draws the electrons of the am-bient plasma into the jet, where they mix with the jet’s electrons. The magnetic piston induces an electric field in the outer cocoon, which pulls jet electrons into it.

The magnetic piston, together with the electric fields around it, distributes the carriers of positive charge into a form that is practically identical to that in Ref. [6] apart from the dilute outflow of energetic positrons at the jet’s head. We thus observe a jet in a collisionless plasma on a microscopic scale that resembles its macro-scopic hydrodynamic counterpart in a collisional medium even though the mechanisms, which shape the jet and its internal structure, are quite different.

Our paper is structured as follows. Section 2 discusses the numerical scheme of our PIC code and the initial conditions of the aforementioned simulation, which we refer to as the main simulation. The filamentation insta-bility between the pair cloud, which moves at a mildly relativistic speed and has a mildly relativistic tempera-ture, and the ambient plasma plays an important role in the formation of the jet in our main simulation.

Sec-tion 3 shows how the filamentaSec-tioninstability grows and

saturates in a reduced geometry and for the initial condi-tions of the plasma close to the symmetry axis of the pair cloud in the main simulation. Section 3 also discusses how the boundary conditions affect this instability. Sec-tion 4 shows how the jet grows out of the interaction

between the pair cloud and the ambient plasma in the

main simulation. Its structure and its evolution are ex-plained in terms of the filamentation instability. Section 5 summarizes our results.

THE CODE AND THE INITIAL CONDITIONS OF THE MAIN SIMULATION

A PIC code represents the electric field E and the magnetic field B on a numerical grid and evolves them with discretized forms of Amp`ere’s law and Faraday’s

law. The PIC code EPOCH we use solves Gauss’ law and the magnetic divergence law to round-off precision. Amp`ere’s law requires the plasma current J in order to update E in time.

PIC codes are based on the kinetic equations and they represent the phase space density distribution of each plasma species by an ensemble of computational parti-cles (CPs). Each particle carries with it a charge and a mass and the charge-to-mass ratio must match that of the plasma species it represents. The current contribution of each CP is deposited on the numerical grid and the global sum over all current contributions yields J, which goes into Amp`ere’s law. After the update of the fields, their values are interpolated to the positions of individual CPs and their momentum is updated with a discretized form of the relativistic Lorentz force equation.

Our main simulation resolves x and y and the jet

we observe will propagate along y. A spatially

uni-form ambient plasma at rest fills the simulation box at the time t = 0. The density and temperature of its electrons with the mass me and protons with the mass mp = 1836me are n0 and T0 = 2 keV. The plasma fre-quency ωp = (n0e2/ϵ0me)

1/2

(e, ϵ0 : elementary charge and vacuum permittivity) and c define the electron skin depth λs= c/ωp that we use to normalize space.

The temperature T0 exceeds that of the ISM. How-ever, the ISM will be heated up by radiation from an ap-proaching jet. Choosing a high temperature for the am-bient plasma yields a large Debye length λD = vth,e/ωp (vth,e = (kBT0/me)

1/2

, kB : electron thermal speed and Boltzmann constant), which allows us to resolve space with larger grid cells and hence reduce the computational cost. The ambient plasma will be heated to a tempera-ture≫ T0 and we do not expect that a reduction of T0 would affect significantly the jet evolution.

We normalize B to ωpme/e, E to ωpmec/e and den-sities to n0. The simulation box resolves the interval 0≤ x ≤ Lxwith Lx= 600 and−Ly/3≤ y ≤ 2Ly/3 with Ly = 3400 using a grid with 3500× 20000 grid cells. The electrons and protons are resolved by 28 CPs per cell, re-spectively. A pair cloud is superimposed on the ambient plasma. The densities of its electrons and positrons are nc(x, y) = 4− (x/94)

2

for x ≤ 188 and y ≤ 0 and zero otherwise as shown in Fig. 1(a). Each species has the temperature 400 keV and is resolved by 8.4× 108 _CPs. The mean speed of the electrons and positrons along y is V0= 0.9c. The plasma is charge- and current neutral at t = 0 and we set E = 0 and B = 0 everywhere.

Figure 1(b) illustrates the structure of a hydrodynamic collimated jet close to the jet’s head in Ref. [6]. The jet consists of light material, that moves through material with a larger mass density. A collimation shock, which is not depicted in Fig. 1(b), has shocked the pair plasma close to the base of the jet.

Our simulation setup exploits the symmetry of the jet

0 1 0 2 -1000 nc 3 4 0 200 (a) X Y 1000 ₂₀₀₀ 400 600 -1.5 -1 -0.5 0 0.5 1 1.5 (b) X 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 Y Shocked Jet Inner Cocoon Outer Cocoon Contact discontinuity Head Ambient medium Reverse shock Forward shock

Figure 1. Panel(a) shows the initial density distributions of theelectrons and positrons of the pair cloud. The boundary condition at x = 0 is reflecting and forms the symmetry axis of the pair cloud. Panel (b) shows the structure of a relativistic and collimated pair jet in a hydrodynamic model. The dashed line shows the jet’s symmetry axis. The pair cloud in (a) and the jet in (b) move along the y-axis.

by slicing it along the dashed line in Fig. 1(b), which results in the distribution of the pair cloud shown in Fig. 1(a). We use reflecting boundary conditions, which cuts the computing time in half and helps us to identify the physical processes at work. Numerical artifacts intro-duced by the reflecting boundary at x = 0 are limited to an interval with a width λsclose to the boundaryas we

show in the next section. The boundary fixes the phase of

thefilamentation modes that grow closest to the

bound-ary. We stop the main simulation at ts = 3150 (unit : ω−1p ) when the fastest particles, which were reflected at the boundary y = 2Ly/3, return to the region of interest.

FILAMENTATION INSTABILITY

We discuss here the growth and saturation of the fil-amentation instability between the pair cloud and the

ambient plasmaand how it is affected by the boundary

conditions. We compare the results of one PIC simulation with periodic boundary conditions, which we denote as Simp, with Simrthat uses reflecting boundary conditions. We argue that the main conclusions of our paper are un-affected by the usage of reflecting boundary conditions, although some differences between both simulations can be spotted on close inspection.

We initialize the plasma in both simulations with the plasma parameters found in the main simulation close to the boundary x = 0 for values y ≤ 0. We resolve only the x-direction, which is the direction that is orthogonal to the expansion direction of the pair cloud. The length

86 of the simulation box is resolved by 500 grid cells.

All plasma species have a spatially uniform tempera-ture and density to start with. The ambient electrons and protons have the temperatureT0= 2 keVand density n0 and their mean speed is zero everywhere. They are rep-resented by 2000 particles per cell each. The electrons and positrons of the pair cloud have the density 4n0, the temperature200 T0 (400 keV)and they propagate with

the speed V0 = 0.9c along y. Each cloud species is re-solved by 3000 particles per cell. All fields are set to zero

at t = 0 and we evolve the instability during0≤ t ≤ 525

with 104 steps. This time interval is smaller by a factor

6 than that covered by the main simulation.

The pair cloud moves with V0 to increasing values of y and a filamentation instability can only displace par-ticles along x. Such a displacement results in a current contribution Jy(x) along y of each cloud species. A mag-netic field Bz(x) grows if the current contributions of both cloud species do not cancel each other out. An elec-trostatic field Ex(x) grows if the instability redistributes the plasma species such that space charge is created.

At t = 0 the CPs were placed at the same positions in Simr and Simp and the same sequence of random num-bers was used to initialize their momenta. We expect at least initially and far from the boundaries an identical evolution of the instability in both simulations. A com-parison of the field and particle data provided by both simulations reveals effects introduced by the boundary conditions.

Figure 2 compares the time evolution of Bz(x, t) and Ex(x, t) in both simulations as well as the phase space density distribution of the protons in the cut plane (x, vx) at the time tc = 200. Strong fields have developed in Figs. 2(a-d) after the initial exponential growth phase of

the instability finished at t ≈ 15. We can not observe

the fields during their initial exponential growth phase since their low amplitudes are not resolved well by a lin-ear color scale. We refer to Ref. [29] for a discussion of the growth and saturation of a similar instability be-tween counterstreaming pair beams. The protons do not show a strong reaction to the electromagnetic fields at this time (not shown). The field amplitudes remain

ap-proximately constant until t≈ 100 and they continue to

grow afterwards on a slower timescale.

Thephase of themagnetic fields Bz(x, t) in Figs. 2(a,

b) differ for the filament closest to x = 0 but not for the one at larger x, which indicates that effects due to the boundary remain localized for these non-propagating

waves and for the short time scales we consider. An

elec-tric field Ex(x, t) grows in Figs. 2(c, d) thathas the same

wavelength as the corresponding oscillation of Bz(x, t)

and is in antiphase with it. Broad electric field bands

surround the spatial intervalsx≈ 2 and x ≈ 9, where a drastic change with x of the proton’s mean velocity is ob-served in Figs. 2(e, f). An integration of the phase space density distributions of the protons along vxyields large density peaks at x =2 and 9in Fig. 2(g). These protons have been pushed aside by the electric field. They form

phase space vortices [33] in the x, vxplane at later times.

Figure 2(d) reveals the growth of an electric field in

the grid cell next to the boundary at x = 0 until t≈ tc.

Itremains constant after that. Such a band is not found in Fig. 2(c)and it is thus a consequence of the reflecting boundary. It remains localized and it has no strong ef-fect on the proton distribution in Fig. 2(f). Effects due to the reflecting boundary other than its adjustment of

Figure 2. The time evolution of Bz(x, t) of Simp is shown

in panel (a) and that of Simr is shown in (b). Panel (c)

shows Ex(x, t) of Simp as a function of time and (d) that of

Simr. A 10-logarithmic time scale is used in these 4 panels.

Panels (e) and (f) show the phase space density distribution

f (x, vx) of the protons in Simpand Simr at the time tc= 200

(red horizontal lines in (a, b)), respectively. All color scales are linear and those of the two figures in the same row are identical. Panel (g) shows the proton densities.

the phase of the filamentation modes are thus weak and limited to a narrow interval.

Figures 3 and 4 show the phase space density distri-butions of the leptons in Simp and Simr at t = tc. The projection of the phase space density distribution on the plane spanned by x and by the momentum component px is f (x, px), while the projection f (x, py) involves the mo-mentum component along y. The lepton distributions in both simulations are practically identicalfor x > 2. They differ close to the boundary, which confines the filament

in Simr but not that in Simp. The confinement of the

filament closest to x = 0 in Simr does not change

quali-tatively its lepton distributions. In what follows we only

discuss the distributions computed by Simr.

The ambient electrons in Fig. 4(a) cluster in spatial intervals where the proton density in Fig. 2(g) is low. Their mean momentum along y in Fig. 4(b) reaches a high value close to the boundaries and close to x = 6 and it decreases to zero at the boundaries of the phase space cloud of the ambient electrons.

Figure 3. Lepton distribution at the time tc= 200 in the

sim-ulation Simp. The upper row shows the distributions f (x, px)

(a) and f (x, py) (b) of the ambient electrons. The second

row shows the distributions f (x, px) (c) and f (x, py) (d) of

the cloud electrons. The bottom row shows the distributions

f (x, px) (e) and f (x, py) (f) of the positrons. All distributions are normalized to the maximum value in the top panel of the respective column and displayed on a linear color range.

Figure 4. Lepton distribution at the time tc= 200 in the

sim-ulation Simr. The upper row shows the distributions f (x, px)

(a) and f (x, py) (b) of the ambient electrons. The second

row shows the distributions f (x, px) (c) and f (x, py) (d) of

the cloud electrons. The bottom row shows the distributions

f (x, px) (e) and f (x, py) (f) of the positrons. All distributions are normalized to the maximum value in the top panel of the respective column and displayed on a linear color range.

Figures 4(c, d) demonstrate that the cloud electrons accumulate mainly in those spatial intervals with a large proton density. The apparent broadening of the distri-bution in Fig. 4(c) along vx at x ≈ 2 and x ≈ 9 is

at least partially caused by the increased density of the

cloud electrons. The distribution of the cloud electrons

in Fig. 4(d) reveals that their mean momentum along y is reduced to a non-relativistic valuein the intervals where

they accumulate. Some of the cloud electrons even

re-verse their momentum at x = 2 and 9, which implies that the current density of the cloud electrons is close to the Alfv´en limit (See page 139 in [34]). The positrons in Figs. 4(e, f) have been expelled from the spatial intervals with a high proton density.

Figures 4(b, d, f) show how electrons and positrons are distributed after the instability saturated in order to drive the current Jy(x, t) that sustains Bz(x, t). We discuss only the distribution close to x = 0 due to the periodicity of the wave. The positron density has a max-imum at x = 0, while the density of the jet electrons has a minimum at this position. Ambient electrons provide an additional negative current contribution. The mean momentum along y of the positrons increases and their density decreases as we go away from x = 0. The oppo-site is true for the electrons. The positron density close to x = 2 is low and their mean velocity is largest close to this position. The magnetic field is primarily sustained by a spatial separation of both cloud speciesand its po-larity in Fig. 2(b) indicates that the current due to the

cloud electrons dominates at x ≈ 2 while that of the

positron cloud is stronger at x≈ 0.

The mean momenta in Figs. 4(d, f) remain positive and mildly relativistic. Let us assume that both cloud species drift at the mean speed vD= (0, vD, 0) with vD> 0 and that this net drift and the magnetic field B = (0, 0, Bz(x, t)) drive the electric field ED≈ −vD×B that is observed in Figs. 2(d). The peak values of Bz(x, tc) and Ex(x, tc) in Fig. 2(b, d) are0.8and 0.5, which gives a resonable value vD ≈ 0.6c, and the phase of the drift electric field matches that observed in the simulations. The drift electric field is responsible for the compression of the protons and the filamentation instability is thus not purely magnetic. This has to be expectednot only because the center of momentum of the plasma is drifting

along y but alsobecause the counter-streaming beams are

not symmetric [35] .

MAIN SIMULATION

Here we present the results of the main simulation with the setup we discussed in Section 2 and we show how parts of the pair cloud start to form a collisionless jet that moves along the boundary x = 0.

Figure 5 shows the density distributions of the elec-trons, positrons and protons and of the magnetic Bz com-ponentin a subinterval of the simulation boxat the final simulation time ts. Figures 5(a, b) reveal high-density bands in the distributions of the electrons and protons with the density≈ 3 that start at x ≈ 60 and y ≈ 700 and end at x≈ 10 and y ≈ 1150. No high-density band is present in the positron distribution in Fig. 5(c). The electron and proton densities are close to 1 for values of x that are larger than those of their high-density band and only few positrons are present.

Figure 5. The distributions of the plasma densities and of the magnetic Bzcomponent at the time ts: the electron distribution

is shown in (a), the proton distribution in (b), the positron distribution in (c) and that of Bz in (d). The color scale in (d) is

clamped to values between -0.35 and 0.35. The data has been smoothed with a filter that averages the quantity over 5× 5 cells (Multimedia view).

A magnetic band with the amplitude -0.35 follows in Fig. 5(d) this high-density band and is located at lower values of x. It extends to y < 700 and maintains its distance x ≈ 30 from the boundary at x = 0 for all 450≤ y ≤ 700. Oscillatory magnetic fields are observed at x≈ 200 and y < 450. These fields mark the front of the electrons and positrons that expand thermally along x. The protons inside the magnetic band have been evacuated from the two intervals 500 ≤ y ≤ 700 and 800 ≤ y ≤ 1150 and were replaced by positrons with a density between 0.7 and 1. Two density peaks are ob-served at x < 20 for 600≤ y ≤ 800 in Fig. 5(b). These structures are composed of electrons and protons.

The distributions of the plasma species and of the mag-netic field in Fig. 5 for 500≤ y ≤ 1200 resemble a jet, which is propagating to increasing values of y and is cen-tred on the axis x = 0. In what follows we will refer to the dominant structures in this y-interval with the terms that are defined for the hydrodynamic jet in Fig. 1(b).

The high-density band in the proton distribution in Fig. 5(b) has no counterpart in the positron distribution and hence it corresponds to piled-up ambient plasma; it is the outer cocoon of the jet. The high-density band ob-served in the electron distribution in Fig. 5(a) neutral-izes its charge. Protons are practically evacuated from the interval to the left of the outer cocoon. This inter-val is filled with pairs, which constitute the shocked jet material in Fig. 1(b).

The absence of binary collisions in the PIC simula-tion implies that contact discontinuities can not be sharp. The outer cocoon is trailed by the magnetic band in Fig. 5(d), which acts like a contact discontinuity. We refer to this discontinuity as the magnetic piston.

We will now address the processes that resulted in the

Figure 6. Evolution of the plasma density and magnetic field: The proton density distribution at t = 840 is shown in (a) and that of Bz in (b). The distributions of the proton density and Bz at t = 2000 are shown in panels (c) and (d) (Multimedia

view).

formation of the jet and interpret them in terms of the fil-amentation instability discussed in the previous section. Figure 6 examines the distributions of the proton den-sity and of Bzat the times t=840 and 2000. Figure 6(a) shows filaments with a thickness≈ 1 and separation ≈ 8, which are followed by magnetic stripes in Fig. 6(b). We focus on the stripe with the value Bz ≈ −1 at x ≈ 8 and 230 ≤ y ≤ 380. The magnetic fields in Fig. 6(b) and in Fig. 2(b) have a negative amplitude close to the boundary, which implies that a positive current is flowing along y between them and x = 0. The proton density in

the filament in Fig. 6(a) is low close to the boundary at x = 0, which matches the density profile in Fig. 2(g), and the peak densities are comparable.

The proton density close to the boundary in Fig. 6(a) is however lower over a broader spatial interval than that in Fig. 2(g). We understand this difference as follows. The size of the filaments along x is constrained in Simrby the small box sizeand by the presence of equally strong filaments at larger x, which is a consequence of the par-ticle distributions that were initially spatially uniform. Filaments can only grow via mergers in the 1D geometry [29]. The filaments can expand along x in Fig. 6 because

the pair cloud has a limited extent along x. Figures 6(c,

d) confirm that indeed the filaments grow in time. Figure 6(b) demonstrates that the magnetic field is confined to a narrower interval along x compared to that in Fig. 2(b). A filamentation instability and the cur-rents it drives is constrained to spatial intervals in which positrons and protons coexist since only there a filamen-tation instability between the ambient plasma and the cloud plasma can develop. The expulsion of protons from extended spatial intervals in Fig. 6(a) suppresses the

fil-amentation instability in these intervals. The

filamen-tation instability could not evacuate all protons in Fig. 2(g) and hencethe ambient plasma and the cloud plasma

could interactover a broad spatial interval.

The interval, in which the instability unfolds, has moved to larger y at t = 2000. Figure 6(c) shows two intervals with a proton density ∼ 0. The largest deple-tion is marked with F1. Its right boundary is charac-terized by a density band with the peak value ∼ 3 at x≈ 30 and a magnetic field Bz with a negative polarity at lower x in Fig. 6(d). The protons have been swept out by the magnetic field structure, which we identify as the magnetic piston in Fig. 5(d). A second similar-sized in-terval, where protons have been evacuated and replaced by positrons (not shown), is marked with F2 in Fig. 6(c). It demonstrates that the growth of the filament F1 along the boundary is not a numerical artifact. The filament F1 is the dominant one because the ram- and thermal pressures of the pair cloud are largest at low x. Indeed

Figs. 5(b, c) reveal smaller jets for example at x ≈ 200

and y≈ 250 where protons have been expelled.

Figure 7 shows the proton density, the magnetic Bz and the electric Ex component in a sub interval of the

simulation box close to the largest jet that flows along

the boundary x = 0 at the time t = ts. The other field

components are at noise levels. The proton density dis-tribution in Fig. 7(a) shows a low-density region at low x with a diagonal front for 800≤ y ≤ 1150. This region is bounded to the right by the outer cocoon, which has a thickness that increases approximately linearly along x with decreasing 800≤ y ≤ 1150.

The electric field distribution in Fig. 7(b) shows a narrow electric field band E1 with a positive amplitude and a weaker one E2 with a negative amplitude. Their

contours are overplotted in Fig. 7(a) and they match the borders of the outer cocoon. The band E2 follows the minimum of the magnetic field amplitude in the magnetic piston in Fig. 7(c) in the interval 800≤ y ≤ 1150.

The magnetic piston extends up to y ≈ 1180, which

marks the jet’s head, and it is separated from the reflect-ing boundary by an interval along x that is comparable to the size of a filament in Fig. 2(b). The magnetic piston does not lead to a sharp outer cocoon at the jet’s head. The observed steady increase of the proton density with increasing y across the head is in line with what we ex-pect from a hydrodynamic model. The latter states that the ambient material is piled up and deflected around the head by the forward shock so that it does not enter the jet’s interior. The density distributions of the elec-trons and posielec-trons in the shocked jet material do not reveal jumps and there is no sharp magnetic field band that could mediate a shock between the jet material and the inner cocoon; no well-defined shock has formed that bounds the inner cocoon.

The phase space density distribution of the protons in Fig. 8 reveals the cause of both electric field bands as well as of the proton density spikes in the shocked jet material, which we observed in in Fig. 5(b). The dense spiky structures at x≈ 0, y ≈ 700 and |v| ≤ 5 · 106 m/s are ion acoustic solitary waves [36]. These structures form when electrons and positrons stream over protons at rest [30]. Here they form in the protons that were located behind the magnetic piston when it formed and could therefore not be swept out by it.

The distribution at large x, y in Fig. 8 corresponds to the ambient protons with the temperature T0. The elec-tric field band E1 in Fig. 7(b) sustains an electrostatic shock that starts at y = 1150 and x≈ 10 and extends to increasing x and decreasing y. The overturn of the pro-ton distribution into the direction of large x, y at large |v| = (v2

x+ vy2) 1/2

is typical for a shock-reflected proton beam and a unique feature of collisionless shocks. The shock normal and the direction of the shock-reflected ion beam are almost parallel to the x-axis. We do not ob-serve a shock at the head of the jet, which explains the

slow change of the proton density in Fig. 7(a).

If the shock reflects the protons specularly then the peak value |v| ≈ 1.8 · 107 m/s implies a shock speed vs≈ 9 · 106m/s. Figure 2(f) demonstrates that the fila-mentation instability can indeed accelerate a large num-ber of protons to the speed vs. The ion acoustic speed in the ambient plasma is cs = (kB(γeT0+ γiT0)/mp)

1/2

with the adiabatic constants γe = 5/3 for electrons and γi = 3 for protons and cs ≈ 1.7 · 106 m/s. The Mach number of the shock is≈ 5.

Figure 8 does not show a shock that could be associ-ated with the electric field band E2 in Fig. 7(b). That electric field band coincides in space in Fig. 7(a) with a proton density gradient. Its polarity is such that it

ac-Figure 7. The density and field distributions at t = ts: (a) shows the proton density. The electric field component Exis shown

in panel (b) while (c) shows the magnetic Bz component. The contours of the electric field modulus|Ex| = 0.1 are overplotted

in red in (a,c). The band with Ex= 0.1 is E1 and that with Ex=−0.1 is E2.

Figure 8. The 10-logarithmic phase space density distribution of the protons with|v| = (vx2+ vy2)

1/2

at the time ts.

celerates protons to lower values of x and y, which would erode the density gradient. This is the ambipolar elec-tric field that is driven by thermal diffusion of electrons across a density gradient. The only field structure that is located close to E2 in Fig. 7(b) and strong enough to balance the electric field and, hence, counteract the ero-sion of the proton density gradient is the magnetic piston with its large amplitude for Bz.

Given the normalization of B and Bz≈ −0.35 we get the electron gyrofrequency ωce = eBz/me ≈ ωp/3. The proton gyro-frequency ωci= ωceme/mi≈ 1.8 · 10−4 and tsωci ≈ 0.6. The proton’s gyromotion is too slow to be responsible for the plasma’s pile-up.

A large value of Bz also corresponds to a large

mag-netic pressure. Let us assume that the outer cocoon is at rest and that its electrons and protons have the tem-perature T0 and density 3. Its thermal pressure is PT = 6n0kBT0 in SI units. The magnetic pressure is PB =

B2_/2µ

0 (µ0: vacuum permeability) in SI units. Only Bz grows in our simulation and we obtain from it the mag-netic pressure in SI units as PB = ω2pm2eBz2/(2µ0e2). We equate the thermal and magnetic pressures and arrive at the condition B2

z = 12kBT0/mec2 or |Bz| ≈ 0.22; the magnetic pressure due to Bz = −0.35 can balance the thermal pressure of the outer cocoon.

The shock at the front of the outer cocoon forms a straight diagonal line in Fig. 7(a). The shock in 7(a) forms a straight line between (x, y) = (10, 1100) and

(40, 900), which gives the line y = 1100− 6.5 · (x − 10).

We have previously estimated that the speed of the shock,

which bounds the outer cocoon, is about vs≈ 9·106m/s.

The shock propagates perpendicularly to the front of the outer cocoon in Fig. 7(a) and we assume for simplicity that the shock propagates along x.

We have determined the speed of the jet along y by comparing the proton density distribution along y and

across the external shock at x = 25 at the time 0.93ts

(not shown) with that at ts. The shock is located at

y≈ 980 at x = 25 in Fig. 7(a). We obtained a jet speed

≈ 5vs along y. Figure 8 shows that the jet’s head at

y ≈ 1150 and x < 10 launches the electrostatic shock,

which then propagates along x. If the external shock is launched by the head of the jet, which moves about 5 times faster along y than the shock moves along x, then we would expect a diagonal front of the outer cocoon with a slope that is comparable to that in Fig. 7(a).

We examine in more detail the mechanism that sus-tains the magnetic piston and how well it can separate the ambient plasma from the shocked jet plasma.

Figure 9 compares the density distributions of the am-bient electrons and the cloud particles with the contour line |Bz(x, y)| = 0.25. The contour line marks the posi-tion of the magnetic piston. The outer cocoon is bounded at low y by the magnetic piston. Figure 9(a) demon-strates that most of its electrons are ambient electrons. The ambient electrons also contribute significantly to the shocked jet plasma, which would not be possible in a hy-drodynamic model.

The electrons of the pair cloud in Fig. 9(b) are de-pleted close to x = 0 and their density is elevated in the outer cocoon with its increased proton density like in Figs. 4(c, d). Most positrons in Fig. 9(c) are confined to the left of the magnetic piston, which is in agreement with Figs. 4(e, f). We note in this context that only a small fraction of the positrons participated in the jet formation. Most propagated away from the region close to x = 0 before the jet formed. The confinement of the positrons by the magnetic piston is thus better than that suggested by Fig. 9(c).

We turn to the structure of the jet’s head, which is located in the interval 0 ≤ x ≤ 4 and y ≈ 1175. Its width along x matches the width of the positron cloud close to x = 0 in Fig. 4(e, f), which suggests that the jet’s head is a filamentation mode. The electromagnetic fields of the saturated filamentation mode had a strong effect on the particle momentum distributions in Fig. 4. Figure 10 shows the momentum distributions of the electrons and positrons in the phase space plane defined by y and py and the velocity distribution of the protons in the plane spanned by y and vy. We have integrated them across the interval 0≤ x ≤ 7.5, which includes the magnetic piston.

A comparison of the phase space density distributions of the ambient and cloud electrons along y and pyin Figs. 10(a, b) shows that both agree apart from the density change with y; the jet electrons are gradually replaced by ambient electrons with increasing y. Both species have the same mean momentum and momentum spread along the jet propagation direction. They have thermalized and form one population within the jet. The density of the combined electron distributions is comparable to that of the positrons (See also Fig. 9).

The positrons in Fig. 10(c) are hotter than the elec-trons up to y≈ 1100 and the fastest ones reach a Lorentz factor that is higher than that of the electrons by 2-3. The mean momentum of the positrons increases in the interval 1100 ≤ y ≤ 1200 where the head of the jet is located. Positrons, which have escaped from the jet via the head, reach a peak momentum≈ 8mec while that of the electrons is≈ 5mec.

Some positrons in the interval y≤ 103 _{stream back at} a speed that yields the Lorentz factor 3. We also find backstreaming electrons but they move at a lower speed. Jet particles, which can keep up with the jet, must move at least with the speed 0.15c along y. Backstreaming

par-ticles must have interacted with the jet’s electromagnetic fields. Backstreaming particles could only flow along the inner cocoon in a hydrodynamic model. This is not nec-essarily the case in a collisionless plasma, where no strong constraints exist for the shape of the phase space density distribution along the velocity direction. Figures 10(a-c) evidence backstreaming particles even at y≈ 800, where the integration interval 0 ≤ x ≤ 7.5 ends far from the magnetic piston. Here they mix with the material of the shocked jet in Fig. 1(b).

The acceleration of positrons in the interval 1100 ≤ y ≤ 1200 can only be accomplished by an electric field. The protons also react to this electric field. Those at y≈ 1125 are faster than those upstream with y > 1250. The mean velocity of the protons decreases with increas-ing y > 1125 while that of the positrons increases with increasing 1100≤ y ≤ 1200. This implies that the elec-tric field structure moves to increasing values of y at a speed that is small compared to that of the positrons. We had previously determined that the head of the jet moves at a speed 0.15 c, which fulfills this criterion. The electric field is thus tied to the jet’s head.

We attribute the electric field Ey > 0 to a lossy jet front. Amp`ere’s law is∇ × B = µ0J + µ0ϵ0∂E/∂t. The magnetic term and the current term will balance each other if no current is dissipated. However, current can be dissipated for example by beam instabilities or by the loss of positrons along the normal of the jet’s head. An electric field grows in this case which tries to increase the current to the value needed to balance ∇ × B. If the positive current to the left of the magnetic piston is dissipated away then an electric field Ey > 0 grows that will accelerate the jet positrons and the ambient protons in the way we observed in Fig. 10. An electric field Ey > 0 at the jet’s head will also accelerate ambient electrons into the jet. Figure 9 shows that indeed we find mostly ambient electrons close to the jet’s head.

Figure 10(d) shows that, apart from the slight acceler-ation at the head, the protons maintain their mean speed in the displayed y-interval. We do not observe an elec-trostatic shock at the front of the jet’s head. If there is no shock, then the proton density gradient across the head can only be explained by a proton flow around the head and into the outer cocoon. Such a proton deflection by the jet’s head matches the deflection of the ambient material by the head of a hydrodynamic jet.

The proton deflection can be accomplished in three ways. If the current, which sustains the piston, is dissi-pative also at its front then the electric field will deflect positive charges around the head and into the outer co-coon. A second contribution comes from the convectional electric field associated with the magnetic field of the pis-ton that moves at 0.15c into the ambient plasma ahead of the jet. A third contribution may also arise from the electric field, which is associated with the magnetic field of the filamentation mode (See. Fig. 2(b, d)).

Figure 9. The density distributions of the ambient electrons (a), of the electrons of the pair cloud (b) and of the positrons (c) at the time t = ts. The color scales are linear. The red curve denotes the contour|Bz(x, y)| = 0.25.

Figure 10. The phase space density distributions of the plasma species along the jet propagation direction at the time t = ts: The momentum distributions of the ambient

electrons (a), pair cloud electrons (b), and positrons (c) are normalized to the maximum value in (a) and displayed on a 10-logarithmic scale. Panel (d) shows the proton velocity distribution on a linear scale. All distributions have been in-tegrated over 0≤ x ≤ 7.5.

SUMMARY

We examined with a PIC simulation the expansion of a pair cloud into an electron-proton plasma. Their large temperature implied that the electrons and positrons of the cloud expanded rapidly, which decreased the den-sity of the cloud. A filamentation instability developed between the ambient plasma and the pair cloud in the

interval where the latter was still dense. This instability expelled the protons from large areas, which were then filled with positrons. Magnetic fields grew only in those locations where protons and rapidly streaming jet

par-ticleswere present, which confined the magnetic field to

small spatial intervals. We observed in the simplified one-dimensional study that the filamentation instability drives an electric field. The effect of the electric field is to push protons away from the positron filament. The instability and the magnetic field it drives follows the protons and, hence, the filament grows in size.

The largest filament grew along the reflecting bound-ary of our simulation and the magnetic field that swept the protons out became a stable magnetic piston. This filament was the largest one because the density of the cloud was largest close to the boundary and because it was aligned with the flow direction of the pair cloud. It had available the largest pool of directed flow energy, which is converted into magnetic energy by the filamen-tation instability.

The filament evolved into a pair jet that was sepa-rated magnetically from the expelled and shocked am-bient plasma. The front of the jet propagated with the speed 0.15c along the boundary and expanded laterally at a speed that amounted to 0.03c. The growth of the fil-ament was limited by our simulation box size and by the limited cloud size; a decrease of the ram pressure would inevitably lead to a weakening of the filamentation insta-bility and to a collapse of the jet. But it appears that, as long as the pair cloud has enough ram pressure, the filaments can grow to arbitrarily large sizes if the fila-mentation instability develops between a pair cloud and an electron-proton plasma at least for plasma parameters similar to those we used here.

We compared the structure of our jet (Fig. 7) with

that of a hydrodynamic jet (Fig. 1(b)) and we found

several similarities. The magnetic piston acts as a dis-continuity between the outer cocoon of the jet and the inner one. Hydrodynamic models attribute this role to a contact discontinuity. Sharp contact discontinuities can only develop in a collisional plasma and hence we could not observe such a discontinuity here. The magnetic pis-ton balanced thermal against magnetic pressure and it was thus similar to a tangential discontinuity.

Hydrodynamic jets in the model proposed by Ref. [6] and the jet in collisionless plasma we observed here have a flat head. We note here that not all jet models predict a flat head [10]. The size of the head was comparable to one wavelength of the waves driven by the filamentation instability. It amounted in our simulation to about 4-5 electron skin depths. If we take into account the fact that our simulation resolved only one half of the head, its width will be of the order of 10 electron skin depths. It is likely that this width will vary with the ratio between the densities of the pair cloud and the ambient plasma and the temperatures of both. Like for hydrodynamic jets the jet in our simulation deflected the ambient protons around its head and into the outer cocoon.

The inner cocoon in a hydrodynamic model is defined as the region, where the jet converts kinetic energy into thermal pressure. The thermal pressure pushes the con-tact discontinuity and the shocked ambient material away from the center of the jet. This was accomplished by our jet in the interval where protons and streaming positrons coexisted and hence we can interpret this region as the inner cocoon. We could not detect a shock between the inner cocoon and the jet plasma. The pair plasma was probably too hot to yield a strong shock and the back-streaming pairs mixed with the pair plasma in the center of the jet instead.

The distributions of the carriers of positive charge re-sembled those of the light and heavy material in leptonic hydrodynamic jets. We found only one difference: a di-lute population of positrons, which was accelerated by the electric field of the jet’s head, escaped from the jet into its upstream region.

Our test simulations demonstrated that the reflecting boundary condition does not lead to artifacts apart from fixing the phase of the filamentation modes. It cut the simulation time in one half and it helped us to determine the various parts of the jet and connect the particle and field distributions to those of the basic filamentation in-stability. We expect that instabilities may twist the jet in the absence of the rigid spine. Future work should thus address the stability of a jet that does not have a rigid jet spine in the form of a reflecting boundary condition. Fu-ture work should also determine if an inner cocoon forms if the pair cloud is colder than the one here.

Acknowledgements: M. E. D. acknowledges

finan-cial support by a visiting fellowship from the ´Ecole

Na-tionale Sup´erieure de Lyon, Universit´e de Lyon. DF and RW acknowledge support from the French Na-tional Program of High Energy (PNHE). GS and MB wish to acknowledge support from EPSRC (grant No: EP/N027175/1). The simulations were performed with the EPOCH code financed by the grant EP/P02212X/1 on resources provided by the French supercomputing fa-cilities GENCI through the grant A0030406960.

[1] M. Diaz Trigo, J. C. A. Miller-Jones, S. Migliari S., J. W. Broderick, and T. Tzioumis, Nature, 504, 260 (2013). [2] H. Falcke, and P. L. Biermann Astron. Astrophys., 308,

321 (1996).

[3] I. F. Mirabel, and L. F. Rodriguez, Ann. Rev. Astron. Astrophys., 37, 409 (1999).

[4] A. M. Stirling, R. E. Spencer, C. J. de la Force, M. A. Garrett, R. P. Fender, R. N. Ogley, Mon. Not. R. Astron. Soc., 327, 1273 (2001).

[5] G. Madejski, and M. Sikora, Ann. Rev. Astron. Astro-phys., 54, 725 (2016).

[6] O. Bromberg, E. Nakar, T. Piran, and R. Sari, Astro-phys. J., 740, 100 (2011).

[7] P. Kumar, and B. Zhang B., Phys. Rep., 561, 1 (2015). [8] T. Leismann, L. Anton, M. A. Aloy, E. Muller, J. M.

Marti, J. A. Miralles, and J. M. Ibanez, Astron. Astro-phys., 436, 503 (2005).

[9] J. M. Marti, E. Muller, J. A. Font, J. M. Z. Ibanez, A. Marquina, Astrophys. J., 479, 151 (1997).

[10] J. M. Marti, and E. Muller, Living Rev. Relat., 6, 7 (2003).

[11] K. M. Ferriere, Rev. Mod. Phys., 73, 1031 (2001). [12] P. Jean, W. Gillard, A. Marcowith, and K. Ferriere,

As-tron. & Astrophys., 508, 1099 (2009).

[13] K. I. Nishikawa, J. T. Frederiksen, A. Nordlund, Y. Mizuno, P. E. Hardee, J. Niemiec, J. L. Gomez, A. Pe’er, I. Dutan, A. Meli, H. Sol, M. Pohl, and D. H. Hartmann, Astrophys. J., 820, 94 (2016).

[14] M. Lyutikov, V. I. Pariev, D. C. Gabudza, Mon. Not. R. Astron. Soc., 360, 869 (2005).

[15] R. Ruffini, G. Vereshchagin, and S. Xue, Phys. Rep., 487, 1 (2004).

[16] R. Moll R., Astron. Astrohys., 512, A5 (2010).

[17] Y. Mizuno, M. Pohl, J. Niemiec, B. Zhang, K. I. Nishikawa, and P. E. Hardee, Astrophys. J., 726, 62 (2011).

[18] E. S. Weibel, Phys. Rev. Lett., 2, 38 (1959).

[19] M. E. Dieckmann, and A. Bret, Mon. Not. R. Astron. Soc., 473, 198 (2018).

[20] F. Califano, F. Pegoraro, and S. V. Bulanov, Phys. Rev. E, 56, 963 (1997).

[21] M. V. Medvedev, and A. Loeb, Astrophys. J., 526, 697 (1999).

[22] M. Lemoine, and G. Pelletier, Mon. Not. R. Astron. Soc.,

402, 321 (2010).

[23] M. Lemoine, G. Pelletier, L. Gremillet, and I. Plotnikov, Mon. Not. R. Astron. Soc., 440, 1365 (2014).

[24] M. Milosavljevic, E. Nakar, and A. Spitkovsky, Astro-phys., J. 637, 765 (2006).

Astron. Soc., 477, 5238 (2018).

[26] L. O. Silva, R. A. Fonseca, J. W. Tonge, J. M. Dawson, W. B. Mori, and M. V. Medvedev, Astrophys. J., 596, L121 (2003).

[27] L. Sironi, and A. Spitkovsky, Astrophys. J., 698, 1523 (2009).

[28] A. Yalinewich, and M. Gedalin M., Phys. Plasmas, 17, 062101 (2010).

[29] A. Vanthieghem, M. Lemoine, and L. Gremillet, Phys. Plasmas, 25, 072115 (2018).

[30] M. E. Dieckmann, A. Alejo, G. Sarri, D. Folini, and R. Walder, Phys. Plasmas, 25, 064502 (2018).

[31] M. E. Dieckmann, A. Alejo, and G. Sarri, Phys. Plasmas,

25, 062122 (2018).

[32] T. D. Arber, K. Bennett, C. S. Brady, A. Lawrence-Douglas, M. G. Ramsay, N. J. Sircombe, P. Gillies, R. G. Evans, H. Schmitz, A. R. Bell, and C. P. Ridgers, Plasma Phys. Controll. Fusion, 57, 113001 (2015). [33] B. Eliasson, and P. K. Shukla, Phys. Rep., 422, 225

(2006).

[34] R. B. Miller, An Introduction to the Physics of Intense

Charged Particle Beams (Springer, Boston MA, 1982).

[35] M. Fiore, L. O. Silva, C. Ren, M. A. Tzoufras, and W. B. Mori, Mon. Not. R. Astron. Soc., 372, 1851 (2006). [36] Hosseini J. S. M. Hosseini, and F. Spanier, Phys. Rev. E,

0 1 0 2 -1000 n c 3 4 0 200 (a) X Y 1000 400 2000 600 _{-1.5 -1 -0.5 0 0.5 1 1.5} (b) X 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Y Shocked Jet Inner Cocoon Outer Cocoon Contact discontinuity Head Ambient medium Reverse shock Forward shock