The baryon loading effect on relativistic astrophysical jet transport in the interstellar medium

Astrophysical jets exist widely in the Universe [1–3]. Among them, relativistic jets are often connected with the observed exciting high-energy astrophysical phenomena, such as gamma-ray burst [4–6], active galactic nucleus (AGNs) [7–9], tidal disruption events [10, 11] and compact binary stars [12, 13]. Recently, the Fermi Gamma-Ray Space Telescope has detected a typical AGN [14], in which the jet bulk Lorentz factor is about 10 ∼ 20. However, despite the growing interest in relativistic astrophysical jets, our understanding of them is still limited. Two main reasons are as follows.

On the one hand, because these high-energy astrophysical phenomena are so far away from Earth, in situ probes are unavailable. As a result, many fundamental parameters of relativistic jets are still under debate, especially the jet composition. For simplification purposes, purely leptonic jet components (electrons and positrons) are often assumed. However, recent studies [15, 16] have identified that both leptons and baryons are presented in the relativistic jets. In fact, the influence of the baryon components on the relativistic jet transport in interstellar medium (ISM) is uncertain yet.

On the other hand, the observed morphology of the relativistic jets during their transport in the ISM has been extensively investigated. Specifically speaking, various dichotomies have been given, such as Fanaroff–Riley (FR) type-I/type-II [17–20] and bottom-wide/center-wide [21–23], etc. It is generally thought that the macroscopic hydrodynamic instabilities, e.g. Rayleigh–Taylor [24] and Kelvin–Helmholtz [25, 26] instabilities (KHI), result in the formation of these different morphology. However, to give a more consistent physical explanation of the jet morphology, the microscopic particle dynamics need to be seriously considered [27], especially in regions such as the X-point of magnetic reconnection [28–31] and the discontinuity surface of collisionless shock [32]. From this point of view, kinetic instabilities, with the spatial scale length at the order of the plasma skin depth, have been identified as playing important roles, e.g. the Buneman instability (BI) [33–37] and the Weibel instability (WBI) [38, 39].

In this paper, the effect of the baryon component, or the so-called baryon loading effect (BLE) [40], on the relativistic astrophysical jet transport dynamics in the ISM has been investigated theoretically and by fully kinetic particle-in-cell (PIC) simulations with a novel particle injection module. From a series of PIC simulations of relativistic jets transport in the ISM with various initial parameters including densities, drifting velocities, compositions, and distributions, we show clearly that with the BLE, the jet can transport in a much longer distance, because the jet electrons draw a significant amount of energy from the jet baryons via the Buneman-induced ES waves and the Weibel-mediated shock waves. Moreover, by increasing the BLE, the jet phase space distribution transforms from a bottom-wide-single-peak structure to a center-wide-multiple-peak one. As a wider spread in the transverse phase space distribution indicates a wider jet radius, the BLE can largely influence the observed jet morphology.

The paper is organized as follows. Section 2 introduces the numerical simulation setup. In section 3, representative simulation results demonstrate the BLE on both the relativistic jet transport dynamics and their morphology. To explain, kinetic instabilities and momentum phase space distributions are carefully analyzed in section 4. In section 5, parametric studies systematically show the robustness of our findings. Section 6 gives the conclusions and discussions.

