Fast Magnetic Reconnection Structures in Poynting Flux-dominated Jets

Relativistic jets are observed in different classes of astrophysical sources, from mildly- or highly relativistic in microquasars (black hole X-ray binaries) and active galactic nuclei (AGNs), to ultra-relativistic ones in gamma-ray bursts (GRBs). The observation of polarized nonthermal radiation (from radio to gamma-rays) is suggestive of strong magnetization in these jets, particularly near the nuclear region (e.g., Laurent et al. 2011; Doeleman et al. 2012; Martí-Vidal et al. 2015). There are even observational hints of helical magnetic field structures in these regions (as observed, e.g., in M87, Harris et al. 2003).

The most accepted models explaining the origin of jets combine magnetic processes with rotation. In the seminal work of Blandford & Znajek (1977), jets could be powered by the black hole spin transferred to the surrounding magnetic flow. Alternatively, jets could be driven by magneto-centrifugal acceleration in helical magnetic fields arising from the accretion disk (Blandford & Payne 1982). In both cases, the prediction is that jets should be born as magnetized flows. Observations also indicate that at distances large enough from the source, corresponding to scales of several orders of magnitude of the Schwarzschild radius, these jets should already be kinetically dominated (see, e.g., Nakamura & Asada 2013; Zamaninasab et al. 2014; Christie et al. 2019; Giannios & Uzdensky 2019; Zhang & Giannios 2021), so there should be some efficient conversion (or dissipation) of magnetic into kinetic energy. It has been argued that magnetic reconnection could be an important mechanism in allowing such conversion (e.g., Giannios 2010; de Gouveia Dal Pino & Kowal 2015; de Gouveia Dal Pino et al. 2018; Werner et al. 2018; Giannios & Uzdensky 2019, and references therein). Magnetic reconnection is, in fact, now regarded as a fundamental plasma process in a variety of flaring phenomena in the universe, and not only in the solar system in association with solar flares and the Earth magnetotail storms, where reconnection has been directly observed, but also beyond that. Reconnection is now invoked to explain energy dissipation and specially high-energy nonthermal variable emission in systems like pulsar wind nebulae (e.g., Lyubarsky & Kirk 2001; Clausen-Brown & Lyutikov 2012; Cerutti et al. 2014), jets and accretion disks in microquasars, and AGNs (de Gouveia Dal Pino & Lazarian 2005; Giannios et al. 2009; de Gouveia Dal Pino et al. 2010; Giannios 2010, 2011; Nalewajko et al. 2011; McKinney & Uzdensky 2012; Kadowaki et al. 2015; Khiali et al. 2015; Singh et al. 2015, 2016; Sironi et al. 2015; de Gouveia Dal Pino et al. 2018; Kadowaki et al. 2018a, 2018b; Christie et al. 2019; Fowler et al. 2019; Nishikawa et al. 2020), and GRBs (e.g., Drenkhahn & Spruit 2002; Giannios & Spruit 2007; Zhang & Yan 2011; McKinney & Uzdensky 2012).

Magnetic reconnection occurs when two magnetic fluxes of opposite polarity encounter, then partially break and rearrange their configuration at a velocity V_rec, which is a substantial fraction of the local Alfvén speed if the reconnection is fast (e.g., Priest et al. 2003; Zweibel & Yamada 2009; Lazarian et al. 2020). Different processes, such as plasma instabilities, anomalous resistivity, and turbulence, may trigger a fast reconnection. The latter process, in particular, is very efficient and probably the main driving mechanism of fast reconnection in collisional astrophysical flows. Since turbulence causes an efficient field-fluid slippage and stochasticity of the magnetic field lines, bringing initially distant lines into close separations through Richardson diffusion (Eyink et al. 2013; Jafari et al. 2020), this leads to many patches reconnecting simultaneously, making the reconnection rate very fast (and actually independent of the intrinsic plasma microscopic magnetic resistivity; see, e.g., Lazarian & Vishniac 1999; Eyink & Aluie 2006; Kowal et al. 2009; Eyink et al. 2011; Eyink 2015; Takamoto et al. 2015). The rearrangement of the configuration of the magnetic field may convert magnetic into thermal and kinetic energies (e.g., Yamada et al. 2016), and allows for efficient particle acceleration in a Fermi-like stochastic process (de Gouveia Dal Pino & Lazarian 2005; see also the reviews on this process, e.g., in de Gouveia Dal Pino & Kowal 2015; Matthews et al. 2020). Particle acceleration in reconnection regions has been successfully tested in several numerical works employing multidimensional (2D and 3D) magnetohydrodynamic (MHD) simulations with test particles (e.g., Kowal et al. 2011, 2012; Beresnyak & Li 2016; del Valle et al. 2016; de Gouveia Dal Pino et al. 2018, 2020), and particle-in-cell simulations (e.g., Drake et al. 2006; Zenitani & Hoshino 2007, 2008; Lyubarsky & Liverts 2008; Drake et al. 2010; Cerutti et al. 2012; Clausen-Brown & Lyutikov 2012; Sironi & Spitkovsky 2014; Guo et al. 2015, 2016; Li et al. 2015; Lyutikov et al. 2017; Werner et al. 2018, 2019).

As remarked, this mechanism can operate in magnetized astrophysical flows in general, and particularly in relativistic jets, in the regions near the jet base where they are possibly magnetically dominated (see, e.g., Zamaninasab et al. 2014). Exploring reconnection in these objects is the main focus of this work. Magnetic reconnection seems to be especially helpful in solving puzzles related to the very-high-energy (VHE) emission of blazar jets. Blazars are AGNs with highly beamed relativistic jets that point closely to the line of sight (Blandford & Rees 1978; Urry & Padovani 1995). They are the most common extragalactic sources of gamma-rays, both in the gigaelectronvolt (GeV; e.g., Acero et al. 2015) and teraelectronvolt bands (TeV; e.g., Rieger et al. 2013; Tavecchio et al. 2020). Their powerful nonthermal emission, spanning the entire electromagnetic spectrum has a power-law distribution function with variability timescales ranging from days (as in Mrk 421, Kushwaha et al. 2017) to minutes (as, e.g., in PKS 2155-304, Mrk 501, 3C 279, and 3C 54.3; see Aharonian et al. 2007; Albert et al. 2007b; Ackermann et al. 2016; Britto et al. 2016). This multiwavelength emission is usually attributed to relativistic particles accelerated stochastically in recollimation shocks along the jet and in their head (e.g., Mizuno et al. 2015). However, these shocks may be not powerful enough to probe emission in the magnetically dominated regions of these jets (e.g., Sironi et al. 2013; de Gouveia Dal Pino & Kowal 2015; Bell et al. 2018). Striking examples are some of the blazars mentioned above with very short duration TeV flares (of only a few minutes). These imply explosive and extremely compact acceleration/emission regions (<R_S /c) with Lorentz factors much larger than the typical jet bulk values (which are Γ ≃ 5–10) in order to avoid electron–positron pair creation, and thus, entire gamma-ray absorption within the source (e.g., Begelman et al. 2008). A strong candidate mechanism (and possibly the only one) that seems to be able to circumvent this problem and explain this high variability and compactness of the TeV emission is fast magnetic reconnection involving misaligned current sheets inside the jet (see, e.g., Giannios et al. 2009; Giannios 2013; Kushwaha et al. 2017). A similar process has been also proposed for GRBs (e.g., Giannios 2008; Zhang & Yan 2011).

Instabilities occurring in the jet can drive turbulence, and thus, as described above, to substantial randomness and diffusion of the field lines, leading to the production of current sheets where fast turbulent reconnection is triggered (e.g., Spruit et al. 2001; Giannios & Spruit 2006; Singh et al. 2016; Tchekhovskoy & Bromberg 2016; Barniol Duran et al. 2017; de Gouveia Dal Pino et al. 2018; Gill et al. 2018; Nishikawa et al. 2020). Jets with helical magnetic field structure, in particular, are susceptible to current-driven kink (CDK) instability (Begelman 1998; Giannios & Spruit 2006; Mizuno et al. 2009, 2011, 2012, 2014; Alves et al. 2018; Das & Begelman 2019). MHD simulations of the launching and propagation of relativistic jets indicate that this instability can be triggered where the jet recollimates (Bromberg & Tchekhovskoy 2016; Nishikawa et al. 2020). Several concomitant works (see Porth & Komissarov 2015; Singh et al. 2016; Striani et al. 2016; Tchekhovskoy & Bromberg 2016) have also revealed that this instability can operate in the jet spine over limited regions without disrupting the entire jet structure, which is compatible with the observations that these jets can propagate over very large distances and remain stable. Besides, these works have also confirmed that CDK instability converts substantial magnetic energy into kinetic energy and also drives magnetic reconnection. In particular, the three-dimensional special relativistic magnetohydrodynamic (3D-SRMHD) simulations of rotating, Poynting flux-dominated jets with helical fields of Singh et al. (2016) have identified CDK instability-induced turbulence and the formation of current sheets with fast reconnection rates of ∼0.05 V_A.

