Mushroom-instability-driven Magnetic Reconnections in Collisionless Relativistic Jets

The relativistic jets launched from the supermassive black holes (Blandford et al., 2019) show or imply the enigmatic sub-structure inside them. The elliptical galaxy M87 is one of the best laboratories to study the relativistic jet because it possesses one of the closest jets, so that its jet structure is best resolved (Hada et al., 2011; Walker et al., 2018). Recently, very long baseline interferometer (VLBI) observation in radio band revealed the triple-ridge structure in M87 jet at a distance of $r{\gtrsim}10^{3}r_{\rm g}$ from the central supermassive black hole (Asada et al., 2016; Mertens et al., 2016; Hada, 2017), where the distance is projected on the sky plane. Here, $r_{\rm g}=GM/c^{2}$ is the gravitational radius of the central black hole with its mass $M$, $G$ is the gravitational constant, and $c$ is the speed of light. The triple-ridge structure can be interpreted as the observational feature of the spine-sheath structure of the jet, which is composed of the fast flowing spine and the slow flowing sheath (Laing & Bridle, 2002). The spine-sheath structure is also proposed to explain the short-term variation of TeV $\gamma$-ray emission in blazars (e.g., Mrk 501 and PKS 2155-304) and may be ubiquitous in the relativistic jet (Ghisellini et al., 2005).

In order to understand the physics and the observed features of the relativistic jets, a number of magnetohydrodynamic (MHD) simulations have been performed (Koide et al., 1999; McKinney, 2006; Mizuno et al., 2007; Tchekhovskoy et al., 2011; Porth & Komissarov, 2015; Gourgouliatos & Komissarov, 2018; Nakamura et al., 2018; Davelaar et al., 2019; Matsumoto & Masada, 2019). The general relativistic radiation MHD (GRRMHD) simulations qualitatively reproduced the sheath brightened structure (Chael et al., 2019). In addition, a possible formation of the plasmoids in the sheath is proposed by general relativistic MHD simulations (Nakamura et al., 2018) and general relativistic resistive MHD simulations (Ripperda et al., 2020). However, GR(R)MHD simulations could not explain the formation of the jet spine up to today¹¹1Alternatively, it has long been known that the natural place for under expanded hypersonic jets (in rocket jets, in particular) to dissipate is at the center line. Particles are energized downstream of Mach disks (shocks) that form along the center-line due to compressive waves emanating from the boundary (Sanders, 1983; Edgington-Mitchell et al., 2014)..

Meanwhile, a limited number of particle-in-cell (PIC) simulations of jets have been carried out as a strong tool to explore the kinetic plasma processes of jets (Nishikawa et al., 2003, 2016; Alves et al., 2018; Davelaar et al., 2020). In the 2010s, an important electron-scale shear-instability called Mushroom instability (MI; Alves et al., 2012) was found in the plasma with relativistic bulk shear flow. Importantly, the MI dominates the electron-scale Kelvin-Helmholtz instability when relative bulk speed is greater than ${\sim}0.3c$ (Alves et al., 2015). Based on kinetic plasma simulations (Alves et al., 2012; Grismayer et al., 2013; Liang et al., 2013; Alves et al., 2014, 2015), the mechanism of MI is the following: At first, the interchanging thermal motion of electrons across the shear surface, which is faster than proton motions, results in the electric currents along the jet bulk-motion near the shear surface. Then, the magnetic field is generated mainly parallel to the shear surface. The generated magnetic field further induces the interchanging electron motion across the velocity-shear surface, and the perturbation grows into a number of electron channel-flows.

In this paper, we study the kinetic plasma dynamics triggered by the MI, namely Mushroom-instability-driven magnetic reconnection (MI-driven MR), by carrying out two-dimensional particle-in-cell (PIC) simulations on the transverse plane of the relativistic jets with velocity-shear surface. Here, MR is the topological rearrangement of magnetic fields with dissipation of magnetic energy and resultant heating and acceleration of plasma particles (Parker, 1979; Horiuchi & Sato, 1999; Biskamp, 2000; Zenitani & Hoshino, 2001; Drake et al., 2006; Yamada et al., 2010; Sironi & Spitkovsky, 2014). We demonstrate that the intermittent MI-driven MR results in the electron concentration to the center (i.e., spine). We propose that this new kinetic plasma dynamics will lead to the formation of the jet spine in the relativistic jets.