To self-consistently study the relativistic jet–ISM interaction from the microscopic particle dynamics, a series of two-dimensional (2D) PIC simulations are performed. Similar to the simulation model in previous studies of jet-ambient interaction [41, 42], the relativistic jet is injected from the left into the simulation box, which is filled up with ISM. As figure 1 shows, the simulation box has a length of L_x = 5000λ_e and a width of L_y = 50λ_e, where ${\lambda }_{e}=c/{\omega }_{{pe}}$ is the electron skin depth. Here, c is the light speed, ${\omega }_{{pe}}={({n}_{m}{q}_{e}^{2}/{\varepsilon }_{0}{m}_{e})}^{1/2}$ is the electron plasma frequency (n_m is the ISM electron number density, ${\varepsilon }_{0}$, m_e and q_e are the vacuum permittivity, electron mass and electron charge, respectively). The simulation box is divided into 25000 × 250 cells along x and y, respectively, with the Debye length ${\lambda }_{{\rm{De}}}={({\varepsilon }_{0}{{kT}}_{m0}/{n}_{m}{q}_{e}^{2})}^{1/2}$ resolved, where ${T}_{m0}=1\mathrm{keV}$ is the initial ISM plasma temperature and k is the Boltzmann constant. In order to clearly demonstrate how the jet evolves during its transport, the simulation time is kept rather long (${t}_{{\rm{sim}}}=5000{\omega }_{{pe}}^{-1}$). The ISM is composed of electrons and protons, each specie is represented by 8 macro-particles per cell. We use a reduced proton mass ${m}_{p}/{m}_{e}=64$ due to the limitation of the computational resources, while it manages to ensure that the scale separation between electron and proton dynamics can be modeled in fully kinetic simulations [43]. The ISM species follow the Maxwell distribution ${f}_{s}\propto \exp ({m}_{s}{c}^{2}/{T}_{s})(s\,=\,p,e)$, where the initial electron and proton temperatures are both 1 keV.

Zoom In Zoom Out Reset image size

A novel injection module has been implemented in the code to generate a jet flowing into the simulation box from the left boundary at x = 0. Here the jet is charge-neutralized without the requirement of adding an additional field damping layer at the injection boundary, which is different from injecting an electron beam into plasmas [44]. Also, our injection module is different from the 'reflection method' commonly adopted in collisionless shock studies, which take advantage of the reflect boundary condition [39, 45]. Relativistic jets with various initial parameters including densities, drifting velocities, compositions, and distributions can be loaded from our injection module. Furthermore, in order to suppress the numerical Cherenkov noise [46] caused by the relativistic jet (with drifting velocity close to c), the 3rd-order B-Spline shape function, 5th-order weighting, and 6th-order finite different scheme for solving Maxwell's equations have been adopted. The jet is composed of electrons (${e}^{-}$), positrons (${e}^{+}$) and baryons (loading with protons (${p}^{+}$)), which satisfies ${n}_{{e}^{+}}+{n}_{{p}^{+}}={n}_{{e}^{-}}={n}_{j}$ (n_j represents the jet density). The ratio of the relativistic jet density n_j to that of the ISM n_m is chosen to be $\alpha ={n}_{j}/{n}_{m}=1.1$ here, and the jet keeps injecting into the ISM continuously till the simulation ends. All species of the jet follow a drifted Maxwell–Jüttner distribution [47, 48]

$$\begin{eqnarray*}&&{f}_{s}({\boldsymbol{p}})\propto \exp \left[-\displaystyle \frac{{\gamma }_{d}}{{T}_{s}}(\sqrt{1+{{\boldsymbol{p}}}^{2}}-{V}_{d}{p}_{x})\right](s={p}^{+},{e}^{+},{e}^{-}),\end{eqnarray*}$$

where the temperature ${T}_{s}=1\,\mathrm{keV}$ is taken. ${\boldsymbol{p}}$ is the momentum and V_d is the drift velocity along x-direction. The corresponding Lorentz factor is ${\gamma }_{d}={(1-{({V}_{d}/c)}^{2})}^{-1/2}=15$. The Sobol method [49] is used to load relativistic particles at each step.

In the following sections, three typical simulation cases are demonstrated:

• (i)  
  ${e}^{\pm }$ case, the jet is composed of electrons and positrons only, i.e. the leptonic jets;
• (ii)  
  ${e}^{-}{e}^{+}{p}^{+}$ case, the jet is composed of electrons, positrons and protons with densities ${n}_{{p}^{+}}={n}_{{e}^{+}}=0.5{n}_{{e}^{-}}$, i.e. the lepton/baryon mixture jets;
• (iii)  
  ${e}^{-}{p}^{+}$ case, the jet is composed of electrons and protons only, i.e. the baryonic jets.