Despite the extensive study of CDK, none of these previous works have performed a detailed and systematic identification of the location of all the fast reconnecting current sheets driven by the CDK instability over the entire jet domain, nor a precise determination of the reconnection rates, or the magnetic power released in these sites. This determination is very important particularly because these are the potential sites for particle acceleration and nonthermal emission, as stressed above. In this work, we expand upon the work of Singh et al. (2016), and perform a systematic analysis of reconnection sites inside the simulated jet, in a similar way to that employed before in MHD simulations of turbulent environments (Zhdankin et al. 2013), in shearing-box accretion flows Kadowaki et al. (2018a), or turbulent current sheets with a slab geometry (Kowal et al. 2009, 2020). To this aim, we employ here a modified version of the magnetic reconnection search algorithm developed in Kadowaki et al. (2018a) to incorporate relativistic effects. We apply this search algorithm to several snapshots of the simulated jet which allow us to obtain robust averages of the reconnection rates and the magnetic power of multiple reconnection events, as well as the time-variability pattern of them. As an example to show the robustness of the method, we scale our results to the blazar jet Mrk 421 (Kushwaha et al. 2017) and find that the magnetic reconnection power is compatible with the observed gamma-ray power, and the time variability driven by reconnection is consistent with the observed flaring in gamma-rays. We should remark that, in order to obtain a full understanding of how this magnetic energy released by the CDK instability is channeled into energetic nonthermal particles, in a companion paper we have shown results of particle acceleration, injecting thousands of test particles in this 3D-SRMHD jet simulation (Medina-Torrejon et al. 2021, hearafter MGK21). Preliminary results of this study have been presented in de Gouveia Dal Pino et al. (2018), Kadowaki et al. (2018b), and de Gouveia Dal Pino et al. (2020).

The paper is organized as follows: In Section 2, we describe the numerical methodology. In Section 3, we present our results. In Section 4, we show the applications of our results to the Mrk 421 source. Finally, in Section 5 we draw our conclusions and discuss our findings.

The numerical solution of a Poynting flux-dominated jet under the action of the CDK instability has been obtained from the ideal SRMHD equations, which describe the macroscopic behavior of a relativistic magnetized fluid. To this aim, we have used the RAISHIN code based on a 3+1 formalism of the general relativistic conservation laws of continuity, energy momentum, and Maxwell's equations in a curved spacetime (see Mizuno et al. 2006). In the present work, we have adopted an ideal gas equation of state with p = (Γ − 1)ρ e, where p is the thermal pressure, ρ the rest-mass density, ρ e the specific internal energy density, and Γ = 5/3 the adiabatic index. We have also adopted the flat Minkowski spacetime as an appropriate metric to perform the simulations in the special relativistic regime.

We have employed a Harten–Lax–van Leer–Einfeldt approximate Riemann solver to compute the intercell fluxes of the computational grid (see Einfeldt 1988), a flux-interpolated constrained transport method to maintain a divergence-free condition for the magnetic field (flux CT; Tóth 2000), and a third-order Runge–Kutta scheme to advance the equations in time. Furthermore, we have used a second-order monotonized central-difference slope limiter for the reconstruction step, and since the RAISHIN code uses a conservative scheme, an inversion procedure based on the method of Mignone & McKinney (2007) has been used to transform the conserved variables into the primitive ones (for more details, see also Mizuno et al. 2006, 2011). Finally, the equations have been solved in dimensionless code units (c.u.) without a factor of 4π and considering the light speed as the velocity unit (i.e., c = 1). ⁸

2.1. Initial Conditions

Following the work of Mizuno et al. (2012; see also MGK21), we have used as an initial condition a force-free helical magnetic field profile that decreases as a function of radius and with constant pitch (${ \mathcal P }={{RB}}_{z}{B}_{\phi }$). In a cylindrical coordinate system (R, ϕ, z), these equations are given by

$$\begin{eqnarray}&&{B}_{z}=\displaystyle \frac{{B}_{0}}{1+{\left(R/{R}_{0}\right)}^{2}}\,\mathrm{and}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{B}_{\phi }=-\displaystyle \frac{{B}_{0}(R/{R}_{0}){\left[1+{\left({\rm{\Omega }}{R}_{0}\right)}^{2}\right]}^{1/2}}{1+{\left(R/{R}_{0}\right)}^{2}},\end{eqnarray} \tag{ 2 }$$

where B₀ = 0.7 is the magnetic field amplitude, and R₀ = 0.25 the jet core radius (all in code units). The angular velocity of the jet is given by