We perform fully kinetic PIC simulations of the relativistic jets of electron-proton plasma. The simulations are carried out by using a two-dimensional relativistic PIC simulation code PASTEL (Moritaka et al., 2016) in the Cartesian coordinate system with $x$–$y$ plane, which can explicitly solve the Maxwell equations and the Newton-Lorentz equation in a self-consistent way (e.g., Birdsall & Langdon, 1991; Hockney & Eastwood, 1988; Esirkepov, 2001). We set the $x$–$y$ plane in the transverse direction of the jet, i.e., the propagation direction of the jet is perpendicular to the computational domain. Although the computational domain is in two-dimension, the physical quantities have three dimension in space, i.e., $z$-component of physical quantities of the plasma and electromagnetic fields are also time integrated.

The size of the simulation domain is $(409.6\times 409.6)$, where the unit of the length is the electron skin depth $c/\omega_{\rm pe}$, where $\omega_{\rm pe}$ is the plasma frequency of electrons. This simulation domain is resolved by $(8192{\times}8192)$ grid cells, and the cell size is one-half of the Debye length of the electrons. The periodic boundary condition is imposed at each boundary of the computational domain.

The initial condition of the simulation is as follows. We define the cylindrical radius $r=(x^{2}+y^{2})^{1/2}$ and that of the initial velocity-shear surface $r_{\rm shear}=50c/\omega_{\rm pe}$. We set the high-bulk velocity plasma with the speed $v_{\rm bulk}$ whose corresponding Lorentz factor is $\Gamma_{\rm bulk}$ in the region $r\leq r_{\rm shear}$. We note that, in this work, we assume the bulk speed to be zero in the region $r>r_{\rm shear}$ to focus on kinetic plasma phenomena in a simplified model setup, rather than setting a structured jet with velocity shear inside a static ambient plasma, which may be realized in a more realistic situation. We examine a moderately relativistic bulk velocity model with $v_{\rm bulk}=0.9c$ ($\Gamma_{\rm bulk}\simeq 2.3$) and a highly relativistic bulk velocity model with $v_{\rm bulk}\simeq 0.99c$ (i.e., $\Gamma_{\rm bulk}=10$). The temperatures of the electrons and protons are constant in space, and their thermal velocities equal to $0.1c$ and $0.025c$, respectively. The particles are also uniformly distributed and neither magnetic field nor electric current is assumed. We add no artificial initial perturbation to the plasma. The adequacy of our simple assumption on no initial magnetic field, which is expected to be consistent for the jet plasma far from the central black hole, will be discussed in §4.

Since the injection mechanism of the plasmas into the jet is an open question, we simply assumed the initially uniform density plasma composed of protons and electrons as described above. The protons and electrons might be injected into the jet base as a result of the MRs triggered at the interface of jet and turbulent accretion flow (Tchekhovskoy et al., 2011) or the decay of the neutrons generated via the accelerated proton in the underlying accretion flows (Toma & Takahara, 2012).

We carried out two dimensional PIC simulations to study the kinetic plasma dynamics in the transverse direction of the relativistic jet with velocity-shear surface. Although we start the simulation with no initial magnetic field, the magnetic fields are generated by MI near the velocity-shear surface, and it is discovered that intermittent MRs induced by MI, which we refer to as ”MI-driven MRs”, occur in the transverse plane of the jet. As a result, the number density of high energy electrons increases with time in the jet-center region accompanying with the generation and amplification of magnetic field by MI. The maximum Lorentz factor of the accelerated electrons increases with the initial bulk Lorentz factor inside the velocity-shear surface.