Certainly, one of the most concerned issues is how far the relativistic jet can stably transport in the ISM. Figure 2 shows a global picture of the jet transport dynamics for the above cases with/without the BLE. From the time evolution of electron densities in figures 2(a)–(c), when the ratio of the baryon component increases, the jet can transport in much longer distances. Specifically, without the BLE in ${e}^{\pm }$ case, the leptonic jet only transports for a distance of x = 2100λ_e (see also in figure 5(a)); while adding 50% baryons in ${e}^{-}{e}^{+}{p}^{+}$ case, the jet transport distance is extended to about x = 2700λ_e (see also in figure 5(b)). Moreover, for a purely baryonic jet in ${e}^{-}{p}^{+}$ case, the jet transports to a much longer distance of nearly x = 3500λ_e (see also in figure 5(c)). From 2D distributions of jet electron density in figures 2(d)–(f) at $t=500{\omega }_{{pe}}^{-1}$, when the BLE increases, we can see that the density modulations of the jet become more violent. Specifically, at the jet bottom region ($0\lt x\lt 100{\lambda }_{e}$), transverse filamentation clearly develops, especially for the ${e}^{-}{p}^{+}$ case. By contrast, at the jet head region ($400{\lambda }_{e}\lt x\lt 500{\lambda }_{e}$), pronounced longitudinal modulations can been seen.

Zoom In Zoom Out Reset image size

These density modulations are caused by kinetic instabilities. From 2D electromagnetic (EM) field distributions in figures 3(a)–(c) around 0 < x < 400λ_e at $t=500{\omega }_{{pe}}^{-1}$, transverse EM modes occur for all cases with the spatial scales close to the ion skin depth ${\lambda }_{i}=c/{\omega }_{{pi}}$, where the WBI dominates and lead to the jet transverse density modulations [41, 42, 50, 51]. When the BLE plays role, B_z becomes stronger and the WBI dominated region (WDR) significantly expands, corresponding to the extension of jet transport distance. Moreover, from 2D electrostatic (ES) field distributions in figures 3(d)–(f) around $400{\lambda }_{e}\lt x\lt 500{\lambda }_{e}$ at the same time, the ES waves can be seen in all cases with the wave length in electron skin depth scale (${\lambda }_{e}=c/{\omega }_{{pe}}$), where the BI dominates and leads to jet longitudinal density modulations [52–54]. In essence, this longitudinal electric field and density modulation are caused by the ion density gradient at the leading edge of the jet, together with the different mobilities of electrons and ions induced by BLE [55]. To clearly see the BLE on the BI-induced ES wave, the integrals of the ES fields ($\phi =\int {E}_{x}{d}_{x}$) are shown in figures 3(g)–(i) at $380{\lambda }_{e}\lt x\lt 480{\lambda }_{e}$. With the increase of the BLE, the Δϕ is also increased, indicating that when electrons transport through this region, they gain more energy from the baryons (see also in figures 5(d)–(f)). However, in e^± case, the jet electrons and positrons cannot draw energy from each other. Instead, they gradually lose energy to the ISM. Thus, the leptonic jet transport in the ISM is a pure stopping process, resulting in a much limited transport distance. Note that during the relativistic jet–ISM interaction, the BI-induced ES mode is coexisting with the WBI-induced EM mode [39, 56, 57].

Zoom In Zoom Out Reset image size

From another perspective, we also see that the WBI-mediated collisionless shocks form and propagate in the ISM from figures 2(a)–(c). By tracing the slope of the piled-up electron density surface of ${n}_{j}+{n}_{m}$, the shock velocity can be estimated as ${v}_{{\rm{sh}},{ee}}\approx 0.42c$, ${v}_{{\rm{sh}},{eep}}\approx 0.54c$, and ${v}_{{\rm{sh}},{ep}}\approx 0.70c$, respectively. Assuming the upstream electron temperature of all cases are the same (T_e = 1 keV), the corresponding ion sound speed will be ${C}_{s}=\sqrt{{{kT}}_{e}/{m}_{p}}\approx 0.0055c$. Thus, their Mach numbers $M={v}_{{\rm{sh}}}/{C}_{s}$ can be calculated as M_ee ≈ 76.02, M_eep ≈ 97.74 and M_ep ≈ 126.70, which are significantly enhanced with the increase of the BLE. In fact, compared with leptonic ones, jets with baryon components provide more free energy for the WBI when injecting with the same drift speed, because the initial energy of each specie is ${\varepsilon }_{s0}={\gamma }_{d}{m}_{s}{V}_{d}^{2}$ ($s={p}^{+},{e}^{+},{e}^{-}$). As a result, with the increase of the BLE, electrons gain more energy from the baryons, leading to stronger collisionless shocks.