$$\begin{eqnarray}{\rm{\Omega }}=\left\{\begin{array}{cc}{{\rm{\Omega }}}_{0} & \mathrm{if}\,\,R\leqslant {R}_{0}\\ {{\rm{\Omega }}}_{0}({R}_{0}/R) & \mathrm{if}\,R\gt {R}_{0}\end{array}\right.,\end{eqnarray} \tag{ 3 }$$

where Ω₀ = 2.0 is the angular velocity amplitude.

We have also used a decreasing rest-mass density and pressure profiles with radius given by

$$\begin{eqnarray}&&\rho ={\rho }_{1}\sqrt{\tfrac{{B}^{2}}{{B}_{0}^{2}}}\,\mathrm{and}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}p=\left\{\begin{array}{cc}{p}_{0} & \mathrm{if}\,R\leqslant {R}_{p}\\ {p}_{0}({R}_{p}/R) & \mathrm{if}\,R\gt {R}_{p}\end{array}\right.,\end{eqnarray} \tag{ 5 }$$

where ρ₁ = 0.8 is the density at the jet spine (R = 0), and p = 0.02 is a constant pressure inside a characteristic radius R_p = 0.5, which is used to keep a nondominant thermal pressure force in the initial force-free configuration.

The initial conditions for the vertical and azimuthal components of the velocity field have been defined from the drift velocity v = c E × B /B² and are given by

$$\begin{eqnarray}&&{v}_{z}=-\displaystyle \frac{{B}_{\phi }{B}_{z}}{{B}^{2}}{\rm{\Omega }}R\,\mathrm{and}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{v}_{\phi }=\left(1-\displaystyle \frac{{B}_{\phi }^{2}}{{B}^{2}}\right){\rm{\Omega }}R.\end{eqnarray} \tag{ 7 }$$

Finally, an unstable perturbed profile has been used in the radial velocity component to trigger the CDK instability (see, e.g., Mizuno et al. 2011, 2012), so that

$$\begin{eqnarray}&&{v}_{R}=\displaystyle \frac{\delta v}{N}\exp \left(-\displaystyle \frac{R}{{R}_{0}}\right)\displaystyle \sum _{n=1}^{N}\cos (m\phi )\sin \left(\displaystyle \frac{\pi {nz}}{{L}_{z}}\right),\end{eqnarray} \tag{ 8 }$$

where δ v = 0.01 is the perturbation amplitude, N = 8 is the total number of wavelengths used to excite the kink mode (m = 1) in the system, and L_z = 6 is the size of the computational domain in the z direction.

2.2. Computational Domain and Boundary Conditions

2.3. Magnetic Reconnection Search Algorithm and Reference Frames

In order to search magnetic reconnection events in our Poynting flux-dominated jet simulations in the relativistic regime, we have used a modified version of the algorithm developed by Kadowaki et al. (2018a; see also, Zhdankin et al. 2013). As in the previous work, we select a sample of cells with a current density value ( J = ∇ × B ) five times higher than the average one taken in the whole system ($| {{\boldsymbol{J}}}_{\max }| \gt \epsilon \langle | {\boldsymbol{J}}| \rangle$, for = 5), and choose those cells where $| {{\boldsymbol{J}}}_{\max }|$ is a local maximum within a subarray data cube of size 3 × 3 × 3 cells. Then, we evaluate the eigenvalues and eigenvectors of the current density Hessian matrix around each local maximum cell to obtain a new local coordinate system (e₁, e₂, e₃; see Figure 1) centered in the magnetic reconnection site since it may not be necessarily aligned with one of the axes of the Cartesian coordinate system (for more details, see Kadowaki et al. 2018a; Zhdankin et al. 2013). The edges of the reconnection (or diffusion) regions are defined when the current density decays to half of the maximum value along the e₁, e₂ and e₃ directions ($| {{\boldsymbol{J}}}_{\mathrm{edge}}| =0.5| {{\boldsymbol{J}}}_{\max }|$).

Zoom In Zoom Out Reset image size

These latter steps are applied in the coordinate frame (hereafter, jet frame), where the relativistic jet simulations are performed (see the previous sections). However, since each local maximum cell can move relative to the jet frame with relativistic velocities, we have introduced the reconnection co-moving frame (tilded variables) to quantitatively evaluate the reconnection events. We then transformed all the variables inside each subarray data cube from the jet to the reconnection frame, via generalized Lorentz transformations. In order to remove false-positive events, we also introduced additional criteria for the magnetic and velocity fields (at the reconnection frame), where we selected only those regions with opposite magnetic field and inflow velocity components at the upper and lower edges of the reconnection sites. Furthermore, all the identified regions in the final sample have at least one outflow velocity component at the left or right edges (see more details in Section 3.2).

With this algorithm, we can evaluate several features of the diffusion regions, such as the reconnection structures, the reconnected components of the magnetic field, and the magnetic reconnection rate, showing the advantages of this analysis. In Section 3, we will focus on the study of these features.

Finally, besides the jet and reconnection reference frames mentioned above, we will also have to deal with the observer's frame (primed variables). We assume that the entire computational domain is moving relative to this frame with a bulk Lorentz factor ${{\rm{\Gamma }}}_{j}^{{\prime} }$, and an angle ${\theta }_{j}^{{\prime} }$ to the line of sight (see Figure 1 for a schematic view of the scenario studied in the present work). This is reasonable, since the initial vertical and azimuthal velocity components (Equations (6) and (7), for Ω₀ = 2.0) are mildly relativistic and defined by the drift velocity (similar assumption was used by Mizuno et al. 2012). In this work, the observer's frame will be used to build a synthetic light curve based on the observations of the Blazar Mrk 421, considering reconnection events as the primary source for dissipation and variable emission at high energies (see more details in Section 4).

2.4. Simulation and Algorithm Parameters

As mentioned previously, the simulations were performed using a computational domain of size 6 × 6 × 6 (c.u.), so that a small portion of the magnetized relativistic jet can be reproduced. Besides, we also considered three different resolutions (120³, 240³, and 480³ cells) for convergence test purposes. For the reconnection search algorithm parameters, we identify cells with $| {{\boldsymbol{J}}}_{\max }| \gt \epsilon \langle | {\boldsymbol{J}}| \rangle$ (for = 5), and consider a subarray data cube of size 3 × 3 × 3 cells to select only the local maximum current densities. Finally, the diffusion region's edges were defined by a current density criterion $| {{\boldsymbol{J}}}_{\mathrm{edge}}| \,=0.5| {{\boldsymbol{J}}}_{\max }|$ for our reference model (m240ep0.5). However, in addition to this, we have also considered extended diffusion regions with $0.1| {{\boldsymbol{J}}}_{\max }|$ and $0.05| {{\boldsymbol{J}}}_{\max }|$. The parameters of the simulations are shown in Table 1, and each model name is composed of the resolution (m120, m240, and m480) plus the edge position values (ep0.5, ep0.1, and ep0.05).

Table 1. Parameters

Model	Simulation		Algorithm
	Computational Domain (c.u.)	Resolution (cells)	Subarray Size (cells)		Edge Position
m120ep0.5	[6 , 6 , 6]	[120, 120, 120]	[3, 3, 3]	5.0	$0.5| {{\boldsymbol{J}}}_{\max }|$
m240ep0.5	[6 , 6 , 6]	[240, 240, 240]	[3, 3, 3]	5.0	$0.5| {{\boldsymbol{J}}}_{\max }|$
m240ep0.1	[6 , 6 , 6]	[240, 240, 240]	[3, 3, 3]	5.0	$0.1| {{\boldsymbol{J}}}_{\max }|$
m240ep0.05	[6 , 6 , 6]	[240, 240, 240]	[3, 3, 3]	5.0	$0.05| {{\boldsymbol{J}}}_{\max }|$
m480ep0.5	[6 , 6 , 6]	[480, 480, 480]	[3, 3, 3]	5.0	$0.5| {{\boldsymbol{J}}}_{\max }|$

Download table as:  ASCIITypeset image

Singh et al. (2016) performed 3D-SRMHD simulations of Poynting flux-dominated relativistic jets with an initial helical magnetic field structure suitable for jets near the launching region (see also Mizuno et al. 2009, 2011, 2012, 2014). Considering models with moderate magnetization (σ ∼ 0.1−1), and different density profiles of the jet and the environment, the authors induced small precession perturbations to trigger the CDK instability that in turn excites turbulent fast magnetic reconnection, as stressed in Section 1. Investigating regions of maximum current density in the jet domain, they derived a reconnection rate of ∼0.05 V_A (see also Takamoto et al. 2015). Extending the work of Singh et al. (2016), we employ here the method developed by Kadowaki et al. (2018a), described in Section 2.3, to search and quantify every reconnection events in these relativistic jets.

We start our study considering the reference model m240ep0.5 (see Table 1). The diagrams in Figure 2 show the current density isosurfaces of half of the maximum ∣ J ∣ (orange color), the magnetic field topology (black lines), the density profile at the middle of the box (y−z plane at x = 0), and the magnetic reconnection sites (or diffusion regions) identified by the algorithm (represented by colored circles along the distorted jet spine by the CDK instability). The different colors of the circles correspond to different values of the current density magnitude at t = 0, 30, 40, and 62 (all in code units) in the jet frame.

Zoom In Zoom Out Reset image size

We note that there are some differences between the simulations presented here and those by Singh et al. (2016), including the reduced size of the box and the use of periodic boundaries in the z direction, rather than outflow boundaries, as in Mizuno et al. (2012). This implies that the CDK instability is continuously driven in the system, rather than propagating downstream and eventually leaving the system as seen in Singh et al. (2016). Therefore, as the instability grows and magnetic energy is transformed into kinetic energy, the amplitude of the helical distortions of the jet spine increases up to a complete disruption of the jet, as we see at t = 60 in Figure 2. In MGK21 (see their Figure 2), we analyzed in detail the evolution of the CDK instability and found that its nonlinear growth starts around t ∼ 30 and saturates with the formation of a plateau in the kinetic energy density around t ∼ 40. The instability drives turbulence, which develops completely after this time inducing magnetic reconnection and the formation of several current sheets. Figure 2 shows that even during the growth of the CDK instability, before saturation, at t = 30, there are already sites of reconnection (see also MGK21). The analysis performed in the next sections will show that only part of them are real locci of fast reconnection. Furthermore, in the same snapshot t = 30, despite the high value of ∣ J ∣ for each site, the associated current density isosurfaces show thick structures implying that the magnetic field lines have not accumulated as much. This reflects in a low average reconnection rate as we will see in Section 3.3. In t = 40 and beyond, on the other hand, there are several reconnection sites following the growing wiggling amplitude all along the jet spine, and the associated current density isosurfaces become thinner, characterizing fast reconnection current sheets (Section 3.3). After t = 62, the helical spine structure is completely disrupted, and the maximum value of the mass density reduces to ∼36% of its initial value.

In all diagrams in Figure 2, several reconnection sites close to each other were identified by the algorithm, due to our choice of the size of the cubic subarray of data (of 3 × 3 × 3 cells; see Section 2.3). With this choice, it has been possible to obtain not only the primary site (with the maximum value of ∣ J ∣), but also the secondary ones along the current sheet produced by the encounter of magnetic fields of opposite polarities. Due to this approach, as can be seen in the diagrams, as time evolves, the fragmentation of the large-scale current isosurfaces (and the diffusion regions) that can be associated with turbulence and fast magnetic reconnection events (after t = 40) as will be seen in the subsequent sections. Similar behavior was obtained in resistive MHD simulations of braided coronal loops conducted by Pontin et al. (2011), where the magnetic reconnection leads to a cascade of multiple small-scale events with turbulent-like behavior (see also Kowal et al. 2020), which is in agreement with the predictions of the turbulent magnetic reconnection mechanism (Lazarian & Vishniac 1999; Eyink et al. 2011; Santos-Lima et al. 2021).

3.1. Magnetic Reconnection Structures

We have analyzed the magnetic reconnection structures for a sample of sites identified by the search algorithm in order to verify how reconnection events can behave along the relativistic jet. Figure 3 shows three candidate sites in the jet frame, selected at t = 50. The 2D maps have been produced using the line integral convolution (LIC) method (Cabral & Leedom 1993) to represent, in the same diagram, the magnetic streamlines, and the magnetic field (top diagrams) and current density (bottom diagrams) magnitudes. The diagrams were obtained by a cubic interpolation of the data (for visualization purposes) and correspond to an arbitrary slice in the e₁ − e₂ plane (the local coordinate system at the jet frame, see Section 2.3). The diagrams show clearly the anticorrelation between the magnetic field and the current density, as we should expect for reconnection events. Besides, the bottom diagrams show the high concentration of the magnetic field around the reconnection sheets. The first column shows an elongated "bow-shaped" magnetic island with an X point-like structure right below, whereas the middle diagrams show a complex topology with at least three magnetic islands. The third column shows similar behavior with two islands (already reconnected) separated by an X point-like structure and accumulated magnetic field lines (in reconnection) below them. This study is useful to follow the time evolution of these structures since each example may represent reconnection events at different stages. Moreover, such analysis allows us to identify false-positives events in our sample (i.e., regions identified by the search algorithm that are actually not associated with reconnection events; see Section 2.3). However, for a complete analysis, it would be necessary to check about 200 sites per snapshot separately (in a total of 660 snapshots for model m240ep0.5), but that is beyond the scope of the present study. For future works, machine-learning techniques can be applied to recognize and select automatically such events from these diagrams.

Zoom In Zoom Out Reset image size

Despite the information provided by the LIC diagrams (Figure 3), the 2D projection of the magnetic field components is not enough to reveal the complex 3D topology of a reconnection event. Besides, 3D reconnection does not always imply that the three components of the magnetic field are annihilating at a time as we see several cases in our work, where only one of the components of the magnetic field is annihilated (see also, e.g., Priest et al. 2003; Parnell et al. 2010; Yamada et al. 2010; Pontin et al. 2011, and references therein). The top diagram in Figure 4 shows a zoom-in plot of the magnetic field streamlines (black and magenta lines) around a sample of identified reconnection sites (green and black circles) in the jet frame at t = 50. The colored arrows represent the axes e₁, e₂, and e₃ (yellow, blue, and red) of the local coordinate system, the orange isosurface corresponds to the associated current sheet with half of the maximum value of ∣ J ∣ at the primary reconnection site (green circle), and the black circles correspond to secondary reconnection events. The black lines represent the asymptotic (nonreconnected) magnetic field far away from the diffusion region, whereas the magenta lines represent the twisted and braided magnetic field lines that produce a thin and strong current density isosurface (as in the diagrams in Figure 2, after t = 40) with the maximum value at the green circle's position. As we expected, the e₁-axis is perpendicular to the current sheet, whereas the e₂–e₃ plane is aligned with it. The e₃-axis matches the direction of the local magnetic guide field (and the current sheet), proving the efficiency of the algorithm in separating the dominant nonreconnected magnetic component from those in reconnection.

Zoom In Zoom Out Reset image size

Finally, the 2D LIC maps in the middle and bottom of Figure 4 show the magnetic and current density magnitudes (as in Figure 3) around the primary reconnection site (green circle in the top diagram). Both maps show a magnetic island topology as a result of the projection of the reconnected components in the e₁−e₂ plane. The local magnetic guide field is still present, but hidden by the projection. Therefore, despite the useful information obtained from these 2D maps, the overall 3D scenario and the real nature of the magnetic islands, which are 2D projections of reconnected flux tubes (Kowal et al. 2011, 2012), should be analyzed carefully. This discussion is important since the formation of current sheets along the jet due to reconnection events is expected (see, e.g., Giannios et al. 2009; Christie et al. 2019), but the real magnetic topology can be far more complex than island-like structures.

3.2. Profiles at the Edges of the Diffusion Regions

In the previous section, we have shown the capability of the search algorithm to recognize magnetic reconnection sites and described qualitatively the 2D and 3D topologies of such events along the relativistic jet. In this section and the next ones, we will present a quantitative analysis of the diffusion region in the reconnection reference frame (tilded quantities), as described in Section 2.3.

Figure 5 shows, from the left to the right, the normalized 2D histograms (with 200 bins in each direction, and considering all events during the system's evolution) of the reconnected magnetic field components, the Alfvén speed, and the inflow and outflow velocities (all in units of c) at the upper and lower edges of the diffusion region (see Figure 1 for the reference model m240ep0.5 ⁹ ). The axes correspond to the values of each variable at the upper and lower edges, and the color bars correspond to the normalized counts. The shapes of the first and third histograms are due to selection effects, because we have constrained the sample to characterize the opposite magnetic field and inflow velocity components on each side of the current sheets in the reconnection events (see Section 2.3). In particular, the signs of the inflow velocity at the upper edge must be negative, and positive at the lower edge (as can be seen in the third diagram) due to the direction of the eigenvectors of the local coordinate system obtained from the Hessian matrix. Also, considering this local system of coordinates, in the evaluation of outflow velocity in the fourth histogram, we removed from the sample, sites with positive velocity at the left edge and negative at the right edge. This is because in such cases the plasma is moving back into the diffusion region instead of going out. Nevertheless, we have kept in the histogram the cases with outflow velocities with the same sign both on the left and right edges since at least one of the edges shows an outflow profile. The second diagram from the left shows that the Alfvén speed values are between 0.07c and 0.71c, with more dominant cases in a range of 0.17c–0.38 c. The third and fourth diagrams indicate mildly relativistic velocities with the highest inflow velocity ${\tilde{V}}_{\mathrm{in}}\sim 0.075\,c$, and the highest outflow velocity ${\tilde{V}}_{\mathrm{out}}\sim 0.28\,c$ (for all sites, the outflow velocity is always higher than the inflow velocity, as one might expect). All the diagrams in Figure 5 show high asymmetric events at the diffusion region's edges, which should be taken into account in the evaluation of the reconnection rate (see, e.g., Cassak & Shay 2007), as will be seen in the next section.

Zoom In Zoom Out Reset image size

Finally, we note that, as can be seen in the second diagram in Figure 5, the high values of the Alfvén speed at the reconnection frame justify the relativistic correction for this quantity. However, the diffusion regions velocities are mildly relativistic with respect to the jet frame. Figure 6 shows that the time distribution of the Lorentz factor (Γ) in the diffusion regions is not higher than 1.09. Thus, although we have used the Lorentz transformation to obtain the velocity profiles around the diffusion region in the reconnection reference frame, such correction will not change significantly with respect to the Galilean transformation.

Zoom In Zoom Out Reset image size

3.3. Magnetic Reconnection Rate

We evaluate the reconnection rate as the ratio between the inflow velocity and the Alfvén speed (${\tilde{V}}_{\mathrm{rec}}={\tilde{V}}_{\mathrm{in}}/{\tilde{V}}_{{\rm{A}}}$) at the upper and lower edges of the diffusion region (in the reconnection reference frame). Furthermore, since the reconnection sites in this turbulent system are highly asymmetric events, as we saw in the previous section, a slight correction has been performed where a weighted arithmetic mean is used to obtain a single value for each site, so that

$$\begin{eqnarray}&&\langle {\tilde{V}}_{\mathrm{rec}}\rangle =\displaystyle \frac{{\left[{\tilde{h}}_{r}| {\tilde{V}}_{{e}_{1}}| /{\tilde{V}}_{{\rm{A}}}\right]}_{\mathrm{lower}}+{\left[{\tilde{h}}_{r}| {\tilde{V}}_{{e}_{1}}| /{\tilde{V}}_{{\rm{A}}}\right]}_{\mathrm{upper}}}{{\left[{\tilde{h}}_{r}\right]}_{\mathrm{lower}}+{\left[{\tilde{h}}_{r}\right]}_{\mathrm{upper}}},\end{eqnarray} \tag{ 9 }$$

where ${\tilde{h}}_{r}$ is distance of the (lower or upper) edges to the center of the diffusion region (see the last sketch in Figure 1), and $| {\tilde{V}}_{{e}_{1}}|$ the absolute inflow velocity in the ${\tilde{e}}_{1}$ direction at the reconnection edges.

Figure 7 shows different histograms of $\langle {\tilde{V}}_{\mathrm{rec}}\rangle$ for model m240ep0.5. The left diagrams correspond to the time distributions obtained during the entire evolution of the simulation ¹⁰ (with 200 bins in each direction), and the right diagrams correspond to 1D distributions obtained between t = 50 and 66 (with 140 bins), when the system is already in a quasi-steady turbulent regime (see also MGK21). The top left diagram shows sporadic and slow reconnection events ($\langle {\tilde{V}}_{\mathrm{rec}}\rangle \lt 0.01$) around t = 20. After t ∼ 35, the number of events and the value of $\langle {\tilde{V}}_{\mathrm{rec}}\rangle$ increase as the system reaches the saturation of the exponential growth of the CDK instability (see Section 3.1 in MGK21). Before that, we note that there is a lack of events, particularly between t = 33 and 34. This is because at these times the CDK is still growing and the magnetic field is essentially getting distorted, but with almost no reconnection. More suspicious are those events found around t = 20, which are very slow and possibly are not real reconnection layers. After t ∼ 35, $\langle {\tilde{V}}_{\mathrm{rec}}\rangle$ increases, having the fastest rate of $\langle {\tilde{V}}_{\mathrm{rec}}\rangle \sim 0.23$ at t = 56.

Zoom In Zoom Out Reset image size

The top right histogram in Figure 7 shows that the distribution of $\langle {\tilde{V}}_{\mathrm{rec}}\rangle$ does not resemble a normal distribution, showing a long tail on the side of the fastest rates (similar to the results obtained by Kadowaki et al. 2018a for nonrelativistic accretion disk systems). This skewed feature is characteristic of a log-normal distribution, and to test this hypothesis we have performed fits (black lines in the right histograms) and the results are presented in Table 2. We have also compared the results of the fit with the four statistical moments of the sample (i.e., mean, variance, skewness, and kurtosis moments). In both histograms, we obtained for the skewness and kurtosis nonzero (positive) values, characteristic of a skewed distribution with a peaked shape near the mean, as we expected. Furthermore, we applied a reduced chi-square statistic to evaluate the quality of the fit, obtaining the poor value of ${\chi }_{\mathrm{red}}^{2}\sim 17$ (see Table 2). Based on this result, we constrained the sample in order to obtain the reconnection rate only for the most symmetric profiles (bottom histograms in Figure 7). ¹¹ This constraint reduces to one-fourth the sample size, but also implies in a reduced chi-square to ${\chi }_{\mathrm{red}}^{2}\sim 2.8$ (see Table 2). Despite this difference, we obtained from the fit of the constrained sample an average reconnection rate $\langle {\tilde{V}}_{\mathrm{rec}}{\rangle }_{f}$ of the order of 0.050 ± 0.021, which does not differ from the value of the original sample ($\langle {\tilde{V}}_{\mathrm{rec}}{\rangle }_{f}=0.051\pm 0.026$, see Table 2). Furthermore, this value is compatible with that found in Singh et al. (2016).

Table 2. Statistics of the ${\tilde{V}}_{\mathrm{rec}}$ Distribution for Model m240ep0.5 (t = 50−66)

Statistics	Symmetric and Nonsymmetric Profiles		Symmetric Profiles
	Sample	Fit (log-normal)	Sample	Fit (log-normal)
Mean	0.047	0.05085 ± 0.00060	0.048	0.04969 ± 0.00038
$\sigma (\sqrt{\mathrm{Variance}})$	0.024	0.02612 ± 0.00082	0.022	0.02146 ± 0.00049
Skewness	0.99	1.676 ± 0.048	1.26	1.376 ± 0.029
Kurtosis	1.99	5.38 ± 0.33	4.17	3.55 ± 0.16
${\chi }_{\mathrm{red}}^{2}$	⋯	17	⋯	2.8
Degrees of freedom	⋯	78	⋯	76

Download table as:  ASCIITypeset image

We have also obtained the time evolution of the average magnetic reconnection rate, $\langle {\tilde{V}}_{\mathrm{rec}}{\rangle }_{s}$ (taken from all identified sites in each snapshot), considering different criteria for the position of the edges of the diffusion region (0.5, 0.1, and $0.05| {{\boldsymbol{J}}}_{\max }| ;$ see Section 2.3). This is shown in the top diagram in Figure 8, where each line corresponds to models m240ep0.5 (blue line), m240ep0.1 (red line), and m240ep0.05 (green line; see Table 1). In general, the evolution of $\langle {\tilde{V}}_{\mathrm{rec}}{\rangle }_{s}$ does not change significantly between t = 16 and 34, when the CDK instability grows exponentially (see MGK21, Figure 2). All models show few and slow events with a peak around t = 32, and a gap between t = 33 and 34. After this time, models m240ep0.1 and m240ep0.05 show convergence, but with $\langle {\tilde{V}}_{\mathrm{rec}}{\rangle }_{s}$ values larger than those obtained for the reference model m240ep0.5 (at the limit of 1σ uncertainty, indicated by the colored shades in Figure 8), particularly after the CDK instability has achieved the saturation and quasi-steady-state turbulence is settled in the system, beyond t = 40. The reconnection regions in models m240ep0.1 and m240ep0.05 are larger than in model m240ep0.5 since the asymptotic magnetic field and the velocity have been measured at locations where the current density decays to 0.1 and 0.05 of the maximum, respectively. Despite the convergence of models m240ep0.1 and m240ep0.05, we adopted for our reference model the conservative criterion of $0.5| {{\boldsymbol{J}}}_{\max }|$ (model m240ep0.5), also used in previous works (see, e.g., Kowal et al. 2009; Zhdankin et al. 2013; Kadowaki et al. 2018a). In fact, we do not expect a convergence of $\langle {\tilde{V}}_{\mathrm{rec}}{\rangle }_{s}$ in such comparisons since we can extrapolate the size of a single identified region.

Zoom In Zoom Out Reset image size

The bottom diagram in Figure 8 compares the time evolution of $\langle {\tilde{V}}_{\mathrm{rec}}{\rangle }_{s}$ obtained for the reference jet model with resolution of 240³ cells (model m240ep0.5, blue line), and with two lower and higher resolution models, i.e., 120³ (model m120ep0.5, red line) and 480³ (model m480ep0.5, green line) cells, respectively. The lowest resolution model (m120ep0.5) has the largest differences, with a smooth increase of $\langle {\tilde{V}}_{\mathrm{rec}}{\rangle }_{s}$ (with no gaps) between t = 23 and 36. After this time, model m120ep0.5 achieves a peak around t ∼ 40, i.e., when the CDK instability reaches a plateau, followed by high variability, but with average fast rates over the quasi-steady-state turbulent regime, as in models m240ep0.5 and m480ep0.5. The highest resolution model (m480ep0.5) ¹² shows the smallest variability, and a gap between t = 31 and 32 (earlier than the one obtained in the reference model m240ep0.5). Furthermore, it converges approximately to model m240ep0.5 after t = 40 (within 1σ uncertainty). The latter results are in agreement with the earlier works of Mizuno et al. (2009, 2012, 2014), where the numerical convergence was carefully tested and found to be between the intermediate and high-resolution models, but not for the lowest one. The absence of significant differences between the models in the bottom diagram are also compatible with the fact that fast magnetic reconnection driven by turbulence is independent of the numerical resistivity (as obtained by Kowal et al. 2009, 2012; Kadowaki et al. 2018a). This is also compatible with the turbulence-induced fast reconnection theory of Lazarian & Vishniac (1999) that predicts that the turbulence speeds up the reconnection independently of the ohmic resistivity of the environment.

In order to verify the resolution effects in the estimates of the width and height of the magnetic reconnection sites, we have also evaluated the time evolution of the average values of these quantities. Figure 9 compares these quantities, where we included the standard deviations of the average values for each model (colored shades; as in the diagrams in Figure 8). As we might expect, the intermediate (m240ep0.5) and high-resolution (m480ep0.5) models also show convergence within 1σ uncertainty after t = 40, which indicates that the $\langle {\tilde{V}}_{\mathrm{rec}}{\rangle }_{s}$ behavior shown in the top diagram in Figure 8 for different edge criteria will not change significantly in the case of the higher resolution simulations.

Zoom In Zoom Out Reset image size

Finally, as mentioned before, the log-normal distribution can be associated with the turbulence. Lately, authors have been using the four statistical moments to characterize compressible MHD turbulence from simulations, evaluating its relation with different sonic and Alfvénic Mach numbers, especially for the interstellar medium (see, e.g., Kowal et al. 2007; Burkhart et al. 2009, 2017; Barreto-Mota et al. 2021, and references therein). In these works, the density distribution of molecular clouds in the presence of turbulence naturally converges to a log-normal profile. It is interesting to note that the $\langle {\tilde{V}}_{\mathrm{rec}}\rangle$ distribution also follows this profile (Figure 7), at the same time that the averaged reconnection rate becomes fast (Figure 8), i.e., in the quasi-steady-state turbulent regime (as can be seen in the last diagram in Figure 2). Such behavior also occurs in accretion disks where turbulence induced by magneto-rotational and Parker–Rayleigh–Taylor instabilities drive fast reconnection (Kadowaki et al. 2018a). Likewise, the results obtained in this section indicate that the turbulence is the process driving fast reconnection events (as predicted by Lazarian & Vishniac 1999) in the relativistic jet, under the action of CDK instability. Other systems where turbulence is excited by distinct physical mechanisms such as Kelvin–Helmholtz and Weibel instabilities also experience fast reconnection (e.g., Kowal et al. 2020; Nishikawa et al. 2020).

3.4. Turbulence Power Spectra

In the previous section, we showed that the 1D histograms of $\langle {\tilde{V}}_{\mathrm{rec}}\rangle$ resemble a log-normal distribution when evaluated at times between t ≃ 50 and 66 (Figure 7), which can be associated with turbulence (see, e.g., Burkhart et al. 2009, 2017; Barreto-Mota et al. 2021). In order to further quantify the presence of turbulence, we have performed a Fourier analysis and evaluated the three-dimensional power spectra of the magnetic and kinetic energy densities. These have been taken from the average of spherical shells between k and k+dk (where $k=\sqrt{{k}_{x}^{2}+{k}_{y}^{2}+{k}_{z}^{2}}$) in the Fourier space, and are given by (see, e.g., Simon et al. 2009)

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}| {\boldsymbol{B}}({\boldsymbol{k}}){| }^{2}=\displaystyle \frac{1}{2}\left[| {B}_{x}({\boldsymbol{k}}){| }^{2}+| {B}_{y}({\boldsymbol{k}}){| }^{2}+| {B}_{z}({\boldsymbol{k}}){| }^{2}\right],\mathrm{and}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}| \sqrt{\rho }{\boldsymbol{v}}({\boldsymbol{k}}){| }^{2}=\displaystyle \frac{1}{2}\left[| \sqrt{\rho }{v}_{x}({\boldsymbol{k}}){| }^{2}+| \sqrt{\rho }{v}_{y}({\boldsymbol{k}}){| }^{2}+| \sqrt{\rho }{v}_{z}({\boldsymbol{k}}){| }^{2}\right].\end{eqnarray} \tag{ 11 }$$