We discuss the consistency of our numerical setup for the dynamics of MI-driven MRs with no initial magnetic field. The necessary condition is that the Larmor radius of the electrons in the jets in which the MI-driven MRs are expected to take place is larger than the length of our simulation domain. In this case, the effects of the initial magnetic field can be reasonably neglected. The Larmor radius of the electrons before the acceleration via the MI-driven MRs (i.e., ${\gamma}_{e}\sim 1$) is expected to be $\gamma_{e}m_{e}c^{2}/eB\sim 10^{7}\gamma_{e}(B/1{\rm mG})^{-1}$ cm, as the magnetic field could be dissipated and its strength would be $\sim 1$ mG (e.g., MAGIC Collaboration et al., 2020) while that near the black hole is indicated to be $\gtrsim 1$ G (Event Horizon Telescope Collaboration, 2021). On the other hand, the simulation box size is $\sim 100c/\omega_{\rm pe}\gtrsim 10^{6}$cm. This is because the electron skin depth in the jet can be roughly estimated to be $c/\omega_{{\rm p}e}=c/\sqrt{4\pi e^{2}n_{\rm e}/m_{\rm e}}\gtrsim 10^{4}\,(n_{\rm e}/10^{4}{\rm cm}^{-3})^{-1/2}\,{\rm cm}$, where the number density of electrons of the distant jet will be less than that near the black hole $n_{\rm e}~{}\lesssim 10^{4}\,{\rm cm}^{-3}$ (see, e.g., Event Horizon Telescope Collaboration, 2019; Kawashima et al., 2021). Our simulation setup with no initial magnetic field will be reasonable, at least to study the dynamics of MI-driven MR itself, because of $r_{{\rm L}e}\gtrsim 100\,(c/\omega_{{\rm p}e})$, while the PIC simulation with initial magnetic field will be performed in forthcoming papers to confirm the results.

In the highly relativistic bulk speed model, the stronger magnetic field surrounding the jet-center region is generated, which results in less concentration of the accelerated electrons in the jet center by suppressing the electron motion across the strong magnetic field via the Lorentz force. However, the strong magnetic field surrounding the high-bulk-Lorentz-factor region can be unstable against the kink instability. Here, we consider the Kruskal–Shafranov stability criterion which proposes that the jet is kink unstable if $q=(2\pi r/L)(B_{z}/B_{\phi})<1$ is satisfied, where $r$, $L$, $B_{z}$ and $B_{\phi}$ are the cylindrical radius of the jet, typical length scale of the jet in the $z$ direction, magnetic field strength in the $z$ direction, and the magnetic field strength in the azimuthal direction of the cylindrical jet, respectively. Here, $\phi$ is the azimuthal angle in the cylindrical coordinate system defined as $\tan\phi=y/x$. Figure 9 shows the distribution of $2\pi B_{z}/B_{\phi}$, which means that the plasma jet with $L>r$, which can be easily realized because of the narrow opening angle of the relativistic jets observed in active galactic nuclei, can be unstable in the region with $\log(2\pi B_{z}/B_{\phi})<0$ (i.e., red colored region). Figure 9 represents that the jetted plasma can be kink unstable for the highly relativistic bulk speed model (bottom panel), while it will be stable for the mildly relativistic bulk speed model (top panel). In the highly relativistic bulk Lorentz factor model, strong MI growth results in the generation of intense magnetic field in the $\phi$-direction, which suppresses the transverse motion of electrons (i.e., suppression of the generation of the magnetic field in $r$-direction). As a consequence, the electric current in $r$-direction is suppressed and the generation of the magnetic field in the $z$-direction is reduced. It can be proposed that the highly relativistic shear plasma flow will generate strong magnetic fields with accompanying the MI-driven MRs, and as a consequence, the large scale MRs with drastic deformation of jet structure triggered by the kink-instability and accelerate the electrons more. After the large scale MRs, it is expected that the MI-driven MRs will restart to transport the electrons towards the jet center because a part of the strong magnetic field surrounding the jet center will disappear just after the large scale MRs induced by the kink instability, although 3D PIC simulations are needed to confirm the expected behavior of the plasma including kink instability.