Additionally, to show the BLE on electron acceleration during the relativistic jet–ISM interaction, energy evolution and spectra are given in figure 4. Specifically, figure 4(a) shows the energy exchange of leptons and EM fields in e^± and ${e}^{-}{p}^{+}$ case. As the jet continuously injecting into the ISM, in ${e}^{-}{p}^{+}$ case, more free energy are provided with baryons to the WBI. When the WBI saturates and collisonless shock forms, more energy is converted into magnetic field energy (see green lines) and then drawn by electrons (see red lines), compared with that in e^± case. Similarly, the BI also gets more free energy and is converted into more electric field energy (see blue lines). While in figure 4(b), the energy spectra of jet electrons at $t=500{\omega }_{{pe}}^{-1}$ and $t=5000{\omega }_{{pe}}^{-1}$ are shown. In e^± case, few electron acceleration is realized, where the maximum electron energy (normalized by ${m}_{e}{c}^{2}$ ) increases only from 90 to 150 (dashed arrow) during the transport. However, for ${e}^{-}{p}^{+}$ case, the maximum electron energy increased from about 100 to 450 (horizontal solid arrow). Moreover, the particle number of the high-energy electrons also substantially increases (vertical solid arrows). Note that electrons will keep gaining energy until magnetic field energy reach the sub-equipartition level [58].

Zoom In Zoom Out Reset image size

According to the morphology of relativistic astrophysical jets, they have been sorted into typical dichotomies, as mentioned in section 1. Generally speaking, the FR-I jet tends to decelerate in the ISM and is comparatively weak, while the FR-II jet is more powerful and spatially extended. Besides, the most luminous part of the bottom-wide jet lies at the bottom, while that of the center-wide one lies at the center. Owing to the lack of in situ probes, the exact composition of these relativistic jets remains unknown. By kinetic PIC simulations, our findings can provide several implications for jet compositions. For example, from the transport dynamics perspective in section 3, our results imply that FR-I jets should be dominated by leptons, while FR-II ones by baryons.

Similarly, the jet electron momentum phase space distributions also help interpreting the jet composition, as they largely determine the jet morphology. Although a 3D "global jet" model [26, 51] is required to properly simulate the jet morphology, here in our 2D simulation, a wide spread in the transverse phase space distribution (or, a large p_y) indicates wide jet radius. Specifically speaking, in e^± case, without the BLE, jet electrons lose forward momentum during the transport (figure 5(d)), and the WDR remains at the jet bottom region (0 < x < 2000λ_e), where the electrons are heated and scattered via the WBI (figure 5(g)), forming a bottom-wide structure. However, in ${e}^{-}{e}^{+}{p}^{+}$ case, with a certain BLE, jet electrons gain momentums in both directions from jet baryons during the transport (figures 5(e) and (h)), which moves the WDR to the jet body region (2000λ_e < x < 4000λ_e). Moreover, in ${e}^{-}{p}^{+}$ case, jet electrons gain much more momentums during the transport both longitudinally and transversely (figures 5(f) and (i)), and the WDR extends completely to the jet body, forming a center-wide structure.

Zoom In Zoom Out Reset image size

Note that at the leading edge of the jet with baryons, narrow longitudinal momentum peaks develop due to the BI-induced ES waves (figures 5(e) and (f)). Because those electrons already have a relativistic injection speed, they co-move with the ES wave towards upstream, leading to a long dephasing length and continuous acceleration. As a result, the jet phase space distribution transforms from a bottom-wide-single-peak structure to a center-wide-multiple-peak one, where we can interpret that the former should be dominated by leptons, while the latter by baryons. This also implies that the single-zone model [59] (commonly adopted in the jet emission calculation) may only suitable for leptonic jets. Moreover, apart from the diffusive shock acceleration [60–62], this acceleration at the leading edge of the relativistic jet may serve as another mechanism in generating the ultra-high-energy cosmic ray (UHECR) [58, 63].