Figure 10 shows the power spectra for the magnetic (top diagram) and kinetic (bottom diagram) energy densities for different times (from t = 10–60; colored continuous lines). A 3D-Kolmogorov spectrum (∝ k^−11/3; red dashed line) was included for comparison. The diagrams show inertial ranges with slopes (k^ν ) between −3.7 and −3.6 for $| \sqrt{\rho }{\boldsymbol{v}}({\boldsymbol{k}}){| }^{2}$ (from ${\mathrm{log}}_{10}k=0.7$ to ${\mathrm{log}}_{10}k=1.4$), and −5.1 < ν < −3.3 for ∣ B ( k )∣² ($1.1\lt {\mathrm{log}}_{10}k\lt 1.4$), both in agreement with a Kolmogorov-like spectrum after t ≃ 50 (green and blue lines), indicating a turbulent energy cascade. The magnetic energy spectrum shows steeper slopes, probably due to the strong guiding magnetic field, which is still present at later times (see, e.g., Kowal et al. 2007). We also note that at the time that the CDK instability achieves the saturation (t = 40), a nearly Kolmogorov-like shape is already present in the inertial range of the spectrum, though it is not yet as well defined as in later times, thus indicating that the turbulence is still not fully developed at this time. These results strengthen those shown in Section 3.3 and Figure 7 at the same time interval. Furthermore, We also notice that $\langle {\tilde{V}}_{\mathrm{rec}}{\rangle }_{s}$ achieves a quasi-steady-state behavior with the fastest rates between t = 50 and 66 (see Figure 8), demonstrating the correlation between turbulence and fast magnetic reconnection events.

