Relativistic Jet Simulations of the Weibel Instability in the Slab Model to Cylindrical Jets with Helical Magnetic Fields

: The particle-in-cell (PIC) method was developed to investigate microscopic phenomena, and with the advances in computing power, newly developed codes have been used for several ﬁelds, such as astrophysical, magnetospheric, and solar plasmas. PIC applications have grown extensively, with large computing powers available on supercomputers such as Pleiades and Blue Waters in the US. For astrophysical plasma research, PIC methods have been utilized for several topics, such as reconnection, pulsar dynamics, non-relativistic shocks, relativistic shocks, and relativistic jets. PIC simulations of relativistic jets have been reviewed with emphasis placed on the physics involved in the simulations. This review summarizes PIC simulations, starting with the Weibel instability in slab models of jets, and then focuses on global jet evolution in helical magnetic ﬁeld geometry. In particular, we address kinetic Kelvin-Helmholtz instabilities and mushroom instabilities.


Introduction
Relativistic jets are collimated outflows of ionized matter powered by black holes.Sites for such jets include the collapse of the core of a massive star forming a neutron star or a black hole, the merger of binary neutron stars, supermassive black holes associated with active galactic nuclei (AGN), gamma-ray bursts (GRBs), and pulsars (e.g., [1]).GRBs and blazars produce the brightest electromagnetic phenomena in the universe (e.g., [2]).Despite extensive observational, theoretical, and simulation studies, the understanding of their formation, their interaction with interstellar mediums, and consequently, their observable properties, such as spectra, variability, and polarization (e.g., [3]), remain quite limited.
Astrophysical jets are ubiquitous and exhibit a wide range of plasma phenomena, such as propagation in the interstellar medium, generation/decay of magnetic fields, magnetic reconnection, and turbulence.In these dynamic environments, particle acceleration may be able to achieve the highest level of energies observed in cosmic rays.Many of the processes that determine the evolution of global relativistic jets are very complex, and they occur on small spatial and short temporal scales associated with plasma kinetic effects.It is especially challenging to integrate microscopic physics into global, large-scale dynamics, which is crucial to understand the full dynamics of the jets.Kinetic plasma simulations are traditionally performed using particle-in-cell (PIC) codes, with the intent of addressing particle acceleration and kinetic magnetic reconnection, which cannot be investigated with fluid models (i.e., relativistic magnetohydrodynamic (RMHD) simulations).In particular, PIC simulations indicate that particle acceleration occurs due to kinetic instabilities, such as electron and ion Weibel instabilities (e.g., ).
In general, these simulations confirm that the Weibel instability is dominant among kinetic instabilities develped in weak or nonmagnetized plasma [27].These instabilities, which develop in relativistic outflows, also lead to multiple shock structures.Dynamically changing current filaments and magnetic fields (e.g., [28]) accelerate electrons (e.g., [12]) and cosmic rays, which affect the pre-shock medium [29].In order to model a shock, a relativistic plasma flow is injected from one end of the computational grid and reflected from a boundary at the opposite end.Such simulations are performed by the following: 1D simulations by [30,31], 2D simulations by (e.g., [14,15,19,25,26,32]), and 3D simulations by [33,34].This method creates two identical counter-streaming beams which collide and interact.This approach also simplifies the numerical method, but leads to the drawbackwhere only one forward-moving shock (FS) is generated.In these settings, the backward (reverse) shock (RS) is indistinguishable from FS.There is another method where a jet is injected into an ambient plasma where FS and RS shock structures are fully modelled.Contact discontinuity (CD) is generated due to deceleration of the jet flow by the ambient plasma.The CD is the location where the electromagnetic field, the velocity of the jet, and the ambient plasmas are similar, but the density changes.FS and RS propagate away from the CD into the jet and ambient plasmas (in the CD frame) [18,[21][22][23].Ardaneh et al. [22] showed that FS, RS, and CD separate the jet and ambient plasma into four regions: (1) the unshocked ambient, (2) shocked ambient, (3) shocked jet, and (4) unshocked jet.In this way, the jet-to-ambient density ratio was selected as the appropriate plasma conditions of AGN and GRB jets.The shock formation processes can be investigated temporally and spatially.A leading and trailing shock system develops with strong electromagnetic fields accompanying the trailing shock.PIC simulations where jets are reflected at the simulation boundary were reviewed, including the generation of high-energy particles, by [35].
In this review, we briefly summarize our previous studies from the slab jet case to the global cylindrical jet case and present new three-dimensional simulation results for an electron-positron jet injected into an electron-positron plasma using a long simulation grid in the jet-propagation direction.We also present the results of a new study of global relativistic jets containing helical magnetic fields.The global simulation results, including velocity shears (this time) using a small simulation system, validate the use of the simulation code for the research project.

PIC Simulations in a Slab Model
It is natural to start to perform PIC simulations in a slab model where jets are injected into the whole simulation system.Since we use the periodic boundary conditions in the transverse direction to the jets, we are simulating a part of the jets without taking into account the boundary between the jets and the ambient plasmas.The instabilities generated between jets and ambient plasmas are described later.

Simulation of the Weibel Instability
The Weibel instability is a plasma instability which occurs in homogeneous or nearly homogeneous plasmas, where an anisotropy in the momentum (velocity) space exists [27].The Weibel instability is often referred to as a filamentation instability [36].
The mechanisms of Weibel instability growth are explained as the following: Suppose a field B = B z cos k y is spontaneously generated by thermal fluctuation.Here, k is a wave number, the x, y, and z are the coordinates, and electrons travel along the x-direction.The Lorentz force (−ev × B) then bends the electron trajectories (travelling along the x-direction) along the y-direction, resulting in congregation of the electrons.The resultant current j = −env e sheets (filaments) create a magnetic field, which enhances the original field and thus grows perturbation [28].The Weibel instability is also common in astrophysical plasmas, such as collisionless shock formation in jets, supernova remnants, and GRBs.

Simulation Settings
The code used in this study is an MPI-based parallel version of the relativistic particle-in-cell (RPIC) code, TRISTAN [5,37,38].The simulations have been performed using a grid with (L x , L y , L z ) = (4005, 131, 131) cells and a total of ∼1 billion particles (12 particles/cell/species for the ambient plasma) in the active grid.The electron skin depth is λ s = c/ω pe = 10.0∆,where c = 1 is the speed of light and ω pe = (e 2 n a / 0 m e ) 1/2 is the electron plasma frequency, and the electron Debye length λ D is half of the cell size, ∆.This computational system length is six times longer than that used in the previous simulations [12,39].The jet-electron number density in the simulation reference frame is 0.676n a , where n a is the ambient electron density, and the jet Lorentz factor is γ jt = 15.The jet-electron/positron thermal velocity is v j,th = 0.014c in the jet reference frame.The electron/positron thermal velocity in the ambient plasma is v a,th = 0.05c.As in our previous work (e.g., [12]), the jet is injected in a plane across the computational grid located at x = 25∆ in order to eliminate artificial effects associated with the boundary at x = x min .Radiating boundary conditions are used on the planes at x = x min and x = x max and periodic boundary conditions on all transverse boundaries [37].The jet makes contact with the ambient plasma at a two-dimensional interface spanning the whole computational domain in the y − z plane.In this way, only a small portion of whole jets is studied; that is, the simulation includes the spatial development of nonlinear saturation and dissipation from the injection point to the jet front composed of the fastest-moving jet particles.Therefore, the boundary between jets and ambient plasma is not taken into account, which will be described later.

Simulation Results
Figure 1a,b shows the average (in the y − z plane) of (a) the jet (red), the ambient (blue), and the total (black) electron density, and (b) the electromagnetic field energy divided by the total jet kinetic energy (E Here, "e" and "p" denote the electron and positron.Positron density profiles are similar to the electron profiles, as both particles have the same mass.However, for the electron-ion jets, the densities of the electrons and the ions are slightly different, giving rise to double layers in the plasma [21][22][23].As a result, ambient particles are dragged by the motion of the jet particles up to x/∆∼500.By t = 3250ω −1 pe , the ambient density has evolved into a two-step plateau behind the jet front, which is similar to the electron-ion jet cases [21][22][23].
The maximum density in this shocked region is about three times the initial ambient density.The jet-particle density remains nearly constant up to near the front of the jet.Careful comparisons reveal the differences between the pair jets and the electron-ion jets [21][22][23].The differences arise due to the double layers generated in the trailing and leading edges in the electron-ion jets.are similar to electron profiles.Ambient particles become swept up after jet electrons pass x/Δ ∼ 500.By t = 3250ω −1 pe , the density has evolved into a two-step plateau behind the jet front.The maximum density in this shocked region is about three times the initial ambient density.The jet-particle density remains nearly constant up to near the jet front.
Current filaments and strong electromagnetic fields accompany growth of the Weibel instability in the trailing shock region.The electromagnetic fields are about four times larger than that seen previously using a much shorter grid system (L x = 640Δ).At t = 3250ω −1 pe , the electromagnetic fields are largest at x/Δ∼1700, and decline by about one order of magnitude beyond x/Δ = 2300 in the shocked region (Nishikawa et al. 2006;Ramirez-Ruiz et al. 2007).pe .The jet front propagates with the initial jet speed ( c).Sharp RMHD-simulation shock surfaces are not created (e.g., Mizuno et al. 2009).A leading shock region (linear density increase) moves with a speed between the fastest moving jet particles c and a predicted contact discontinuity (CD) speed of ∼0.76 c (see Section 4).A CD region consisting of mixed ambient and jet particles moves at a speed between ∼0.76 c and the trailing density jump speed ∼0.56 c.A trailing shock region moves with speed 0.56 c; note the modest density increase just behind the large trailing density jump.
Figure 2 shows the phase-space distribution of jet (red) and ambient (blue) electrons at t = 3250ω −1 pe and confirms our shock-structure interpretation.The electrons injected with γ j v x ∼ 15 become thermalized due to Weibel instabililtyinduced interactions.The swept-up ambient electrons (blue) are heated by interaction with jet electrons.Some ambient electrons are strongly accelerated.
Figure 3 shows the velocity distribution of all jet and ambient electrons in the simulation frame.The small peak indicates electrons injected at γ j = 15.Jet electrons are accelerated to a nonthermal distribution.Ambient electrons are also accelerated to speeds above the jet injection velocity.The velocity distributions of jet and ambient electrons near the jet front (at x/Δ > 2300) are also plotted.The fastest jet electrons, γ > 20, are located near the jet front.On the other hand, the fastest ambient electrons are located farther behind the jet front (at x/Δ < 2300).Thus, strong acceleration of the ambient electrons accompanies the strong fields associated with the Weibel instability.are similar to electron profiles.Ambient particles become swept up after jet electrons pass x/Δ ∼ 500.By t = 3250ω −1 pe , the density has evolved into a two-step plateau behind the jet front.The maximum density in this shocked region is about three times the initial ambient density.The jet-particle density remains nearly constant up to near the jet front.
Current filaments and strong electromagnetic fields accompany growth of the Weibel instability in the trailing shock region.The electromagnetic fields are about four times larger than that seen previously using a much shorter grid system (Lx = 640Δ).At t = 3250ω −1 pe , the electromagnetic fields are largest at x/Δ∼1700, and decline by about one order of magnitude beyond x/Δ = 2300 in the shocked region (Nishikawa et al. 2006  pe .The jet front pro initial jet speed ( c).Sharp RMHD-simulati are not created (e.g., Mizuno et al. 2009).A lea (linear density increase) moves with a speed b moving jet particles c and a predicted con (CD) speed of ∼0.76 c (see Section 4).A CD of mixed ambient and jet particles moves at ∼0.76 c and the trailing density jump speed ∼ shock region moves with speed 0.56 c; note t increase just behind the large trailing density Figure 2 shows the phase-space distribu and ambient (blue) electrons at t = 3250ω our shock-structure interpretation.The electr γj vx ∼ 15 become thermalized due to W induced interactions.The swept-up ambient el heated by interaction with jet electrons.Some are strongly accelerated.
Figure 3 shows the velocity distributio ambient electrons in the simulation frame.indicates electrons injected at γj = 15.accelerated to a nonthermal distribution.Amb also accelerated to speeds above the jet injec velocity distributions of jet and ambient elec front (at x/Δ > 2300) are also plotted.The fa γ > 20, are located near the jet front.On the fastest ambient electrons are located fa jet front (at x/Δ < 2300).Thus, strong ac ambient electrons accompanies the strong field the Weibel instability.

Figure 1.
Averaged values of (a) the jet (red), the ambient (blue), and the total (black) electron density, and (b) electric (red) and magnetic (blue) field energy divided by the jet kinetic energy at t = 3250ω −1 pe .Panel (c) shows the evolution of the total electron density in time intervals of δt = 250ω −1 pe .Diagonal lines indicate the motion of the jet front (blue: ≤c), the predicted contact discontinuity (CD) speed (green: ∼0.76c), and the trailing density jump (red: ∼0.56c).Adapted from Figure 1 in [18].
The growth of the Weibel instability creates current filaments and strong electromagnetic fields in the trailing shock region.Since the nonlinear stage is formed in this simulation, the electromagnetic fields are about four times larger than those seen previously in simulations with a much shorter grid system (L x = 640∆).At the simulation time t = 3250ω −1 pe , the electromagnetic fields have the highest intensity at x/∆∼1700, which then declines by about one order of magnitude beyond x/∆ = 2300 in the shocked region [12,39].
Figure 1c shows the total electron density plotted at time intervals of δt = 250ω −1 pe .The jet front propagates with the initial jet speed (≤c).Since anomalous resistivity exists in PIC simulations, sharp RMHD-simulation shock surfaces are not generated (e.g., [40]).A leading shock region (where the linear density increases) moves with a speed between that of the fastest moving jet particles ≤c and a predicted CD value of ∼0.76c.A CD region consisting of mixed ambient and jet particles moves at a speed which is between ∼0.76c, and the trailing density jump speed ∼0.56c.The modest density increase just behind the large trailing density jump should be taken note of.Similar shock structures and their velocities for the electron-ion jets are discussed in [21][22][23].
It is important to show the differences between the reflection and the injection models.The shock is set up by reflecting a cold "upstream" flow from a conducting wall located at x = 0 (Figure 1).The interaction between the incoming beam (that propagates along −x) and the reflected beam triggers the formation of a shock, which moves away from the wall along +x [33].This setup is equivalent to the head-on collision of two identical plasma shells, which would form a forward and reverse shock and contact discontinuity as an injection scheme.However, the forward and reverse shocks are not distinguished as in the injection scheme.Furthermore, the conducting wall corresponds to the contact discontinuity.The simulation is performed in the "wall" frame, where the "downstream" plasma behind the shock is stationary-and on the contrary, in the injection scheme, FS, RS, and CD are moving in the same direction.
In 3D, periodic boundary conditions are employed both in yand in z-directions.Each computational cell is initialized with four particles (two per species) in 2D and two particles (one per species) in 3D.They have performed limited number of experiments with a larger number of particles per cell (up to eight per species in 2D), though essentially obtaining the same results.
Their 3D structure is shown in Figure 2, for a relativistic electron-positron shock with magnetization σ = B 2 /n e m e γ jt c 2 = 0 (top panel) and σ = 10 −3 (bottom panel).The background magnetic field B 0 is initially set along the z-direction, in the same way as for our 2D simulations.The yz slice of the magnetic energy fraction in Figure 2c shows that for σ = 10 −3 , the magnetic field ahead of the shock is primarily organized in pancakes stretched in the direction orthogonal to the background magnetic field (i.e., along y).This can be easily understood, considering that the Weibel instability is seeded by focusing the counter-streaming particles into two channels of charge and current.In the absence of a background magnetic field, the currents tend to be organized into cylindrical filaments, as demonstrated by Spitkovsky [41] and shown in the yz slice of the top panel in Figure 2. In the presence of an ordered magnetic field along z, the particles will preferentially move along the magnetic field (rather than orthogonal), so that their currents will more likely be focused at certain locations of constant z, into sheets elongated along the xy plane.This explains the structure of the magnetic turbulence ahead of the shock in the bottom panel of Figure 2, common to all the cases of weakly magnetized shocks they have investigated (i.e., 0 < σ ≤ 10 −1 ).completeness we compare our results with the unmagnetized case σ = 0, where the non-thermal tail is still evolving to higher and higher energies.We find that strongly magnetized electron-positron shocks, with σ 10 −2 , are poor particle accelerators, in agreement with the conclusions of SS09.The post-shock spectrum at late times (see the black solid line for σ = 10 −2 ) is fully consistent with a Maxwellian distribution.This result does not depend on the reduced dimensionality of our 2D computational domain.We have performed a largescale 3D simulation of an electron-positron perpendicular shock with σ = 10 −1 , and we confirm that the post-shock particle spectrum (dotted cyan line in Figure 7(a)) does not show any evidence for particle acceleration.This undoubtedly proves that the absence of accelerated particles in the 2D simulations of perpendicular strongly magnetized shocks performed by SS09 is a physical consequence of the lack of sufficient selfgenerated turbulence, 10 rather than an artifact of the reduced 10 More precisely, the fluctuations that get self-excited in σ 10 −2 shocks (cyclotron modes and their harmonics) have a short path length for emission and absorption, so they constantly enforce the local thermal equilibrium, giving Maxwellian energy spectra.For weakly magnetized shocks, with σ 3 × 10 −3 , we find efficient particle acceleration, with a non-thermal tail of slope p 2.4 (dashed black line in Figure 7(a)) that contains ∼1% of particles and ∼10% of flow energy, regardless of the magnetization.The low-energy end of the non-thermal tail does not significantly depend on the magnetization (γinj 5γ0, or equivalently ηinj 5), but the high-energy cutoff at saturation is systematically higher for lower magnetizations.This is confirmed by the inset of Figure 7(a), where we plot the evolution in time of the maximum Lorentz factor γmax. Regardless of the magnetization, γmax initially grows as γmax ∝ (ωpit) 1/2 , with a coefficient of proportionality that does not significantly depend on σ .At later times, the maximum energy departs from this scaling, and it saturates at a Lorentz factor γsat which is larger for smaller magnetizations.
For relatively high magnetizations (black for σ = 10 −2 and purple for σ = 3 × 10 −3 in the inset of Figure 7(a)), the maximum energy initially grows as ∝ (ωpit) 1/2 , then saturates at γsat, and finally drops to a smaller value.The drop at late 7 Figure 3 shows the phase-space distribution of the jet (red) and the ambient (blue) electrons at t = 3250ω −1 pe and confirms the shock-structure interpretation.The electrons injected with γ jt v x ∼15 become thermalized due to the Weibel instability, which is induced by interactions.The swept-up ambient electrons (blue) are heated by interaction with the jet electrons.Some ambient electrons are strongly accelerated.tron profiles.Ambient particles become swept ons pass x/Δ ∼ 500.By t = 3250ω −1 pe , the d into a two-step plateau behind the jet front.nsity in this shocked region is about three bient density.The jet-particle density remains to near the jet front.ts and strong electromagnetic fields accome Weibel instability in the trailing shock remagnetic fields are about four times larger eviously using a much shorter grid system t = 3250ω −1 pe , the electromagnetic fields are 00, and decline by about one order of mag-Δ = 2300 in the shocked region (   pe .The jet front propagates with the initial jet speed ( c).Sharp RMHD-simulation shock surfaces are not created (e.g., Mizuno et al. 2009).A leading shock region (linear density increase) moves with a speed between the fastest moving jet particles c and a predicted contact discontinuity (CD) speed of ∼0.76 c (see Section 4).A CD region consisting of mixed ambient and jet particles moves at a speed between ∼0.76 c and the trailing density jump speed ∼0.56 c.A trailing shock region moves with speed 0.56 c; note the modest density increase just behind the large trailing density jump.

WEIBEL INSTABILITY WITH STRONG MAGNETIC FIELDS L11
Figure 2 shows the phase-space distribution of jet (red) and ambient (blue) electrons at t = 3250ω −1 pe and confirms our shock-structure interpretation.The electrons injected with γ j v x ∼ 15 become thermalized due to Weibel instabililtyinduced interactions.The swept-up ambient electrons (blue) are heated by interaction with jet electrons.Some ambient electrons are strongly accelerated.
Figure 3 shows the velocity distribution of all jet and ambient electrons in the simulation frame.The small peak indicates electrons injected at γ j = 15.Jet electrons are accelerated to a nonthermal distribution.Ambient electrons are also accelerated to speeds above the jet injection velocity.The velocity distributions of jet and ambient electrons near the jet front (at x/Δ > 2300) are also plotted.The fastest jet electrons, γ > 20, are located near the jet front.On the other hand, the fastest ambient electrons are located farther behind the jet front (at x/Δ < 2300).Thus, strong acceleration of the ambient electrons accompanies the strong fields associated with the Weibel instability.pe .About 18,600 electrons of both species are selected randomly.Adapted from Figure 2 in [18].
This simulation shows that the shocks are excited through the injection of a relativistic jet into ambient plasma, leading to two distinct shocks (referred to as the trailing shock and the leading shock) and contact discontinuity.It should be noted that the simulations where jets are reflected on the simulation boundary do not show the structure of a leading shock, contact discontinuity, and a trailing (reverse) shock.
For the electron-ion jet case, the mass ratio is m i /m e = 16 and, therefore, the evolution of density (shock) structures are different to those in the electron-positron jet (m i /m e = 1) [22,23].Furthermore, the double layers generated in the trailing and leading edges further accelerate the electrons up to the ion kinetic energy [23].

Simulation of Jets with Velocity-Shears
The generation of shocks in slab jet models have been studied extensively; however, the velocity shears between the jet and the ambient medium still need to be taken into account, where the outflow interaction with an ambient medium induces velocity shearing.
Alves et al. [54] presented the shear surface instability that occurs in the plane perpendicular to that of the ESKHI.These new unstable modes explain the transverse dynamics and the plasma parameter structures similar to those observed in the PIC simulations performed by [47,49,52,55].They named this effect "mushroom instability" (MI), due to the mushroom-like structures that emerge in the electron density, and the 2D simulation in particular.In 3D simulations, the shape of mushrooms cannot seen clearly; nevertheless, they grow to be good and strong [56].
Multi-dimensional PIC simulations confirm the analytic results and further show the appearance of mushroom-like electron density structures in the nonlinear stage of the instability, similar to those observed in the Rayleigh Taylor instability, despite the great disparity in scales and different underlying physics [54,56].This transverse electron-scale instability may play an important role in relativistic and supersonic sheared flow scenarios, which are stable at the (magneto)hydrodynamic level.This aspect will be discussed later in the case of a cylindrical relativistic jet.Macroscopic (dimensional scale c/ω pe ) fields are shown to be generated by this microscopic shear instability, which are relevant for the generation of a DC electric field and toroidal magnetic field (B φ ), acceleration of particles, and emission, as well as seeding magnetohydrodynamic processes at long time-scales [54,56].

Spine-Sheath (Two-Components) Jet Setup
Next, we consider the simulation of a jet with a spine-sheath (two-component) plasma jet structure, which was studied using the counter-streaming plasma setup implemented in simulations by [47,[49][50][51][52][53].In the setup, a jet spine (core) with velocity γ core propagates in the positive x-direction in the middle of the computational box.The upper and lower quarters of the numerical grid contain a sheath plasma that can be stationary or moving with velocity v sheath in the positive x-direction [48,55].This model is similar to that used in the RMHD simulations [44] containing a cylindrical jet spine (core).
Nishikawa et al. [55] performed 3D PIC simulations of the kKHI and the MI for both e ± and e − − p + plasmas.The processes studied here are inspired from the jets from AGN and GRBs that are expected to have velocity shears between a faster spine (core) and a slower sheath wind (stationary ambient plasmas).In these simulations, large velocity shears were studied with relative Lorentz factors of 1.5, 5, and 15.
Figure 4a shows the structure of the B y component of the magnetic field in the y − z plane (jet flows out of the page) at the midpoint of the simulation box, where x = 500∆.Figure 4b depicts 1D cuts along the z axis showing the magnitude and direction of the magnetic field components at the midpoint of the simulation box, where x = 500∆ and y = 100∆ for the e − − p + case at the simulation time t = 300ω −1 pe , with γ jt = 15 [55].In the e − − p + case, magnetic fields appear relatively uniform at the velocity shear surfaces along the transverse y-direction, just as it had been at the velocity shear surfaces along the parallel x-direction, with almost no transverse fluctuations visible in the magnetic field structure.Small fluctuations in the y-direction over distances on the order of ∼10∆ are visible in the currents, whereas small longitudinal mode fluctuations in the x-direction occur over distances ∼100∆.This behavior indicates that the MI generates DC fields in the transverse direction, a fact that has also been seen in the results of global jet simulations without helical magnetic fields [56].For the e ± case, the magnetic field alternates in both the yand z-directions, and these transverse fluctuations occur over distances of the order of ∼100∆, whereas longitudinal mode fluctuations in the x-direction occur over distances ∼100∆ [55].The 1D cuts show that (i) the B y field component dominates in the e − − p + case, (ii) the B y field component is about an order of magnitude smaller for the e ± case, and (iii) the B z component is significant for the e ± case.The 1D cuts also show that there is a sign reversal of the magnetic field on either side of the maximum, which is relatively small for the e ± case but much more significant for the e − − p + case.More details are revealed by the enlargement of the region contained in the squares, as it is shown in Figure 4c.For the e − − p + case, the generated relatively uniform DC magnetic field is symmetric about the velocity shear surface-e.g., note that B y > 0 immediately around the shear surface, and B y < 0 in the jet and the ambient plasmas at somewhat larger distances from the shear surface.It should be noted that this DC magnetic field is generated by the MI and saturated at this time.The MI is also generated in the global e − − p + jet, where this instability generates toroidal magnetic fields that pinch the jet plasma [56].On the other hand, for the e ± case, the generated AC magnetic field resides largely on the jet side of the velocity shear surface.This phenomenon is also found in the global jet simulation [56] and the outflow simulation [57].
The strong electric and magnetic fields in the velocity shear zone can also provide the right conditions for particle acceleration.Nevertheless, the simulations are too short for definitive statements on the efficacy of the process and the resulting spectra.Also, the organization of the field in compact regions will complicate the interpretation of emission spectra, and a spatially resolved treatment of particle acceleration and transport would be mandatory for a realistic assessment, which is beyond the scope of this review paper.Relativistic electrons, for example, can suffer little synchrotron energy loss outside the thin layer of the strong magnetic field.Thus, synchrotron emissivity can be dominated by the shear layer, and in general, this emissivity can depend on how efficiently electrons can flow in and out of the shear layer and be accelerated in the regions with strong magnetic fields.An immediate consequence for radiation modeling is that the energy-loss time of electrons cannot be calculated with the same mean magnetic field that is used to compute emission spectra, because the former includes the volume-filling factor of the strong-field regions.

PIC Simulations of Cylindrical Jets
Cylindrical geometry is the simplest form that can be used to model the relativistic jets.Therefore, cylindrical jets have been used to study the shear instabilities that occur at the interface between a jet and its ambient plasma, where the plasma is unmagnetized and composed of either e ± or e − − p + .Moreover the jet was implemented in the ambient plasma along the x-direction (periodic along the x-direction).Figure 5 shows isocontour images of the x-component of the current, along with the magnetic field lines that are generated by the kinetic instabilities for both e ± and e − − p + jets.The isocontour images show that in the e − − p + jet case, currents are generated in sheet-like layers and the magnetic fields are wrapped around the jet generated by the dominant MI.On the other hand, in the e ± jet case, many distinct current filaments are generated near the velocity shear, and the individual current filaments are wrapped by the magnetic field.Since the growth rates of kinetic instabilities depend on the species of jets, dominant growing modes are different.The clear difference in the magnetic field structure between these two cases may make it possible to distinguish different jet compositions via differences in circular and linear polarization, which are seen clearly in the global jets injected into ambient plasmas [56].Alves et al. [59] considered magnetic field profiles of the form B(r) = B 0 (r/R c )e 1−r/R c e φ +B z e z , where R c is the cross-sectional radius of the jet spine.They also demonstrated that the toroidal magnetic field profiles decay as r −α (with α ≥ 1) and determined that their overall findings are not sensitive to the structure of the magnetic field far from R c .Near the black hole, the poloidal and the toroidal magnetic field components (B z and B α , respectively) are comparable to one another [60].However, the ratio B z /B α decreases with the increase of the distance from the source, and it can be very small at a distance-relevant to astrophysical jets-of ∼100 pc.The characteristic magnetic field amplitude (henceforth denoted as B 0 ) at such distances, B 0 ∼mG, is quite strong in the sense that the ratio σ of the magnetic energy density to plasma rest-mass energy density may exceed unity.In this review, we would like to emphasise the importance of the macroscopic-like instabilities (as, for example, the kink instability), since strong helical magnetic fields can suppress the kinetic instabilities (such as the Weibel instability, kKHI, and MI) and a kink-like instability is more likely to occur, as it is shown in [61,62].
Recently, global relativistic PIC simulations have been performed where a cylindrical unmagnetized jet is injected into an ambient plasma in order to investigate shock (Weibel instability) and velocity shear instabilities (the kKHI and the MI) simultaneously [56].Previously, these two processes have been investigated separately.For example, kKHI and MI have been investigated for sharp velocity shear slabs and cylindrical geometries extending across the computational grid (e.g., [55,58,59]).

Simulation Setups of Global Jet Simulations
Recently, global simulations have been performed while involving the injection of a cylindrical unmagnetized jet into an ambient plasma in order to simultaneously investigate shock (Weibel instability) and velocity shear instabilities (kKHI and MI) [56].Previously, these two processes have been investigated separately.For example, kKHI and MI have been investigated for sharp velocity shear slab and cylindrical geometries extending across the computational grid (e.g., [55,58,59]).In this section, we present the results of this new study of global relativistic jets containing helical magnetic fields.
Jets generated from black holes and merging neutron stars, which are then injected into the ambient interstellar medium, are thought (in many cases) to carry helical magnetic fields (e.g., [1]).Since many GRMHD simulations of jet formations show that the generated jets carry helical magnetic fields (e.g., [63]), jets in PIC global simulations are injected into an ambient medium implementing helical magnetic fields near the jet orifice, (e.g., [61,64]).One of the key issues is how the helical magnetic fields affect the growth of the kKHI, the MI, and the Weibel instability.The RMHD simulations demonstrated that jets containing helical magnetic fields develop kink instability (e.g., [65][66][67]).Since the PIC simulations are large enough to include kink instability, a kink-like instability was found in the pair and e − − p + jet cases (e.g., [61,64]).

Helical Magnetic Field Structure
In the simulations of [61,62], cylindrical jets containing a helical magnetic field were injected into an ambient plasma (see Figure 6a).The structure of the helical magnetic field was implemented like that in the RMHD simulations performed by Mizuno et al. [68], where a force-free expression of the field at the jet orifice was used; that is, the magnetic field was not generated self-consistently, e.g., from simulations of jet formation by a rotating black hole.For the initial conditions, the force-free helical magnetic field was used as described in Equations ( 1) and (2) of Mizuno et al. (2014) [65].
The following form was used for the poloidal (B x ) and the toroidal (B φ ) components of the magnetic field determined in the laboratory frame: where r is the radial coordinate in cylindrical geometry, B 0 parameterizes the magnetic field, a is the characteristic radius of the magnetic field (the toroidal field component has a maximum value at a, for a constant magnetic pitch), and α is the pitch profile parameter.
The expressions for describing the helical magnetic field used by [61,62] are written in Cartesian coordinates.Since α = 1 Equation ( 1) was reduced to Equation (2), the magnetic field takes the form: The toroidal component of the magnetic field was created by a current +J x (y, z) in the positive x-direction, and it is defined in Cartesian coordinates as: Here, the center of the jet is located at (y jc , z jc ) and r = (y − y jc ) 2 + (z − z jc ) 2 .The chosen helicity is defined through Equation ( 3), which has a left-handed polarity with positive B 0 .At the jet orifice, the helical magnetic field is implemented without the motional electric fields.This corresponds to a toroidal magnetic field generated by jet particles moving along the +x-direction.
The poloidal (B x : black) and the toroidal (B φ : red) components of the helical magnetic field with a constant pitch (α = 1) are shown in Figure 6b.The toroidal magnetic fields become zero at the center of the jet, as shown by red lines in Figure 6b.To date, simulations with a constant pitch (α = 1) and with b = 200 have been performed using r jet = 20, 40, 80, 120∆ [61,62].Here, b is the dumping factor of the magnetic fields outside the jet.
It should be noted that the structure of the jet formation region is more complicated than what is implemented in the PIC simulations at the present time (e.g., [69,70]).Furthermore, so far these global jet simulations have been performed with the simplest kind of jet structure with a top-hat shape (flat-density profile).A more realistic jet structure needs to be implemented in a future simulation study.

Helically Magnetized Global Jet Simulations with Larger Jet Radii
In this section, we explore how the jet evolution is affected by the helical magnetic field using a short system before performing more large-scale simulations.A schematic of the simulation injection setup is shown in Figure 6b [61,62].The initial jet and ambient (electron and ion) plasma number density measured in the simulation frame is n jt = 8 and n am = 12, respectively.This set of plasma parameters is used for obtaining the simulation results presented in [56,61,64].
In their simulations, the electron skin depth λ s = c/ω pe = 10.0∆,where c is the speed of light (c = 1), ω pe = (e 2 n am / 0 m e ) 1/2 is the electron plasma frequency, and the electron Debye length for the ambient electrons is λ D = 0.5∆.The jet-electron thermal velocity is v jt,th,e = 0.014c in the jet reference frame.The electron thermal velocity in the ambient plasma is v am,th,e = 0.03c, and ion thermal velocities are smaller by (m i /m e ) 1/2 .Simulations were performed using an e ± plasma or an e − −p + (with m p /m e = 1836) plasma for the jet Lorentz factor of 15 and with the ambient plasma at rest (v am = 0).In these short system simulations, a numerical grid with (L x , L y , L z ) = (645∆, 131∆, 131∆) (simulation cell size: ∆ = 1) is used, imposing periodic boundary conditions in transverse directions with a jet radius of r jet = 20∆ [61].In this review, all simulation parameters are maintained as described above except for the jet radius and the size of the simulation grid (which is adjusted based on the jet radius) [62].Therefore, the jet radius is increased from the value r jet = 20∆ up to several values: r jet = 40∆, 80∆, and 120∆, which corresponds to a numerical grid with (L x , L y , L z ) = (645∆, 257∆, 257∆), (645∆, 509∆, 509∆), and (645∆, 761∆, 761∆), respectively.The cylindrical jet with jet radius r jet = 40∆, 80∆, and 120∆ is injected into the middle of the yz plane ((y jc , z jc ) = (129∆, 129∆), (252∆, 252∆), (381∆, 381∆)) at x = 100∆.The largest jet radius (r jet = 120∆) is larger than that (r jet = 100∆) in [56], but the simulation length is much shorter (L x = 2005∆).
Other parameters used in their simulations include the initial magnetic field amplitude parameter B 0 = 0.1c, where σ = B 2 /n e m e γ jt c 2 = 2.8 × 10 −3 is used, and a = 0.25 * r jet .The helical field structure inside the jet is defined by Equations ( 1) and (2).For the magnetic fields outside the jet, a damping function exp [−(r − r jet ) 2 /b] (r ≥ r jet ) is imposed on Equations ( 1) and (2) with the tapering parameter b = 200∆.The final profiles of the helical magnetic field components are similar to those obtained in the case where the jet radius is r jet = 20∆, with the only difference being that a = 0.25 • r jet , as it is shown in Figure 6b.
Figure 7 shows the y-component of the magnetic field (B y ) for two values of the jet radius with r jet = 20∆ and 80∆, respectively.In both cases, the initial helical magnetic field (left-handed; clockwise, viewed from the jet front) is enhanced and disrupted due to the plasma instabilities.Thus, even when shorter simulation systems are used, growing instabilities are affected by the helical magnetic fields.The simple recollimation shock generated in the small jet radius is shown in Figure 7a,b.The currents generated by instabilities in the jets determine these complicated patterns of B y , as it is shown in Figure 7. Using a larger jet radius adds more modes of growing instabilities in the jets, which make the jet structure more complicated.In order to investigate the full development of instabilities in jets with helical magnetic fields, longer simulations are required.
To illustrate the production of acceleration of the particles in the jet, the Lorentz factor of the jet electrons was plotted for the two cases of plasma type used (e − − p + and e ± , respectively) when the jet radius is r jet = 120∆, as it is shown in Figure 8.These observed patterns of the Lorentz factor coincide with the changing directions of the local magnetic fields in the y-direction, which are generated by kinetic instabilities like the Weibel instability, the kKHI, and the MI.The directions of the magnetic fields are indicated by the arrows (black spots) in the x − z plane.(The arrows are better seen when the figure is magnified.)The directions of magnetic fields are determined by the generated instabilities.The structures at the edge of the jets are generated by the kKHI.Moreover, the plots of the Lorentz factor in the y − z plane, which are not presented here, show the production of the MI at the circular edge of the jets.For the jet radii larger than r jet = 80∆, the kKHI and the MI are generated at the jet surface, whereas inside the jet, the Weibel instability is generated together with a kink-like instability, particularly in the case of the e − −p + plasma.Answering the question on how the growth of kink-like instabilities depends on the helical magnetic fields requires further investigations using different parameters, including a, which determines the structure of the helical magnetic field in Equations ( 2) and (3).Furthermore, an imprint on the plasma behavior of different values of the pitch parameter α is also necessary for investigation using Equation (1).
Recently, Dieckmann et al. [57] investigated 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 e − −p + plasma with PIC simulation.They found a mechanism that could collimate the pair cloud into a jet.The electrons and positrons of the cloud expanded rapidly due to their high temperature, which decreased the density of the cloud.A filamentation instability developed between the protons at rest and the moving positrons in the interval, where the latter were still dense.It is noted that it is difficult to distinguish the filament instability from the kKHI, which is shown in the simulation where the electron-positron jet was injected into an electron-positron ambient [56].The 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 positrons were present, which confined the magnetic field to a small spatial span.The effect of the filamentation instability and the resulting magnetic field were to push the protons away from the regions with no protons.The instability and the magnetic field followed the pushed protons and, hence, the filament grew in size.The largest filament grew along the reflecting boundary of their simulation, and the magnetic field that pushed the protons out became a stable magnetic piston.This filament is the largest one because the density of the cloud is largest where it is close to the boundary, and because it was aligned with the flow direction of the pair cloud.The large pool of directed flow energy was converted to magnetic field energy by the filament instability.Similar expansion of electron-positron jet plasmas was observed in the global jets without helical magnetic fields [56].
The filament was generated in a pair jet due to the separation by the generated magnetic field from the expelled and shocked ambient plasma.The front of the jet propagated with the speed 0.15c along the boundary and expanded laterally at a speed that amounted up to about 0.03c.The growth of the filament was limited by their simulation box size and by the limited cloud size; a decrease of the ram pressure would inevitably lead to a weakening of the filamentation instability 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 filamentation instability develops between a pair cloud and an electron-proton plasma, at least for plasma parameters similar to those used here.It should be noted that this simulation study shows the importance of kinetic processes of injected cylindrical plasma clouds using PIC simulation.

Reconnection in Jets with Helical Magnetic Fields
Reconnection is ubiquitous in solar and magnetosphere plasmas, and it is an important additional particle acceleration mechanism for AGN and GRB jets (e.g., [71]).Despite the extensive research on reconnection, most of all reconnection simulations have been performed with the Harris sheet [72].where the unperturbed magnetic fields B are anti-parallel (B = − tanh(x)e y ).The release of energy stored in helical magnetic fields and particle acceleration during reconnection have been proposed as a mechanism for producing high-energy emissions and cosmic rays (e.g., [71,73]).It should be noted that the stored magnetic field energy in anti-parallel magnetic fields in the slab model is not consistent with the helical magnetic fields in the relativistic jets; therefore, a realistic argument on particle acceleration due to reconnection requires consideration of the helical magnetic field in the jets.
The importance of reconnection in jets has been proposed previously, but no kinetic simulation of global jets with helical magnetic fields has been performed before, with the exception of our own simulations [61,62].
Figure 10 shows the vectors of magnetic fields in the quadrant of the front part of the jets.Unfortunately, these vectors do not show the changes in direction which may reveal reconnection sites.In order to find the reconnection region, it is necessary to analyze the critical points (CPs).These CPs or magnetic nulls are the points where the magnitude of the magnetic field vector vanishes [74].These points may be characterized by the behaviour of nearby magnetic field curves or surfaces.The set of curves or surfaces that end on CPs is of special interest because it defines the behaviour of the magnetic field in the neighborhood of CP.
The usual magnetic field configuration satisfies the hyperbolic conditions in which the vector field system has a nonzero real part of eigenvalues.The bifurcation (the topological change) represents the magnetic reconnection in the magnetic field.Thus, the particular sets of CPs, curves, and surfaces can be used to define a skeleton that uniquely characterizes the magnetic field [74].In order to investigate the location of reconnection and its evolution, the method described by Cai, Nishikawa, and Lembege (2007) [74] needs to be employed in future work.

Discussion
In this paper, simulations of relativistic jets have been investigated extensively, starting from the study of the Weibel instability in slab mode, and continuing with simulations of instabilities in velocity-shears.Recently, a cylindrical geometry of the jets has been taken into account to be able to model the jet plasma more realistically.
The global jet simulations performed with large jet radii show the importance of a larger jet radius in PIC simulations for examining the macroscopic processes found in RMHD simulations.Due to the mixing of generated instabilities, the phase space plots of the jet electrons show little or no bunching in comparison to that when the jet radius is smaller, r jet = 20∆.Consequently, recollimation-like shocks occur, rather in the center of the jets.Moreover, the recollimation-like shock structure is dependent on the value of the parameter of the helical magnetic field geometry a.To better understand the production of such recollimation-like shocks, further investigations of PIC simulations performed with even larger radii of the jets are needed.
The Weibel instability is ubiquitous in plasma flows, particularly when the plasma is unmagnetized.However, as shown in one of the simulations with global e − −p + jets without helical magnetic fields, the Weibel instability is suppressed and the MI grows dominantly at the linear stage (see Figure 3a in Nishikawa et al. [56]).On the contrary, for e ± jets, the Weibel instability grows with the kKHI and the MI.
So far, the global jet simulations have been performed only for two values of the ion-to-electron -mass ratio, m i /m e = 1 and 1836.The simulation results obtained even when m p /m e = 1836 indicate that a small grid system is not appropriate for studying the kinetic plasma instabilities altogether in a realistic way.At this time, these two cases will provide us clearer differences between two different cases with the maximum mass ratio.In the simulations performed by Nishikawa et al. [61,62], only a weak magnetization factor was used.Simulations with stronger helical magnetic fields were performed by us, and preliminary results show that MI grows stronger with stronger magnetic fields.However, further investigation is necessary with larger systems.
These simulations show that the excitation of kinetic instabilities like the kKHI, the MI, and the Weibel instability with kink-like instability release the energy stored in the helical magnetic field.Consequently, jet and ambient electrons are accelerated and magnetic fields become turbulent.Furthermore, the accelerating electrons emit radiation, and the turbulent magnetic field induces the polarization of the emitted radiation.
MacDonald & Marscher [3] have developed a radiative transfer scheme that allows the Turbulent Extreme Multi-Zone (TEMZ) code to produce simulated images of the time-dependent linearly and circularly polarized intensity at different radio frequencies.Using the PIC simulation output data as input parameters in the TEMZ code, synthetic polarized emission maps were obtained.These maps highlight the linear and circular polarization expected within the above PIC models.This algorithm is currently being refined to account for slow-light interpolation through the global PIC simulations reviewed here.
Simulations of global jets with helical magnetic fields are promising in regard to providing new insights into jet evolution and associated phenomena.However, at the present time, the length of the system is too small, and a much longer system is required in order to investigate a nonlinear stage.Possibly even when using larger systems, such as a numerical grid with (L x , L y , L z ) = (2005∆, 1005∆, 1005∆), the jet radius 100∆ is not large enough to accommodate the microscopic processes, such as the gyro-motion of electrons and ions.
Therefore, these simulation results only provide some qualitative information which supplements those investigated by RMHD simulations.In the present simulations, jets were injected with a top-hat model.However, jets generated from black holes (either in AGN or in merging systems) have an opening angle and structured shapes.The helical magnetic fields used in the PIC simulations are not formed self-consistently as generated from rotating black holes like those performed in GRMHD simulations, and the initial setup with magnetic fields and the associated jet injection scheme need to be refined in future investigations.Furthermore, simulations of relativistic jets with large Lorentz factor particularly require the inclusion of radiation loss (e.g., [75]).
Since the power of supercomputers is growing rapidly, very large simulations of global jets could be performed, which will provide new insights on jet evolution, including reconnection and associated phenomena such as flares and high-energy particle generation.

Conflicts of Interest:
The authors declare no conflict of interest.

Figure 1 .
Figure 1.Averaged values of (a) jet (red), ambient (blue), and total (black) electron density, and (b) electric (red) and magnetic (blue) field energy divided by the jet kinetic energy at t = 3250ω −1 pe .Panel (c) shows the evolution of the total electron density in time intervals of δt = 250ω −1 pe .Diagonal lines indicate motion of the jet front (blue: c), predicted CD speed (green: ∼0.76 c), and trailing density jump (red: ∼0.56 c).

Figure 2 .
Figure2.Phase-space distribution of jet (red) and ambient (blue) electrons at t = 3250ω −1 pe .About 18,600 electrons of both species are selected randomly.

Figure 1 (
Figure 1(c) shows the total electron density plotted at time intervals of δt = 250ω −1pe .The jet front propagates with the initial jet speed ( c).Sharp RMHD-simulation shock surfaces are not created (e.g.,Mizuno et al. 2009).A leading shock region (linear density increase) moves with a speed between the fastest moving jet particles c and a predicted contact discontinuity (CD) speed of ∼0.76 c (see Section 4).A CD region consisting of mixed ambient and jet particles moves at a speed between ∼0.76 c and the trailing density jump speed ∼0.56 c.A trailing shock region moves with speed 0.56 c; note the modest density increase just behind the large trailing density jump.Figure2shows the phase-space distribution of jet (red) and ambient (blue) electrons at t = 3250ω −1 pe and confirms our shock-structure interpretation.The electrons injected with γ j v x ∼ 15 become thermalized due to Weibel instabililtyinduced interactions.The swept-up ambient electrons (blue) are heated by interaction with jet electrons.Some ambient electrons are strongly accelerated.Figure3shows the velocity distribution of all jet and ambient electrons in the simulation frame.The small peak indicates electrons injected at γ j = 15.Jet electrons are accelerated to a nonthermal distribution.Ambient electrons are also accelerated to speeds above the jet injection velocity.The velocity distributions of jet and ambient electrons near the jet front (at x/Δ > 2300) are also plotted.The fastest jet electrons, γ > 20, are located near the jet front.On the other hand, the fastest ambient electrons are located farther behind the jet front (at x/Δ < 2300).Thus, strong acceleration of the ambient electrons accompanies the strong fields associated with the Weibel instability.

No. 1 Figure 1 .
Figure 1.Averaged values of (a) jet (red), ambient (blue), and total (black) electron density, and (b) electric (red) and magnetic (blue) field energy divided by the jet kinetic energy at t = 3250ω −1 pe .Panel (c) shows the evolution of the total electron density in time intervals of δt = 250ω −1 pe .Diagonal lines indicate motion of the jet front (blue: c), predicted CD speed (green: ∼0.76 c), and trailing density jump (red: ∼0.56 c).

Figure 1 (
Figure 1(c) shows the total electron densi intervals of δt = 250ω −1pe .The jet front pro initial jet speed ( c).Sharp RMHD-simulati are not created (e.g.,Mizuno et al. 2009).A lea (linear density increase) moves with a speed b moving jet particles c and a predicted con (CD) speed of ∼0.76 c (see Section 4).A CD of mixed ambient and jet particles moves at ∼0.76 c and the trailing density jump speed ∼ shock region moves with speed 0.56 c; note t increase just behind the large trailing densityFigure2shows the phase-space distribu and ambient (blue) electrons at t = 3250ω our shock-structure interpretation.The electr γj vx ∼ 15 become thermalized due to W induced interactions.The swept-up ambient el heated by interaction with jet electrons.Some are strongly accelerated.Figure3shows the velocity distributio ambient electrons in the simulation frame.indicates electrons injected at γj = 15.accelerated to a nonthermal distribution.Amb also accelerated to speeds above the jet injec velocity distributions of jet and ambient elec front (at x/Δ > 2300) are also plotted.The fa γ > 20, are located near the jet front.On the fastest ambient electrons are located fa jet front (at x/Δ < 2300).Thus, strong ac ambient electrons accompanies the strong field the Weibel instability.

Figure 5 .
Figure 5. Structure of the flow, from the 3D simulation of an electron-positron shock with magnetization σ = 0 (top) or σ = 10 −3 (bottom).The xy slice shows the particle density (with color scale stretched for clarity), whereas the xz and yz slices show the magnetic energy fraction B (with color scale stretched for clarity).(A color version of this figure is available in the online journal.) dimensionality of the simulation box, as argued by Jones et al. (1998).

Figure 2 .
Figure 2. Structure of the flow, from the 3D simulation of an electron-positron shock with magnetization σ = 0 (top) or σ = 10 −3 (bottom).Panels (b) and (d) show the particle density in the xy slice (with color scale stretched for clarity), whereas Panels (a) and (c)show the the magnetic energy fraction B in the xz and yz slices (with color scale stretched for clarity).Adapted from Figure 5 in [33].

Figure 2 .
Figure2.Phase-space distribution of jet (red) and ambient (blue) electrons at t = 3250ω −1 pe .About 18,600 electrons of both species are selected randomly.

Figure 1 (
Figure 1(c) shows the total electron density plotted at time intervals of δt = 250ω −1pe .The jet front propagates with the initial jet speed ( c).Sharp RMHD-simulation shock surfaces are not created (e.g.,Mizuno et al. 2009).A leading shock region (linear density increase) moves with a speed between the fastest moving jet particles c and a predicted contact discontinuity (CD) speed of ∼0.76 c (see Section 4).A CD region consisting of mixed ambient and jet particles moves at a speed between ∼0.76 c and the trailing density jump speed ∼0.56 c.A trailing shock region moves with speed 0.56 c; note the modest density increase just behind the large trailing density jump.Figure2shows the phase-space distribution of jet (red) and ambient (blue) electrons at t = 3250ω −1 pe and confirms our shock-structure interpretation.The electrons injected with γ j v x ∼ 15 become thermalized due to Weibel instabililtyinduced interactions.The swept-up ambient electrons (blue) are heated by interaction with jet electrons.Some ambient electrons are strongly accelerated.Figure3shows the velocity distribution of all jet and ambient electrons in the simulation frame.The small peak indicates electrons injected at γ j = 15.Jet electrons are accelerated to a nonthermal distribution.Ambient electrons are also accelerated to speeds above the jet injection velocity.The velocity distributions of jet and ambient electrons near the jet front (at x/Δ > 2300) are also plotted.The fastest jet electrons, γ > 20, are located near the jet front.On the other hand, the fastest ambient electrons are located farther behind the jet front (at x/Δ < 2300).Thus, strong acceleration of the ambient electrons accompanies the strong fields associated with the Weibel instability.

Figure 3 .
Figure 3.Phase-space distribution of the jet (red) and the ambient (blue) electrons at t = 3250ω −1 pe .About 18,600 electrons of both species are selected randomly.Adapted from Figure2in[18].

TheFigure 6 .
Figure 6.Magnetic field structure transverse to the flow direction for γjt = 15 is shown in the y-z plane (jet flows out of the page) at the center of the simulation box, x = 500Δ for the e − − p + case (upper row) and the e ± case (lower row) at simulation time t = 300 ω −1 pe .The small arrows show the magnetic field direction in the transverse plane (the arrow length is not scaled to the magnetic field strength).1D cuts along the z axis of magnetic field components Bx (black), By (red), and Bz (blue) are plotted at x = 500Δ and y = 100Δ for (b) the e − − p + case and (e) the e ± case.Note that the magnetic field strength scales in panels (a) (±0.367) and (d) (±0.198) are different.An enlargement of the shear surface structure in the y-z plane contained within the squares in the left panels is shown in the panels (c) and (f) to the right.(A color version of this figure is available in the online journal.)they and z directions and these transverse fluctuations occur over distances on the order of ∼10Δ, whereas longitudinal mode fluctuations in the x direction occur over distances ∼100Δ.

Figure 6 .Figure 4 .
Figure 6.Motion of electrons and/or positrons across the shear surface produces the electric currents shown also in Figure 7 by the arrows.Relativistic jet flow is out of the page and in the − +

Figure 5 .
Figure 5. Isocontour plots of the J x magnitude with magnetic filed lines (one fifth of the jet size) for (a) an e − − p + and (b) an ee ± jet at simulation time t = 300ω −1 pe .The 3D displays are clipped both along and perpendicular to the jet in order to view the interior.Adapted from Figure 4 in Nishikawa et al. [58].

Figure 6 .
Figure 6.Panel (a) shows a schematic simulation setup; a global jet setup.The jet is injected at x = 100∆ with the jet radius r jet at the center of the y − z plane (not scaled).Panel (b) shows the helical magnetic fields, B x (black), B φ (red) with B 0 = 0.01 for the pitch profile α = 1.0 with damping functions outside the jet with b = 200.0.The jet boundary is located at r jet = 20∆ [61].So far, simulations were performed with r jet = 20, 40, 80, 120∆ [62].

Figure 8 .
Figure 8. Panels (a,b) show the 2D plot of the Lorentz factor of jet electrons for e − −p + (a) and e ± (b) jet with r jet = 120∆ at time t = 500ω −1 pe .The arrows (black spots) show the magnetic fields in the x − z plane.Adapted from Figure 3 in Nishikawa et al. [62].

Figure 9
Figure 9 shows the isosurface of the Lorentz factor of the jet electrons for a plasma that is composed of (a) e − −p + and (b) e ± .The 3D isosurface of the averaged jet electron Lorentz factor in a quadrant of the jet front (320 ≤ x/∆ ≤ 620, 381 ≤ y, z/∆ ≤ 531) shows where jet electrons are accelerated (in reddish color) locally.The cross-sections and the surfaces of the jets show complicated patterns that are generated by mixed instabilities, where the fine lines represent the magnetic field lines.(a) (b)

Figure 10 .
Figure 10.Panels show 3D vector plots of the magnetic fields for the (a) e − −p + and (b) e ± jets with r jet = 120∆ at time t = 500ω −1 pe .The colors show the strength of the magnetic fields in the quadrant of the front part of the jets.

Author
Contributions: K.-I.N.: Perform simulations, analyze the data and prepare a manuscript; Y.M.: Understanding PIC simulation results based on RMHD simulations; J.L.G.: Contribute for comparing simulation results to observations; I.D.: Perform PIC simulations; A.M.: Perform simulations, critical reading and discussion on this research; J.N.: Contribute modifying the code for this research; O.K.; Modify the code for this simulation; M.P.: Overlook the simulation results; H.S.: Essential suggestions for this research; N.M.: Perform ray-tracing calculation for polarimetric images; D.H.H.: Useful discussions for this research.

Funding:
This work is supported by NSF AST-0908010, AST-0908040, NASA-NNX09AD16G, NNX12AH06G, NNX13AP-21G, and NNX13AP14G grants.The work of J.N. and O.K. has been supported by Narodowe Centrum Nauki through research project DEC-2013/10/E/ST9/00662. Y.M. is supported by the ERC Synergy Grant "BlackHoleCam-Imaging the Event Horizon of Black Holes" (Grant No. 610058).M.P. acknowledges support through grant PO 1508/1-2 of the Deutsche Forschungsgemeinschaft. Simulations were performed using Pleiades and Endeavor facilities at NASA Advanced Supercomputing (NAS), and using Gordon and Comet at The San Diego Supercomputer Center (SDSC), and Stampede at The Texas Advanced Computing Center, which are supported by the NSF.This research was started during the program "Chirps, Mergers and Explosions: The Final Moments of Coalescing Compact Binaries" at the Kavli Institute for Theoretical Physics, which is supported by the National Science Foundation under grant No. PHY05-51164.The first velocity shear results using an electron−positron plasma were obtained during the Summer Aspen workshop "Astrophysical Mechanisms of Particle Acceleration and Escape from the Accelerators" held at the Aspen Center for Physics (1-15 September 2013).