Shown in Appendix (Figure 11), the MI will dominate the electron-scale Kelvin-Helmholtz instability in our relativistic velocity-shear flow, however, it may not be completely negligible. We should emphasize that the three-dimensional effects can lead to the significant destruction/suppression of the large-scale, coherent plasma structure shown in our two-dimensional simulations, since non-linear behavior of structure formation often depends on the spatial dimension of the system. Three-dimensional PIC simulations to simultaneously solve the MI-driven MRs, electron-scale Kelvin-Helmholtz instability, kink instability, and others are out of the scope of this paper, so that these remain as future works.

It should be mentioned that there is a gap between the length-scales of PIC simulations and realistic jets in the universe. The observed jet-spine seems to have the scale with $10^{2-3}r_{\rm g}$ $\sim$ $10^{16-17}\,(M/6.5\times 10^{9}M_{\odot})$ cm, while the diameter of high-Lorentz-factor region in our simulation is $100c/\omega_{{\rm p}e}$ $\sim$ $10^{6}\,(n_{\rm e}/10^{4}{\rm cm}^{-3})^{-1/2}$cm, i.e., the length-scale of the simulation is much smaller than that of the observed jet-spine. The MI-driven MR may not directly but indirectly form the jet spine via triggering an MHD scale mechanism. Another possible scenario is that the MI-driven MRs continuously occur with moving their reconnection points towards inward until the high energy electrons reaches near the jet-center (see Figure 4). The intermittent MRs would enable the electrons to reach the realistic jet-center because of the acceleration of the electrons and the prompt change of the magnetic field topology during MRs, while the initial and amplified magnetic field may arrest the approach of the electrons towards the jet-center. Alternatively, successive shear-driven acceleration mechanism may work on the electrons after the acceleration via MI-driven MRs, and helps the electrons to penetrate the jet-center region, which is motivated by the recent work on the shear-driven acceleration of electrons after the MRs induced by Kelvin-Helmholtz instability (Sironi et al., 2021). More realistic studies with long term PIC simulations with larger simulation box including a finite initial magnetic field and a finite width of shear-layer, which will give us more insight whether the above scenario solving the scale gap problem is correct, remain as a future work.

Here, in order to demonstrate the possible brightening in the jet-center region as a consequence of the MI-driven MRs, Figure 10 displays the electron density multiplied by the strength of the magnetic field, which qualitatively agrees with the synchrotron emissivity, in the moderately relativistic bulk speed model. One may find that the jet-center region becomes brighter with time. In terms of the timescale, if the straightforward extrapolation of the length scale discussed above is correct and is also applicable to non-linear structure-formation in three-dimensional plasma, the concentration of electrons in the jet spine triggered by MI is not inconsistent with the observation of M87. The typical propagation timescale of the relativistic jet up to the distance with appearance of triple-ridge structure for M87 ($\sim 10^{3}r_{\rm S}$) is $t_{\rm dyn}~{}\gtrsim~{}10^{3}r_{\rm S}/c~{}\sim 10^{8}$s. On the other hand, the growth timescale of MI is much shorter than the jet propagation timescale. In practice, the structure-formation timescale via the MI-driven MR should be estimated by $r_{\rm shear}/v_{\rm trans,e}$, where $r_{\rm shear}$ is the radius of high-bulk-Lorentz-factor region and the transverse electron velocity $v_{\rm trans,e}$ is roughly equal to the initial thermal velocity of electrons. Since the radius of M87 jet at $\sim 10^{3}r_{\rm S}$ from the BH is ${\sim}100r_{\rm S}$ (Asada & Nakamura, 2012; Nakamura & Asada, 2013; Hada, 2017) and the thermal velocity of electrons is expected to be $\gtrsim 0.1c$, the formation timescale of the jet spine is estimated to be $\lesssim$ (several) $\times$ $10^{7}$s.²²2The important point is the ratio of the jet-spine width to the distance from the BH. If this ratio is kept, the same scenario can be applied, even if the jet-spine is found in the region closer to the BH. Hence, with respect to the timescale, our scenario of the jet-spine formation via MI-driven MR is consistent with the observed M87 jet.

If the MI-driven MRs take place in the relativistic jets, the azimuthal component of magnetic fields surrounding the jet spine will be generated via the MI. The linearly polarized radio emission will be detected by VLBI observations. The shear-driven acceleration through the MI-driven MRs may also generate the VHE gamma-ray emissions.