Zoom In Zoom Out Reset image size

In this section, as an example, we consider an application of our general results obtained for magnetized relativistic jets to the blazar source Mrk 421. Based on the study of Kushwaha et al. (2017), who analyzed the light curve of this source and its variability features in the high-energy band between 0.1 and 300 GeV, we try to reproduce it here, assuming that this emission pattern is associated to fast magnetic reconnection events.

As we have stressed in Section 1, regions of fast reconnection are sites of efficient stochastic particle acceleration (e.g., de Gouveia Dal Pino & Lazarian 2005; Kowal et al. 2012; del Valle et al. 2016). Our simulations do not allow us to evaluate the exact fraction of released magnetic energy by reconnection that goes into particle acceleration. Nevertheless, in the companion paper MGK21, we have injected test particles within the same SRMHD background jet model studied here, and found that particles with an initial energy 10⁻⁴ m_p c² are efficiently accelerated in the turbulent fast reconnection current sheets of the system up to VHEs ∼10⁷–10⁹ m_p c², after several hundred hours (for background magnetic fields with ∼0.1–10 G), where m_p is the proton mass. These results indicate that the fast magnetic reconnection events not only accelerate the relativistic particles, but also could account for the observed VHE variable emission that these particles may produce in the magnetically dominated regions of blazars.