Systematic parametric studies of the BLE on relativistic jets transport in the ISM are performed. Firstly, five cases with different baryon densities (${n}_{{p}^{+}}/{n}_{{e}^{-}}=0.00/0.25/0.50/0.75/1.00$) are presented in figure 6. Specifically speaking, with the increase of the baryon amount, the jet transport distance elongates from about 2100λ_e to about 3500λ_e in figure 6(a), and the maximum energy of the jet electrons increase from about 150 to nearly 450 in figure 6(b).

Zoom In Zoom Out Reset image size

Secondly, the mass of the baryon component (${m}_{p}/{m}_{e}$) influences the BLE. With larger ${m}_{p}/{m}_{e}$ (e.g. ${m}_{p}/{m}_{e}=64/100/225/400$), the injected protons have more initial energy (${\varepsilon }_{0}={\gamma }_{d}{m}_{p}{V}_{d}^{2}$) and the difference mobilities in each specie are larger, so that more energy can be drawn by jet electrons from kinetic instabilities, leading to stronger BLE. However, large ${m}_{p}/{m}_{e}$ increases the simulation time needed for the formation of collisionless shock, because it decreases the growth rate of the WBI. Thus, more computational resources are required.

Thirdly, the "hardness" of the jet (i.e. $\alpha ={n}_{j}/{n}_{e}$) also influences the BLE. With smaller α (e.g. α = 0.1/0.5), the WBI growth rate is smaller (compared with the presented cases), and the simulation time needed for WBI saturation becomes longer. While a larger α (e.g. α = 5/10) leads to stronger numerical noise, thus we have to adopt higher simulation spatial resolution. Because the BLE on relativistic jet–ISM interaction is all significant for the above cases (${m}_{p}/{m}_{e}=64\sim 400$, α = 0.1 ∼ 10), only one of them (${m}_{p}/{m}_{e}=64$, α = 1.1, with tolerable computational consumption and numerical noise) is presented here in sections 3 and 4.

In this paper, we have theoretically and numerically investigated the effect of the jet baryon component (or BLE) on the relativistic astrophysical jet transport dynamics in the ISM. Through a series of kinetic PIC simulations with various parameters, we find that, on the one hand, with BLE the jet transports in a much longer distance because the jet electrons/positrons draw a significant amount of energy from the jet baryons; on the other hand, the BLE largely influences the observed jet morphology, where the jet electron/positron phase space distribution transforms from a bottom-wide-single-peak structure to a center-wide-multiple-peak one by increasing the BLE.

These findings provide us several implications on astrophysical studies, especially the jet component, which is still under debate. To begin with, from the transport distance perspective, our findings imply that the FR-I jets should be dominated by leptons, while the FR-II ones by baryons. Moreover, from the jet phase space perspective, our findings imply that the bottom-wide-single-peak jets are likely to be composed of leptons, while the center-wide-multiple-peaks jets may be composed of baryons. Thus, the single-zone model may only be suitable for the leptonic jets. In addition, the continuous electron acceleration process at the jet leading edge may serve as another mechanism in generating UHECR. However, we have to admit that, to investigate the BLE on the global jet structure, not only the working surface (where WBI instability and collisionless shock appear), but also the lateral velocity shear (where KHI and Mushroom instability excite), should be taken into consideration, which requires a 3D global jet model [26, 51].

This work is supported by Science Challenging Project, No. TZ2017005, the NSAF, Grant No. U1630246; the National Key Program of S&T Research and Development, Grant No. 2016YFA0401100; the National Natural Science Foundation of China, Grants No. 11575298, No. 91230205, No. 11575031, and No. 11175026; the National Basic Research 973 Projects No. 2013CBA01500 and No. 2013CB834100, and the National High-Tech 863 Project. BQ acknowledges the support from Thousand Young Talents Program of China. The computational resources are supported by the Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (the second phase). YWP would like to thank ZZQ for valuable discussions about the astrophysical backgrounds.