With this in mind, we here make the assumption that the variability pattern of the reconnection events that we computed can be directly correlated with the variable emission in the high-energy band of sources like the one above, which is expressed in its light curve.

We start by considering that in the dissipation zone of the blazar, roughly half of the magnetic reconnection power is converted into kinetic power accelerating efficiently the particles to VHEs (as found in MGK21), which will then be able to produce high-energy photons. This assumption is compatible with the fact that about 50% of the energy released by reconnection goes into particle acceleration (e.g., Yamada et al. 2016) and thus can be used to constrain the radiative power as well (see also, e.g., Christie et al. 2019).

We should also emphasize that it is out of the scope of this work to discuss the specific radiative nonthermal processes that may lead to photon production by the accelerated particles. This would require a full study involving radiative transfer in the jet background and considerations on the leptonic and/or hadronic composition of the source. Our purpose here is only to benchmark the robustness of our reconnection search method and verify to what extent the reconnection events can be connected with the variability and emission patterns of these relativistic sources. In a forthcoming work, we plan to incorporate radiative transfer effects in our SRMHD code in order to perform a complete reconstruction of observed spectral energy distributions of blazar jets, in a similar way as performed, for instance, in Rodriguez-Ramirez et al. (2019).

Basically, we evaluate here the magnetic power ${\tilde{L}}_{B}({t}_{s})$ released from the reconnection sites (in the reconnection frame) in a snapshot t_s , as follows:

$$\begin{eqnarray}&&{\tilde{L}}_{B}({t}_{s})=\displaystyle \sum _{i=1}^{N({t}_{s})}{\left\langle \displaystyle \frac{{\tilde{B}}_{\mathrm{in}}^{2}}{2}{\tilde{V}}_{\mathrm{in}}{\tilde{w}}_{r}^{2}\right\rangle }_{i},\end{eqnarray} \tag{ 12 }$$

where ${\tilde{B}}_{\mathrm{in}}$ is the reconnected magnetic component, ${\tilde{V}}_{\mathrm{in}}$ is the inflow velocity, and ${\tilde{w}}_{r}^{2}$ is the area of the reconnection site at the edges of the diffusion region, measured in the reconnection frame (see the last sketch in Figure 1). We assume that the emission comes from the entire simulated domain (see Figure 1). To obtain ${\tilde{L}}_{B}({t}_{s})$, we summed the contributions of all N(t_s ) sites identified by the algorithm in the snapshot t_s . Since the velocities of the reconnection sites are mildly relativistic at the jet frame, with Lorentz factors of the order of unity (Γ_r ∼ 1; see Figure 6), then the neglect of a length contraction from Lorentz transformation will not change significantly the results. Therefore, we can estimate the area by the square of the width of the reconnection site measured in the jet frame (${\tilde{w}}_{r}\sim {w}_{r}$). Likewise, we can also evaluate the jet luminosity in the dissipation zone as ${L}_{j}({t}_{s})\sim \eta {\tilde{L}}_{B}({t}_{s})$, where η = 0.5 is the efficiency of conversion of the magnetic power into jet luminosity, as discussed above.

We consider that the simulated jet is moving relative to an observer's reference frame (primed quantities) with a bulk Lorentz factor ${{\rm{\Gamma }}}_{j}^{{\prime} }$, and an angle ${\theta }_{j}^{{\prime} }$ to the line of sight (see the first sketch in Figure 1). This assumption is reasonable as long as the reconnection sites are mildly relativistic at the jet frame. Therefore, we have also considered the observed luminosity ${L}_{j}^{{\prime} }({t}_{s})={\delta }^{{\prime} 4}{L}_{j}({t}_{s})$ related to the co-moving luminosity L_j (t_s ) (assumed to be isotropic at the jet frame), where ${\delta }^{{\prime} }\,={\left[{{\rm{\Gamma }}}_{j}^{{\prime} }(1-{\beta }_{j}^{{\prime} }\cos {\theta }_{j}^{{\prime} })\right]}^{-1}$ is the relativistic Doppler factor, and ${\beta }_{j}^{{\prime} }$ is the jet bulk velocity. We have also evaluated the time intervals measured in the observer frame as ${t}^{{\prime} }={\delta }^{{\prime} -1}t$, where t refers to the jet frame.

4.1. Physical Units and Synthetic Light Curves

As stressed, for the jet simulation we used the 3D general relativistic MHD code RAISHIN (Mizuno et al. 2006) that works with scale-free (code) units. For the comparison with the observations, L_j (t_s ) must be converted into physical units. To this aim, we have specified three physical constant units, velocity, time, and magnetic field (v₀, t₀, B₀), at the jet frame (corresponding to the simulation domain). Other important units are obtained from the previous ones, such as the length unit (r₀ = v₀ t₀), the density unit (${\rho }_{0}={B}_{0}^{2}/(4\pi {v}_{0}^{2})$), and the power unit (${L}_{0}={B}_{0}^{2}{t}_{0}^{2}{v}_{0}^{3}/4\pi$), all in centimeter-gram-second (cgs) units.

We have employed the light speed as the velocity unit (v₀ = c) since the RAISHIN code already works with this normalization. Then, taking the light curve of BL Lac Mrk 421 (see the top diagram in Figure 11; see also Kushwaha et al. 2017) we have defined two other physical units. Assuming that the time step between each snapshot of the simulation is equal to the minimum time variability observed from Mrk 421 at the same frame, the time unit is thus given by

$$\begin{eqnarray}&&{t}_{0}=\displaystyle \frac{{\delta }^{{\prime} }{\rm{\Delta }}{t}_{\mathrm{obs}}^{{\prime} }[\mathrm{cgs}]}{{\rm{\Delta }}t[{\rm{c}}.{\rm{u}}.]},\end{eqnarray} \tag{ 13 }$$

where Δt (at the jet frame) corresponds to the minimum time variability from the simulations, in code units, and ${\rm{\Delta }}{t}_{\mathrm{obs}}^{{\prime} }$ (at the observer's frame) is the observed one for Mrk 421, in cgs units. We have chosen the magnetic field unit B₀ so that the maximum photon flux of Mrk 421 coincides with the one obtained from the simulations (in physical units). Once defined B₀, we have obtained the unit power L₀, and converted L_j (t_s ) to cgs units (L[cgs] = L[c. u. ]L₀).

Zoom In Zoom Out Reset image size

As described previously, we have evaluated the luminosity in the observer's frame by the Doppler factor, and then estimated the synthetic photon flux (i.e., units of photons ${{\rm{c}}{\rm{m}}}^{-2}\,{{\rm{s}}}^{-1}$) in the 0.1–300 GeV energy range. ¹³ Table 3 shows the physical units for each variable obtained from this analysis taking into account different values of Doppler factors. Additionally, we included the corresponding proton density for ρ₀.

Table 3. Physical Constant Units (Centimeter-gram-seconds) in the Jet Frame

${\delta }^{{\prime} }$	t0 (s)	B0 (G)	r0 (cm)	ρ0 (g cm−3)	np ( cm−3)	L0 (erg s−1)
5.0	1.3 × 107	3.7	3.9 × 1017	1.2 × 10−21	717	4.8 × 1045
7.0	1.8 × 107	1.3	5.4 × 1017	1.5 × 10−22	92	1.2 × 1045
14	3.6 × 107	0.16	1.1 × 1018	2.3 × 10−24	1.3	7.2 × 1043
20	5.2 × 107	0.06	1.5 × 1018	3.2 × 10−25	0.19	2.1 × 1043

Download table as:  ASCIITypeset image

We used the jet model m240ep0.5 to build the synthetic light curve (see the bottom diagram in Figure 11) with a typical Doppler factor ${\delta }^{{\prime} }\sim 5$. This model shows an initial peak at t ∼ 30 (∼900 days at the observer frame) followed by a phase of low activity between t = 32 and 39. After this time, the activity increases with three prominent peaks at t ∼ 40, 44, and 63. This behavior is similar to the average reconnection rate profile (see Figure 8) described in Section 3.3. At the first peak, the average reconnection rate is not high (the maximum $\langle {\tilde{V}}_{\mathrm{rec}}{\rangle }_{s}$ value is around 0.02 at the reconnection frame), with a small number of events (see Figure 7), which indicates that each reconnection site identified by the algorithm releases considerable amounts of magnetic energy. Indeed, as can be seen in the second diagram in Figure 2, these events occur near the central axis of the jet where the magnetic field is high. The low activity phase, in turn, is associated to the decreasing number of reconnection events (as can be seen in Figure 7), with an intermediary peak at t ∼ 37. After t ∼ 39, the high activity seen in the synthetic light curve coincides with the nonlinear phase of the kink instability, where several reconnection sites are identified, with $\langle {\tilde{V}}_{\mathrm{rec}}\rangle$ up to 0.02. Model m240ep0.5 is interrupted at t ∼ 66 due to complete disruption of the jet. Therefore, we believe that the high activity phase observed in the light curve of Mrk 421 can be correlated with the nonlinear phase of the kink instability that drives turbulence and fast magnetic reconnection events (see also MGK21). Once the jet is locally disrupted, the activity is expected to decrease. If the kink instability rises again, the high activity can be resumed, starting another cycle. We will discuss in more detail the possible correlation between the time variability of Mrk 421 and reconnection events in the next section.

4.2. Power Spectral Density (PSD) Analysis

It is important to highlight that the aim of the analysis presented in the previous section was not to match exactly the synthetic light curve to the observations of Mrk 421. This because in this work, we did not perform background jet simulations specifically for this source, as they were intended to be more general. Nevertheless, this analysis has shown that we can already obtain reliable values, when converting the results of the simulations in physical units. The choice of Mrk 421 is based on the work of Kushwaha et al. (2017), where a careful study of this source was performed from a light-curve PSD analysis. The authors obtained the statistical proprieties of the time variability of Mrk 421, and concluded that it is broadly consistent with magnetic reconnection events, as described in the jets-in-a-jet model (see, e.g., Giannios et al. 2009; Christie et al. 2019). Therefore, we expect that the observed and synthetic light curves may show common features, such as the time variability.

We have then analyzed the PSD of the synthetic light curve and compared it with the one from Mrk 421 obtained by Kushwaha et al. (2017). We obtained the synthetic PSD from the reference model m240ep0.5 during a time range between t = 39 and 66, where the kink instability achieves the nonlinear phase (see also MGK21). Figure 12 shows such comparison for the observed (red line) and synthetic (black line) PSD in the observer's frame. The latter shows a broken power-law profile ¹⁴ with a frequency break of ${\mathrm{log}}_{10}({\nu }_{0}^{{\prime} })\sim -6.3$, in agreement with one obtained by Kushwaha et al. (2017). We have also obtained a compatible power-law index (a₂ ∼ −0.6) for frequencies higher than ${\nu }_{0}^{{\prime} }$. ¹⁵ On the other hand, the synthetic PSD is steeper (a₁ ∼ −2.0) than the observed one (a₁ ∼ −1.5) for frequencies lower than ${\nu }_{0}^{{\prime} }$. Since the synthetic light curve is a result of turbulent magnetic reconnection events, the similarity found between the PSDs emphasizes the argument that the time variability observed in Mrk 421 can be correlated with such events.

Zoom In Zoom Out Reset image size

In this work, we have evaluated how turbulent magnetic reconnection structures evolve in Poynting flux-dominated relativistic jets subject to the CDK instability. To this aim, we have performed 3D-SRMHD simulations following the work of Mizuno et al. (2012) to obtain suitable initial and boundary conditions for a small portion of a co-moving jet with an initial helical magnetic field and a moderate magnetization (σ ∼ 1.0) near the jet spine. Motivated by the previous work of Singh et al. (2016), we have examined the magnetic reconnection rates in such a system employing the algorithm developed in Kadowaki et al. (2018a), which was modified to take into account relativistic effects. This algorithm was able to perform a systematic search of all reconnection events in the system, allowing for visualization of the local magnetic topology of the current sheets, as well as quantification of the reconnection rate and the magnetic power released in these events, in three different frames (the reconnection, the jet, and the observer's frame). Our results can be summarized as follows:

• 1.  
  The identification of the magnetic reconnection sites reveals that these are initially concentrated in the jet spine. At the beginning of the growth of the twists of the magnetic field driven the CDK instability (between t ∼ 20 and ∼ 30 c. u.), the current density structures are thick. During this phase (after t ∼ 25), mostly sporadic, slow reconnection events with average rates smaller than $0.015{\tilde{V}}_{{\rm{A}}}$ are detected. As time evolves, after the instability grows exponentially up to saturation (around t = 40; see also MGK21), the jet wiggling structure amplitude continues to increase and becomes fully turbulent. The magnetic field lines are distorted and disrupted in several regions, and the associated current density structures become thinner (Figure 2), undergoing fast reconnection with the fastest rate around $0.23{\tilde{V}}_{{\rm{A}}}$ (at t = 56). The system achieves a turbulent, quasi-steady-state regime between t ∼ 50 and 66 with an average rate of $0.051\pm 0.026{\tilde{V}}_{{\rm{A}}}$ in agreement with the previous works of Takamoto et al. (2015) and Singh et al. (2016).
• 2.  
  The reconnection rate values resembles a log-normal distribution especially in the quasi-steady-state phase, after t = 50 (Figure 7). Since log-normal profiles are signatures of turbulence, this is another indication that turbulence is driving the fast reconnection events (in agreement with the theory of Lazarian & Vishniac 1999) along the relativistic jets. A χ² analysis reveals that asymmetric reconnection events produce a deviation from the log-normal distribution, which indicates that such events maybe not be associated with the same process.
• 3.  
  The convergence tests show no major differences in $\langle {\tilde{V}}_{\mathrm{rec}}{\rangle }_{s}$ between the models with the middle and high resolutions (240 and 480 cells in each direction), which implies that the fastest reconnection events are independent of the numerical resistivity, again in agreement with the turbulence-induced fast reconnection theory of Lazarian & Vishniac (1999; see also Kowal et al. 2009; Santos-Lima et al. 2021), which predicts that the reconnection does not depend on the ohmic resistivity of the environment.
• 4.  
  Our qualitative analysis of 2D cuts of the local reconnection events (using the LIC method) evidences X point-like structures of highly concentrated magnetic field lines (in a reconnection process) and magnetic islands (already reconnected) in several regions, where the magnetic field and the current density magnitudes are anticorrelated. These results are encouraging since they have demonstrated the efficiency of the algorithm at localizing events in different stages of reconnection. However, the 3D analysis has revealed that not always those magnetic islands are null point regions, where the magnetic field vanishes. Indeed, in the presence of the strong guide field (the helical magnetic structure), in some situations only one of the projected components is annihilated. This result indicates that analyses of the topology and evolution of blobs produced along relativistic jets due to reconnection events (such as, e.g., in Giannios et al. 2009; Christie et al. 2019) should be considered more carefully.
• 5.  
  A quantitative analysis of the highest velocity values around each reconnection region (in the reconnection frame) indicates an Alfvén speed value of ${\tilde{V}}_{{\rm{A}}}\sim 0.71c$, an inflow velocity of ${\tilde{V}}_{\mathrm{in}}\sim 0.075c$, and an outflow velocity of ${\tilde{V}}_{\mathrm{out}}\sim 0.28c$. Furthermore, the reconnection (diffusion) region velocities are mildly relativistic with respect to the jet frame, where the Lorentz factor is not higher than 1.09. Therefore, as we expected, a relativistic correction is appropriate for the Alfvén speed, but it will not change significantly the velocity field in the jet frame. As we have seen, such correction becomes important only when we consider an observer's frame where the jet frame (or the computational domain) is moving relative to it with high bulk Lorentz factors (${{\rm{\Gamma }}}_{j}^{{\prime} }\gt 10$).
• 6.  
  The synthetic light curve built from the integrated magnetic reconnection power and evaluated from the reference model m240ep0.5 (assuming a typical Doppler factor ${\delta }^{{\prime} }\sim 5$) shows a high activity phase after t ∼ 39 and coincides with the nonlinear phase of the growth of the CDK instability as well as the stage where several fast reconnection events are identified. The synthetic PSD extracted between t = 39 and 66 reveals a broken power-law profile with a frequency break of ${\mathrm{log}}_{10}({\nu }_{0}^{{\prime} })\sim -6.3$, in agreement with the observations of Mrk 421 (in the GeV band) obtained by Kushwaha et al. (2017). Those authors have also obtained a similar profile with PSD power-law indices a₁ ∼ −1.48 ± 0.20 (for ${\nu }^{{\prime} }\leqslant {\nu }_{0}^{{\prime} }$) and a₂ ∼ −0.57 ± 0.48 (for ${\nu }^{{\prime} }\gt {\nu }_{0}^{{\prime} }$). Though we have obtained for the lower frequency range a steeper power-law index for the synthetic PSD (a₁ ∼ −2.0), for the higher frequency range the value is very similar (a₂ ∼ −0.6). These results suggest that the turbulent fast magnetic reconnection driven by CDK instability is a possible process behind the daily time variability observed in the GeV band in the blazar Mrk 421, as suggested by Kushwaha et al. (2017). In future work, we have plans to apply a similar study to observed TeV fast variability emission of other blazars (e.g., PKS 2155-304, Mrk 501, 3C 279, and 3C 54.3), in order to provide hints not only on the variability but also about the broken power law. The origin of this profile is still not well understood in the observations and has been associated to the compact regions of the AGN sources and also in galactic X-ray binaries (see, e.g., McHardy et al. 2005). On the other hand, Kushwaha et al. (2017) have not found any evidence of a break scalability in the sources studied by them (NGC 1275, Mrk 421, B2 1520+31, and PKS 1510-089), which may indicate different origins. Thus, our numerical results can shed some light on this problem, strengthening reconnection models like, e.g., the jets-in-a-jet model (Giannios et al. 2009; Christie et al. 2019), or the kink turbulence-induced reconnection model described here (and in MGK21), or even the stripped jet model (Giannios & Uzdensky 2019; Zhang & Giannios 2021).

The implications of our findings, specially for the VHE emission observed in Poynting flux-dominated relativistic jets are rather important. The reconnection search algorithm employed here can, in principle, be applied to any simulated system, including those with higher Lorentz factors, magnetization, and variability such as the observed blazars with TeV flares with timescales of the order of hundred seconds (see, e.g., Aharonian et al. 2007; Albert et al. 2007a), or GRBs (e.g., Drenkhahn & Spruit 2002; Giannios & Spruit 2007; Zhang & Yan 2011). Employing this algorithm in a systematic study of 3D relativistic MHD jets subject to fast reconnection driven by, e.g., CDK as in this work, or other instabilities (see e.g., Nishikawa et al. 2020), or by mechanisms like the mini jets-in-a-jet model (e.g., Giannios et al. 2009), or the striped jet model (Giannios & Uzdensky 2019), will allow disentangling these different processes and constraining their applicability. However, as stressed, the algorithm has limitations, even with the sample constraints applied in the present work, especially in the identification of structures that are not reconnecting. More sophisticated methods can be used in future works. For instance, as discussed in Section 3.1, the time-evolution analysis of the magnetic structure images is a way to identify reconnection sites at different stages. Despite the high number of images necessary to use this method, machine-learning techniques can be applied to recognize and automatically select such events. Recently, Jafari et al. (2020) presented a new method based on the correlation between the magnetic field-fluid slippage and the system's stochasticity level, which can be evaluated by the scalar field ${\rm{\Psi }}=\displaystyle \frac{1}{2}{{\boldsymbol{B}}}_{l}\cdot {{\boldsymbol{B}}}_{L}$, where B _l and B _L correspond to renormalized magnetic fields at different scales (l and L). This method also represents a powerful tool for the identification of reconnection structures in turbulent environments and future work should confront both techniques.

Finally, as stressed, this study has also provided the basis for us to explore, in a companion work, in situ particle acceleration in the reconnection sites identified along the simulated 3D relativistic jet (see MGK21). We have found that thousands of protons injected with initial energies of 1 MeV are accelerated exponentially, predominantly parallel to the local magnetic field, achieving ultrahigh energies, which could explain the observed VHE and neutrino emissions in relativistic jets. This could be the case, for instance, of blazars like TXS 0506+056 (Aartsen et al. 2018), for which TeV gamma-rays and neutrinos have been simultaneously detected for the first time. In a forthcoming work, we plan to study in detail the radiative losses in these systems by incorporating radiative transfer effects to our model.

The analyses and 3D-SRMHD numerical simulations presented here were performed in the cluster of the Group of Plasmas and High-Energy Astrophysics (GAPAE), acquired with support from the Brazilian funding agency FAPESP (grant No. 2013/10559-5), and in the Blue Gene/Q supercomputer supported by the Center for Research Computing (Rice University) and Superintendência de Tecnologia da Informação da Universidade de São Paulo (USP). This work also made use of the computing facilities of the Laboratory of Astroinformatics (IAG/USP, NAT/Unicsul), whose purchase was also made possible by FAPESP (grant No. 2009/54006-4) and the INCT-A. L.H.S.K. also acknowledges support from FAPESP (grant No. 2016/12320-8), EMdGDP from FAPESP (grant No. 2013/10559-5) and the Brazilian agency CNPq (grant No. 308643/2017-8), and TEMT from CAPES. Y.M. is supported by the ERC Synergy Grant "BlackHoleCam: Imaging the Event Horizon of Black Holes" (grant No. 610058). P.K. acknowledges ARIES Aryabhatta postdoctoral Fellowship (AO/A-PDF/770). Also, we would like to thank the anonymous referees for their useful comments and discussions.