# Rotation and toroidal magnetic field effects on the stability of two-component jets

Dimitrios Millas, Rony Keppens, Zakaria Meliani

► To cite this version:

Dimitrios Millas, Rony Keppens, Zakaria Meliani. Rotation and toroidal magnetic field effects on the stability of two-component jets. Monthly Notices of the Royal Astronomical Society, 2017, 470, pp.592. 10.1093/mnras/stx1288. hal-01582447

HAL Id: hal-01582447

https://hal.science/hal-01582447v1

Submitted on 26 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Rotation and toroidal magnetic field effects on the stability of two-component jets

Dimitrios Millas,<sup>1★</sup> Rony Keppens<sup>1</sup> and Zakaria Meliani<sup>2</sup>

Accepted 2017 May 22. Received 2017 May 21; in original form 2017 March 27

## ABSTRACT

Several observations of astrophysical jets show evidence of a structure in the direction perpendicular to the jet axis, leading to the development of ‘spine and sheath’ models of jets. Most studies focus on a two-component jet consisting of a highly relativistic inner jet and a slower – but still relativistic – outer jet surrounded by an unmagnetized environment. These jets are believed to be susceptible to a relativistic Rayleigh–Taylor-type instability, depending on the effective inertia ratio of the two components. We extend previous studies by taking into account the presence of a non-zero toroidal magnetic field. Different values of magnetization are examined to detect possible differences in the evolution and stability of the jet. We find that the toroidal field, above a certain level of magnetization $\sigma$, roughly equal to 0.01, can stabilize the jet against the previously mentioned instabilities and that there is a clear trend in the behaviour of the average Lorentz factor and the effective radius of the jet when we continuously increase the magnetization. The simulations are performed using the relativistic MHD module from the open source, parallel, grid adaptive, MPI-AMRVAC code.

**Key words:** instabilities – magnetic fields – MHD – methods: numerical – galaxies: active – galaxies: jets.

## 1 INTRODUCTION

Since their discovery in the beginning of the 20th century, astrophysical jets have been extensively observed and studied on multiple scales, ranging from young stellar object jets (YSO jets) and gamma-ray bursts (GRBs) to active galactic nuclei jets (AGN jets). Improved instruments allowed for better observations of these peculiar outflows, resulting in the detection of energetic features (e.g. flares), improved resolution and mapping of interesting small-scale regions and information about their overall structure and evolution. The formation, acceleration and collimation of jets was addressed in the beginning analytically and (relatively) recently via numerical simulations, in different regimes, ranging from simple hydrodynamic (HD) to relativistic and general relativistic magnetohydrodynamics, with contemporary efforts focusing also on the effect and importance of radiation, e.g. Sądowski et al. (2013). AGN and GRB jets are relativistic, with Lorentz factors of $\gamma \sim 10$ and $\gamma \sim 100$, respectively, while YSO jets are non-relativistic, with typical velocities of $\sim 100 \text{ km s}^{-1}$ (Woitas et al. 2002).

Most of the scenarios for the large-scale jets involve some kind of accretion in a massive central object (i.e. a supermassive black hole) and the presence of a magnetic field. In general, the acceleration

of relativistic jets involves the transformation of magnetic energy (Poynting flux) to kinetic, as described with different mechanisms in different scales (Blandford & Znajek 1977; Blandford & Payne 1982), the efficiency of which is often associated with the shape of the poloidal field lines (Vlahakis 2004a,b). The (magnetic, self-) collimation of these outflows is attributed to the hoop stress of the toroidal magnetic field component (Bogovalov 1995), with recent studies focusing on the effect of an external medium, either ambient gas or wind-like outflows from the accretion disc etc., e.g. Globus & Levinson (2016), and references therein. The interaction with the external medium and causality arguments can possibly explain the stability of jets as mentioned in Porth & Komissarov (2015).

Analytical work and simulations, in various frameworks and using different physical modules, have been carried out to test the theoretical predictions concerning the launching, acceleration and collimation of jets. For the AGN case, we mention e.g. Casse & Keppens (2004) and Tchekhovskoy, Narayan & McKinney (2011), concerning the launching of jets and the work of Komissarov et al. (2007), where the MHD acceleration of jets was examined. Other studies focused on the development of instabilities in jets, (e.g. Baty & Keppens 2002; Bodo, Mignone & Rosner 2004; Mizuno et al. 2012; Bodo et al. 2013; Mizuno, Hardee & Nishikawa 2014).

It has also been argued that astrophysical jets are not always homogeneous but rather display a structure in the direction perpendicular to the jet axis. This structure has been observed in terms of

<footnote n="1"><sup>1</sup>Centre for Mathematical Plasma Astrophysics, Department of Mathematics, KU Leuven, 3001, Leuven, Belgium</footnote>

<footnote n="2"><sup>2</sup>LUTH, Observatoire de Paris-Meudon, F-92195, Meudon, France</footnote>

<footnote>* E-mail: dimitrios.millas@kuleuven.be</footnote>
velocity (Boccardi et al. 2016; Mertens & Lobanov 2016) and is believed to exist in different scales, ranging from YSO to AGN jets. Especially for the AGN jet case, the implications on the radiation emitted from specific sources have been investigated (Giroletti et al. 2004; Ghisellini, Tavecchio & Chiaberge 2005). The above studies resulted in the development of ‘spine and sheath’ models (e.g. Gabuzda, Reichstein & O’Neill 2014), where a jet is believed to consist of a fast inner part (spine) and a slower outer part (sheath).

The effects of this structure have also been considered in some analytical studies (e.g. Bogovalov & Tsinganos 2001) and via numerical simulations in 2.5D (e.g. Meliani & Keppens 2007, 2009), concerning the axisymmetric stability, and in 3D (Singh, Mizuno & de Gouveia Dal Pino 2016), concerning the transformation of magnetic to kinetic energy due to kink instability. These papers referred to AGN jets whereas Matsakos et al. (2008, 2009) and Teşileanu et al. (2014) examined a two-component YSO jet.

Our primary task in this paper is to examine the effect of this structure to the transverse jet stability and specifically focus on non-axisymmetric instabilities induced by differential rotation, as described in Meliani & Keppens (2009). We also note that Matsumoto & Masada (2013) argued that in the absence of magnetic field, Rayleigh-Taylor instabilities might still occur in a relativistic jet.

The assumption of a purely poloidal magnetic field will be relaxed, introducing a helical field, while the translational symmetry along the axis of the jet will be maintained for the 2.5D setups. The existence of a toroidal magnetic field component is suggested by Gabuzda et al. (2014) and Gabuzda, Knuettel & Reardon (2015), where different Faraday rotation measures (RM) were examined for a selection of AGN jets. The main part of our 2.5D work will involve a parametric study using different values of magnetization and some extreme scenarios involving a slower inner jet and a slowly rotating jet. The 3D work will only be limited to some extreme cases based on the results obtained by our 2.5D simulations. We assume that the region in question is located sufficiently far from the central source, so the jet is already accelerated and collimated as described in the references above. Therefore, we will only examine cases of kinetically dominated jets. The aim is to see if a helical magnetic field is able to stabilize the jet against instabilities induced by differential rotation.

The results of this research may be of use in the FRI/FRII dichotomy (Fanaroff & Riley 1974), concerning the morphology of radio jets. Our cases will use standard parameters suitable for AGN jets, but in principle one can examine a similar problem in other astrophysical objects known to host jets (e.g. GRBs).

## 2 INITIAL STATE OF THE JET

We will follow the same recipe for the overall configuration of the jet and the normalization as described in Meliani & Keppens (2009). Velocity is normalized in units of c (and we set c = 1), distance is measured in pc and mass is normalized to proton mass $m_p$.

First, we fix the (outer) radius of the jet at $R_{out} = 0.1$ pc, and the inner radius is arbitrarily chosen to be $R_{in} = R_{out}/3$. These values are consistent with Biretta, Junor & Livio (2002), where the opening angle of the M87 jet was found ~10° at a distance of ~4 pc and the jet is almost fully collimated at the same scale. The next step is to define the initial conditions for the density, velocity, magnetic field and pressure profiles of the jet (inner and outer regions) and the surrounding, external medium.

We examine a two-component jet in which the two parts of the outflow are differentially rotating, with an initial toroidal velocity

profile set as:

$$V_{\phi}(R) = \begin{cases} v_{\phi in} \left(\frac{R}{R_{in}}\right)^{\alpha_{in}/2}, & R \le R_{in} \\ v_{\phi out} \left(\frac{R}{R_{in}}\right)^{\alpha_{out}/2}, & R_{in} < R < R_{out} \end{cases}, \quad (1)$$

where we assume $\alpha_{in} = 0.5$ for the inner jet and $\alpha_{out} = -2$ for the outer jet. In contrast with Meliani & Keppens (2007, 2009), where a profile with a jump at the interface of the two components was used, we assume that in our case the toroidal velocity component is continuous at $R = R_{in}$. The continuity of the toroidal velocity profile is implied by Mertens et al. (2016). Initially, we will impose this condition by selecting $v_{\phi in} = v_{\phi out} = 0.01$ for the main part of our work, but later on we will also examine cases in the non-rotating limit, using $v_{\phi in} = v_{\phi out} = 10^{-6}$.

The dominant velocity component of such jet is poloidal and set to a constant value, different for each component. Again the main part of our work will assume $v_z$ with a corresponding Lorentz factor of $\gamma_{z,in} = 30$ for the inner jet and $\gamma_{z,out} = 3$ for the outer jet, which are meaningful values for AGN jets (Giroletti et al. 2004). Later on, we will decrease the Lorentz factor of the inner jet to $\gamma \simeq 10$. As the toroidal velocity component is significantly smaller than the poloidal, the above values are approximately equal to the total Lorentz factor (see Fig. 1c).

For the magnetic field, we assume a uniform poloidal component, constant in each part of the jet ($b_{z,in}, b_{z,out}$) and a toroidal component of the same form as in equation (1):

$$B_{\phi}(R) = \begin{cases} b_{\phi in} \left(\frac{R}{R_{in}}\right)^{\alpha_{in}/2}, & R \le R_{in} \\ b_{\phi out} \left(\frac{R}{R_{in}}\right)^{\alpha_{out}/2}, & R_{in} < R < R_{out} \end{cases}, \quad (2)$$

where the constants $b_{\phi in}, b_{\phi out}$ are defined by fixing the maximum magnetization of the jet, $\sigma = B_{\phi}^2/(\gamma^2 \rho)$ at $R = R_{in}$. Similarly to the definition of the toroidal velocity, we select $b_{\phi in} = b_{\phi out}$, while the poloidal magnetic field will still be discontinuous at $R = R_{in}$. The initial 3D reconstruction of the field lines of this helical magnetic field is shown in Fig. 1(d). For completeness, we include the profiles of $V_z$ and $B_z$ with $R$:

$$V_z(R) = \begin{cases} v_{z,in}, & R \le R_{in} \\ v_{z,out}, & R_{in} < R < R_{out} \end{cases}, \quad (3)$$

$$B_z(R) = \begin{cases} b_{z,in}, & R \le R_{in} \\ b_{z,out}, & R_{in} < R < R_{out} \end{cases}, \quad (4)$$

where $v_{z,in}, v_{z,out}, b_{z,in}, b_{z,out}$ are constants.

The jet is surrounded by a static external medium, with a density roughly of the same order of magnitude as of the inner jet. We set the number density of the surrounding medium to $n_{med} = 10^{-2} \text{ cm}^{-3}$. The jet density values can be calculated using constraints from observations, namely the kinetic luminosity flux for a typical radio loud galaxy (~$10^{46}$ erg s$^{-1}$) and the ratio of the flux carried by each jet component. If the inner part carries only 1 per cent of the flux, then the density values for the inner and outer jet are $\rho_{in} = 6.92\rho_{med}$ and $\rho_{out} = 119.94 10^2 \rho_{med}$, resulting in an initial density ratio of ~$10^4$ between the components. In all cases, the jets are over-dense. This is related to the fact that we are modelling the jet at a sub-parsec region, where the external medium could result from a disc wind and/or a hot medium in the inner region of the AGN. The chosen value of density for the surrounding medium is also reasonable for a jet cocoon (Meliani, Keppens & Giacomazzo 2008).
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Clockwise from top left-hand plot: Initial conditions for density, toroidal velocity, magnetic field lines and Lorentz factor for the inner and outer part of the outflow. Concerning the field lines, the colour scale represents the magnitude of the magnetic field, so the inner jet is displayed in 'warm' colours, the outer jet in 'cold', while the external medium is unmagnetized.</figcaption></figure>

The density is constant in each part of the jet (Fig. 1a and equation 5).

$$ \rho(R) = \begin{cases} \rho_{\text{in}}, & R \le R_{\text{in}} \\ \rho_{\text{out}}, & R_{\text{in}} < R < R_{\text{out}} \\ \rho_{\text{med}}, & R > R_{\text{out}} \end{cases}. \qquad (5) $$

The relativistic equivalent of Rayleigh's criterion for rotational stability is that the angular momentum flux must increase with the radial distance R. The angular momentum flux is given by the formula:

$$ I = \gamma \frac{\rho + \frac{\Gamma}{\Gamma - 1}}{p} \rho V_{\phi} R - \frac{B_p}{\gamma \rho V_p} R B_{\phi}, \qquad (6) $$

and thus if $B_\phi = 0$ and $V_\phi$ as in equation (1), the inner jet is stable (as $\frac{dI}{dR} > 0$) and the outer jet is marginally stable, as $\frac{dI}{dR} = 0$. The interface between the two components is still unstable, as the angular momentum changes significantly at $R = R_{\text{in}}$. In our study, all $B_z$ values still satisfy the criterion.

We assume in addition total pressure equilibrium between the two interfaces (inner and outer jet, outer jet and external medium) at $t = 0$. The steady-state momentum equation is (e.g. Heyvaerts & Norman 2003):

$$ \gamma \rho (v \cdot \nabla)(\xi \gamma v) = -\nabla P + J \times B + \rho_e E, \qquad (7) $$

which is equivalent to equation (23) in Mobarry & Lovelace (1986), ignoring the geometric/gravity terms. Assuming a polytropic equation of state, the differential equation describing the variation of the total pressure $P_{\text{tot}}$ along the radial direction is

$$ \frac{dP_{\text{tot}}}{dR} - \frac{\Gamma}{\Gamma - 1} \frac{\gamma^2 V_{\phi}^2}{R} P_{\text{tot}} = \frac{\gamma^2 V_{\phi}^2}{R} \left( \rho - \frac{\Gamma}{\Gamma - 1} \frac{B_z^2}{2} \right) + \frac{1}{R} \left[ 1 + \frac{\Gamma \gamma^2 V_{\phi}^2}{2(\Gamma - 1)} \right] \left[ -B_{\phi}^2 + (B_{\phi} V_z - V_{\phi} B_z)^2 \right], \qquad (8) $$

where $P_{\text{tot}} = p + \frac{b^2 - \Gamma^2}{2}$ is the total pressure. Using equations (1), (2) and (8), we find that the thermal pressure $p$ varies with $R$ as:

$$ p = \zeta \left[ 1 - \tilde{a} \left( \frac{R}{R_{\text{in}}} \right)^{\alpha} \right]^{-\frac{\Gamma}{\alpha(\Gamma - 1)}} - \frac{\Gamma - 1}{\Gamma} \rho - \frac{(\Gamma - 1)(\alpha + 2)}{2\tilde{a}[\alpha(\Gamma - 1) + \Gamma]} \left[ -b_{\phi}^2 + (b_{\phi} V_z - v_{\phi} B_z)^2 \right] \times \left( 1 - \tilde{a} \left( \frac{R}{R_{\text{in}}} \right)^{\alpha} \right), \qquad (9) $$

where $\tilde{a} = \frac{B_{\phi}^2}{1 - v_z^2}$ and the parameter $\alpha$ is $\alpha_{\text{in}}$ and $\alpha_{\text{out}}$ for the inner and outer jet, respectively. The integration constant $\zeta$ ($\zeta_{\text{in}}$, $\zeta_{\text{out}}$, depending on the component) is obtained from the boundary conditions (the thermal pressure on the axis, $p_o$) and the total pressure matching at the interface of the two components. The full expressions for these constants can be found in the appendix. The effective polytropic index $\Gamma$ is initially approximately 4/3 for the inner jet and the external medium and 5/3 for the outer jet, which correspond to relativistically hot and cold outflows, respectively.

Following Meliani & Keppens (2009), we will later use a Synge equation of state for the simulations, since mixing of matter from both parts of the jet is present. This does not contradict the earlier polytropic approximation since $\Gamma_{\text{eff,in}}$ and $\Gamma_{\text{eff,out}}$ are constant in the initial state.

### 3 NUMERICAL SETUP AND DIFFERENT CASES

We will examine different configurations of two-component jets, surrounded by a static medium, including a toroidal magnetic field component in both parts of the jet. This results in non-zero magnetization for both components, with each sub-case using a different
value for the maximum magnetization. All the chosen values of magnetization, although underestimated, correspond to a kinetically dominated jet, as $\sigma \le 0.1$, which is in agreement with the assumption that far from the central engine, AGN jets have already been accelerated and collimated.

We select values of magnetization ranging from $\sigma = 0.001$ to $\sigma = 0.1$. Every 2.5D run is performed on a Cartesian domain with dimensions $-0.3 < x, y < 0.3$ pc. We use the relativistic MHD module of the MPI-AMRVAC code (Keppens et al. 2012; Porth et al. 2014) selecting an HLLC solver and a third-order limiter (Čada & Torrilhon 2009). The base resolution of each simulation is $128 \times 128$ with three levels of adaptive mesh refinement (AMR) used, achieving an effective resolution of $512 \times 512$ (thus covering $\sim 0.001$ pc per cell in the finest resolution). Higher resolution runs were also performed, where fine structure effects naturally appeared, with the overall image remaining the same. Every simulation covers a time range of $t = 60$ or 195.8 yr, which is roughly equal to three rotations of the inner jet. The boundary conditions are open in every side of the computational box, with restrictions clipping any inflow. We also use two different tracers, one for each jet component, to detect any mixing between the two parts of the outflow.

The 2.5D simulations can be categorized as follows: (i) simulations for different values of $\sigma$ where we chose for the inner jet $\gamma = 30$, (ii) cases using the lowest and highest values of magnetization values of (i) with $\gamma = 10$ for the inner jet, and (iii) Cases using the lowest and highest values of magnetization assuming a very slow rotation for the jet. In all our 2.5D cases, the thermal pressure on the axis is set to $p_o = 2$. When the initial Lorentz factor of the inner jet is reduced to $\gamma = 10$, the total kinetic luminosity flux is reduced by less than $\sim 1$ per cent. Using the same configuration for the outer jet, the density ratio between the components remains constant and thus the main contribution to the kinetic flux is still due to the outer jet, in both cases.

For the 3D simulations, we chose periodic boundary conditions on the $z$ axis, while the computational domain is a Cartesian box with dimensions $-0.3 < x, y < 0.3$ pc, $0 < z < 1$ pc and a resolution of $256^3$, without AMR. We examine the most interesting cases of the 2.5D scenarios, as explained in the relevant section. The thermal pressure on the axis is set to $p_o = 2.3$. For simplicity, we consider no scaling with $z$ for all physical quantities. Although most studies assume $B_\phi \sim 1/z$ and $B_z \sim 1/z^2$, we expect little difference in the small scale examined here. Our 3D simulations therefore correspond to a cylindrical jet.

In some cases, we notice initially the formation of four 'arms', which is due to the fact that we use a Cartesian grid and thus modes with $m = 4$ are being favoured. This can be suppressed by exciting specific modes in the perturbations of the radial velocity. Two different types of perturbations are used, either with $\nabla \cdot v$ or perturbations with selected $m$ (Rossi et al. 2008). The kind of perturbation will be mentioned in each different case.

## 4 INCREASING THE MAGNETIZATION

We first present the results for the cases with a Lorentz factor of $\gamma = 30$ for the inner jet, where we follow the evolution of the density and the Lorentz factor with time. A wide range of different magnetization values is used and the key values are summarized in Table 1. We will present in more detail three representative cases, with $\sigma = 0.001, 0.01$ and $0.1$, while the rest will be used to examine the existence of trends in the evolution of the average Lorentz factor and the effective radius of the jet with $\sigma$. A comparison of the

<table><caption>**Table 1.** Key parameters for the examined cases. The average Lorentz factor for the inner jet and the effective radius of the entire jet after three rotations are given in the last two columns.</caption><thead><tr><th>Case #</th><th>σ<sub>max</sub></th><th>γ<sub>in,0</sub></th><th>ν<sub>φ</sub></th><th>B<sub>z,in</sub></th><th>B<sub>z,out</sub></th><th>γ<sub>in,avg</sub></th><th>R<sub>eff,jet</sub> (pc)</th></tr></thead><tbody><tr><td>1</td><td>0.001</td><td>30</td><td>0.01</td><td>1.228</td><td>1.039</td><td>10.9</td><td>0.144</td></tr><tr><td>2</td><td>0.005</td><td>30</td><td>0.01</td><td>0.549</td><td>1.039</td><td>13.5</td><td>0.134</td></tr><tr><td>3</td><td>0.01</td><td>30</td><td>0.01</td><td>0.388</td><td>1.039</td><td>15.4</td><td>0.129</td></tr><tr><td>4</td><td>0.05</td><td>30</td><td>0.01</td><td>0.174</td><td>1.039</td><td>21.4</td><td>0.116</td></tr><tr><td>5</td><td>0.1</td><td>30</td><td>0.01</td><td>0.123</td><td>1.039</td><td>22.2</td><td>0.116</td></tr><tr><td>6</td><td>0.001</td><td>10</td><td>0.01</td><td>0.429</td><td>1.039</td><td>7.5</td><td>0.125</td></tr><tr><td>7</td><td>0.1</td><td>10</td><td>0.01</td><td>0.043</td><td>1.039</td><td>9.8</td><td>0.12</td></tr><tr><td>8</td><td>0.001</td><td>10</td><td>10<sup>-6</sup></td><td>0.429</td><td>1.039</td><td>27.0</td><td>0.132</td></tr><tr><td>9</td><td>0.1</td><td>10</td><td>10<sup>-6</sup></td><td>0.043</td><td>1.039</td><td>26.7</td><td>0.110</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Toroidal magnetic field magnitude for five different values of magnetization $\sigma$, also mentioned in Table 1.</figcaption></figure>

strength of the toroidal magnetic field, in normalized units for our work, for different values of $\sigma$ is shown in Fig 2. The initial state of the jet in terms of density is the same for all simulations and thus it will be presented only once in the first case.

### 4.1 Case 1: $\sigma = 0.001$, $\gamma_{in} \approx 30$

From a physical point of view, we expect this case to be the most unstable one, since the hoop stress is minimum and the small magnetization value makes it effectively a (relativistic) HD set-up, somewhat similar to Case D studied in Meliani & Keppens (2009) (since $B_z$ only contributes to the total pressure). In this case, we follow equation 3 of Rossi et al. (2008) in the choice of the initial radial velocity perturbation, choosing $m = 3$, while $\omega_l, b_l, l=0$.

Initially, the centrifugal force dominates over the magnetic tension, a fact which is reflected in the quick expansion of the surface of the jet (both for the inner and the outer part), shown in Fig. 3. Even though the toroidal velocity transition is continuous at $R_{in}$ – and not a jump, as in the original study of Meliani & Keppens (2009), some Kelvin–Helmholtz instabilities develop in the early stages of the simulation, but are quickly overcome by a Rayleigh–Taylor-type instability. In the inner part, we notice the formation of ‘arms’ quite early (around half rotation time of the inner jet), a feature which persists throughout the evolution of the jet. Apart from the inner jet material getting mixed in the outer jet by Rayleigh–Taylor formation, there is also inflow from the outer component into the inner jet region. This inflow is of high density and low angular momentum, which contributes to the deceleration of the jet. The mixing starts after $\sim 0.25$ rotations of the inner jet, which becomes clear after $\sim 1$.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Proper density (log) for Case 1 (maximum $\sigma = 0.001$, $\gamma_{in} \simeq 30$). Snapshots after 0, 0.5, one, two and three rotations of the inner jet.</figcaption></figure>

rotation. In the end, a mixed region is formed in the central part, while the high-density outer part of the jet is pushed outwards and the jet as a whole expands. The density ratio between the final state and the initial value for the inner jet can reach the order of ~100.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Lorentz factor for Case 1 (maximum $\sigma = 0.001$, $\gamma_{in} \simeq 30$) after three rotations of the inner jet.</figcaption></figure>

The dominant instability in this simulation is thus a relativistically enhanced Rayleigh–Taylor-type, as discussed in Meliani & Keppens (2009), with Kelvin–Helmholtz effects remaining present mostly early in the simulation and restricted in the boundary between the two jet components. The magnetic tension remains weak throughout the evolution and is not sufficient to constrain the jet expansion.

The distribution of Lorentz factor in the final state is given in Fig. 4. The jet clearly de-collimates after one rotation time, finally reaching an effective radius of $R_{eff} \simeq 0.144$ pc and the average Lorentz factor for the inner jet is decreased to ~10.9.

## 4.2 Case 3: $\sigma = 0.01$, $\gamma_{in} \simeq 30$

We increase the maximum magnetization by a factor of 10 with respect to the first case, to $\sigma = 0.01$. In this case, we do not excite specific modes but rather use a zero-divergence perturbation in the radial velocity. The first stages of the evolution resemble the previous case, at least up to ~0.5 rotations, with low density material from the inner jet moving outwards, due to the centrifugal force and dense, low angular momentum material falling inwards. We now notice the numerical selection of four arms effect, favoured by our Cartesian grid, seen up to ~1 rotation time of the inner jet. Later, after ~ two rotations of the inner jet, the central part of the jet consists of a well-mixed region. After the formation of the mixed region, the ‘high’ and ‘low’ density components of the outflow can be distinguished once more, with different values compared to the initial state. Although what remains of the outer component does not experience any significant change in density, the inner jet remains in a mixed state, with a density of ~50 times greater than the initial inner jet density. The evolution of the density with time can be seen in Fig. 5. We also present the distribution with the radial distance, e.g. along the x-axis. As seen in the previous case, the problem is not axisymmetric, but we can still obtain qualitative information about the mixing.

The increased hoop stress prevents in a sense the development of evident Rayleigh–Taylor-type instabilities in large scale, whereas we still notice the development of Rayleigh–Taylor-type instability between the outer jet and the environment. As before, any Kelvin–Helmholtz phenomena between the two jet components are quickly suppressed.

The Lorentz factor distribution after three rotation times is shown in Fig. 6. Observing the evolution of the Lorentz factor, we notice that even though the effective radius of the inner jet increases, an inner, fast and an outer, slow region can be clearly distinguished in the outflow for the entire evolution shown. The average Lorentz
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Proper density (log) for Case 3 (maximum $\sigma = 0.01$, $\gamma_{in} \simeq 30$). Snapshots after 0.5, one, two and three rotations of the inner jet. The bottom image shows the distribution along the x axis at $t = 3$ (i.e. 195.8 yr).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Lorentz factor for Case 3 (maximum $\sigma = 0.01$, $\gamma_{in} \simeq 30$) after three rotations of the inner jet.</figcaption></figure>

factor after three rotations is 13.5 and the effective radius of the jet $R_{eff} \simeq 0.134$ pc.

### 4.3 Case 5: $\sigma = 0.1$, $\gamma_{in} \simeq 30$

Here, we examine the case with the maximum value of magnetization in our study, $\sigma = 0.1$. Following the evolution of the previous cases, we expect this configuration to be the most stable one, due to the increased magnetic tension. We use the same kind of perturbation as in the case with $\sigma = 0.01$.

As seen earlier in this work, in the early stages we notice the formation of ‘arms’ and some Kelvin-Helmholtz effects between the two jet components. In this case, however, the mixing starts after ~0.5 rotation and is restricted near the interface of the two jet components. The ‘arms’ persist up to ~0.75 rotation of the inner jet and later form a mixed region in the edge of the inner jet. The mixing is now present in a quite smaller region compared to the previous cases, with the final state resulting in a ‘fast jet’ area smaller than the one found in the unstable cases. Any enhanced Rayleigh-Taylor-type phenomenon is quickly (before ~1.5 rotations of the inner jet) suppressed.

The final state is in a sense similar to the previous case due the strong magnetic tension, resulting again in a fast, inner and a slow, outer component. We note that even in this case, the inner jet has a higher density compared to the initial state due to the mixing process, with a numerical value ~20 times greater than the initial value for the inner jet. The distribution of the density over time and the distribution along the x axis can be found in Fig. 7. The Lorentz factor after three rotations is given in Fig. 8.

The outer jet, in terms of density, remains unaffected and after three rotations is still more or less uniform. The effective radius of the inner jet is also in this case increased but the two components are more structured and distinguishable. The average Lorentz factor of the inner jet after three rotations is ~22.2 and the effective radius of the whole outflow is $R_{eff} \simeq 0.116$ pc.

### 4.4 Average lorentz factor, effective radius and mixing

Here, we discuss the evolution of the average Lorentz factor and the effective radius of the jet, inner and outer, for cases with increasing magnetization, some of which were discussed in Sections 4.1–4.3. We also assess the effect of mixing in the effective polytropic index and the density of each jet component.

In the absence of a toroidal magnetic field, it has been shown that the average Lorentz factor is significantly reduced in all unstable cases, while the effective radius of the jet (inner and outer) increases (Meliani & Keppens 2009). This trend persisted, less obvious, also
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Proper density (log) for Case 5 (maximum $\sigma = 0.1$, $\gamma_{in} \simeq 30$). Snapshots after 0.5 and three rotations of the inner jet. The bottom image shows the distribution along the *x* axis at *t* = 3 (i.e. 195.8 yr).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Lorentz factor for Case 5 (maximum $\sigma = 0.1$, $\gamma_{in} \simeq 30$) after three rotations of the inner jet.</figcaption></figure>

in more stable set-ups. We expect a similar behaviour in every case we examined here as well, with a less prominent decrease of the Lorentz factor at least in the very stable cases. The effective radius of each component is expected to increase in every case as well (as radial expansion is observed in every simulation), with the more stable cases resulting in a somewhat smaller increase (at least for the inner part).

First, we define the criteria that we use to distinguish the two components. Specifically, we consider as inner jet the region in which the Lorentz factor has a value of $\gamma > 5$ and the density is $\rho < 1$, whereas the whole jet is distinguished from the surrounding medium if $\gamma > 1.5$ and $B_z > 0.01$. We note however that the precise final result can be quite sensitive to the selection of the above criteria, especially when strong mixing is taking place. The value of the Lorentz factor is sufficient though to mask each region in all our cases so far.

For very low magnetization values, the average Lorentz factor is significantly reduced, as expected since we asymptotically approach the case of a purely poloidal magnetic field. We see a clear trend in the behaviour of the average Lorentz factor of the whole jet with increasing magnetization, as seen in Fig. 9. As $\sigma$ becomes higher, the average Lorentz factor decreases at a slower rate and the effective radius of the whole jet also increases at a slower rate. For the most stable cases, the increase in the effective radius of the jet compared to the initial state is of the order of 10 per cent.

In a similar way, we see a trend also in the evolution of the effective radius of the whole jet with $\sigma$. Increasing $\sigma$ results in a less significant expansion and a more structured outflow, for both components. The effective radius of the inner jet is affected by mixing, which slows down this part of the outflow. Due to the previous criteria, in cases where mixing is stronger (in other words, in the most unstable setups), a smaller region is identified as ‘inner jet’. This deceleration is less prominent in stable cases. We notice a peculiar behaviour in the evolution of the case with $\sigma = 0.05$ (Fig. 9, second panel, left-hand column). This is due to a transient state that results in an oval-like shape of the inner jet, present approximately from ~1.5 to 2.5 rotations.

Last, since the two jet components have different effective polytropic indices, it is interesting to check if the mixing has any effects in the final effective polytropic index of each region. For the unstable case (with $\sigma = 0.001$), we notice that the effective polytropic index is significantly modified after three rotations. In the other two cases, with $\sigma = 0.01$ and 0.1, we note that the outer jet is more or less unaffected, as noted already in terms of density, and the effective polytropic index is the same as in the initial state. The inner jet has a higher value of $\Gamma_{eff}$ compared to the initial state, but still closer to the relativistically hot value (especially for $\sigma = 0.1$).

## 5 INNER JET EVOLUTION AT LOWER LORENTZ FACTORS

In this section, we examine the differences in the evolution of the jet if we assume a lower Lorentz factor for the inner jet (closer to the values suggested by Piner et al. 2003; Giroletti et al. 2004; for radio galaxies). We reduce the poloidal velocity of the inner jet so that $\gamma \simeq 10$, and we examine two values of magnetization, $\sigma = 0.001$ and $\sigma = 0.1$.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Left-hand column: (a) average Lorentz factor of the inner jet, (c) effective radius of the inner jet after three rotations of the inner jet and (g) final average Lorentz factor with $\sigma$. Right-hand column: (b) initial effective polytropic index (for all cases) and effective polytropic index after three rotations for (d) $\sigma = 0.001$, (e) $\sigma = 0.01$, (f) $\sigma = 0.01$ and (h) $\sigma = 0.1$ (i).</figcaption></figure>

## 5.1 Case 6: $\sigma = 0.001$ and $\gamma_{in} = 10$

This is a direct variant of Case 1, with a modified value for the Lorentz factor of the inner jet, $\gamma \simeq 10$.

Initially, the evolution follows what has already been seen in the previous cases, with some notable differences. First, we can mostly distinguish the outer, high density part and the light, inner part of the jet throughout the evolution of the set-up. Kelvin-Helmholtz
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Proper density (log) for Case 6 (maximum $\sigma$ = 0.001). In this case, the initial value of the Lorentz factor for the inner jet is 10. Snapshots after two and three rotations of the inner jet.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Lorentz factor for Cases 6 (left-hand panel, maximum $\sigma$ = 0.001) and 7 (right-hand panel, maximum $\sigma$ = 0.1) after three rotations of the inner jet. Both cases assume an initial Lorentz factor of 10 for the inner jet. We note that in both runs, the final state is much more stable compared to Cases 1 and 5. Especially for Case 7, the change, compared to the initial state, is negligible.</figcaption></figure>

instabilities can first be seen between the two jet components at a very early time (<0.5 rotations of the inner jet and then vanishing), whereas Rayleigh–Taylor modes are mostly suppressed. In terms of Lorentz factor, both components remain structured and collimated after three rotations. The density and Lorentz factor distributions are shown in Figs 10 and 11, respectively.

Although Case 1 was the most unstable one, in terms of deceleration and decollimation, this case is considerably more stable, resembling in a sense the cases of higher magnetization of Section 4. The stability is related to the decrease of the inner jet inertia, since this is proportional to $\gamma^2$, as discussed later in the paper.

## 5.2 Case 7: $\sigma = 0.1$ and $\gamma_{in} = 10$

In a similar fashion, we examine the case with highest magnetization (Case 5, $\sigma$ = 0.1) in a jet where we substitute again $\gamma_{z,in} \simeq 10$.

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Proper density (log) for Case 7 (maximum $\sigma$ = 0.1). In this case the initial value of the Lorentz factor for the inner jet is 10. Snapshots after 2 and 3 rotations of the inner jet.</figcaption></figure>

Case 5 was the most stable, so we expect this feature to persist also here. Initially, there is little interaction between the two jet components.

Even at later times, the inner and outer parts of the jet experience little to almost no mixing and the two components can be clearly distinguished throughout the evolution of the jet, both in terms of Lorentz factor and density. The increased stability is again attributed to the lower effective inertia contrast between the components. For the density and Lorentz factor distributions, we refer to Figs 12 and 11.

We note that in both cases, lower inner jet inertia leads to more stable jets with respect to the Rayleigh–Taylor-type instability, which is in agreement with the (HD or purely poloidal magnetic field) results of Meliani & Keppens (2009). This is not always the case though in other types of instabilities, e.g. Kelvin–Helmholtz, which are stabilized for higher values of jet velocity. This is shown in the case of a single component, matter dominated, cold jet in Bodo et al. (2013).

# 6 THE NON-ROTATING JET LIMIT

In this section, we examine the evolution of the cases with $\sigma$ = 0.001 and 0.1 (keeping $\gamma_{in}$ ≃ 30) in the non-rotating limit, choosing a very slow value for $v_\phi$, namely $v_\phi$ = $10^{-6}$. This effectively cancels the centrifugal force and greatly stabilizes the outflow. Although in 2.5D, we cannot monitor Kelvin–Helmholtz instabilities along the z axis, which is dominant in non-rotating matter dominated jets according to Bodo et al. (2013), we can check if the absence of rotation ceases to induce Rayleigh–Taylor-type instabilities.

## 6.1 Case 8 and 9: $\sigma = 0.001$ and $\sigma = 0.1$

A notable difference in the non-rotating jet limit is the absence of ‘arms’ in the early stages of the evolution (up to ~0.5 rotations).

Even for the lowest value of magnetization, $\sigma$ = 0.001, the centrifugal force is practically zero to play any role in the evolution of
the outflow. A limited inward and outward motion can be noticed between the two parts of the jet, which results in the formation of a small annulus of mixed material that surrounds the inner jet.

In the case of $\sigma = 0.1$ almost no mixing is observed. The two components remain separated throughout the run, clearly seen in the density plots. The average Lorentz factor at the end state is very similar to Case 7. This is an indication that in all cases where we assumed a rotating jet, the centrifugal force plays a key role in the mixing.

The density distributions for each case are given in Figs 12 and 13, respectively. In Fig. 12, we see that for $\sigma = 0.001$, the outer jet still shows interchange/Rayleigh-Taylor-type effects.

## 6.2 Average lorentz factor and effective radius

As stated before, for very stable cases, we expect the average Lorentz factor to be more or less constant throughout the evolution of the jet. The same should hold for the effective radius, where we expect the outflow to expand less than in the unstable cases.

In the non-rotating limit, this holds, with a small increase in the effective radius of the inner jet due to the presence of a small shear region between the components. In principle, using more strict criteria (which may not be appropriate for the previous cases though), one can perform some fine tuning to clearly distinguish all different regions.

The effective radius of the outer jet (and thus the whole outflow) increases due to the diffusion of dense material from the outer part to the static medium. This is more of interchange type, but does not affect significantly the integrity of the outer jet significantly. The trends (Lorentz factor and effective radius over time for different $\sigma$) are presented in Fig. 14.

The increase in the effective radius of the whole jet is smaller compared to Cases 1–5. Even for $\sigma = 0.001$, the increase of the effective jet radius is ~30 per cent compared to ~40 per cent of Case 1. For a more detailed comparison, we refer again to Table 1.

## 7 3D SIMULATIONS

In this section, we examine some of the ‘extreme’ scenarios from the previous sections, namely the 3D equivalent of Cases 1 ($\sigma = 0.001$ and $\gamma_{in} \simeq 30$) and 7 ($\sigma = 0.1$ and $\gamma_{in} \simeq 10$), which are the most unstable and the most stable cases, respectively. In both cases, we refer again to rotating jets, with $v_{\phi} = 0.01$ at $R = R_{in}$.

We remind that for the 3D simulations, we use a uniform, Cartesian grid with a resolution of 256³. While this is lower than the effective resolution used in the 2.5D cases, it is still quite sufficient to capture small-scale phenomena. The perturbation in the radial velocity, in both cases, satisfies $\nabla \cdot \mathbf{v} = 0$.

### 7.1 Case 1b: $\sigma = 0.001$, $\gamma_{in} = 30$

We examine the most unstable case of Section 4, with $\sigma = 0.001$. The aim is to check if the enhanced Rayleigh-Taylor-type instability is still present (e.g. at $z = 0.5$ pc) and if other types of instabilities are dominant.

Initially, on the x, y plane (see Fig. 15), the evolution of the jet displays no major differences compared to the previous cases. Regarding the inner jet, since again no specific modes are excited, we notice the formation of four ‘arms’, a feature which persist up to ~ one rotation. This is in agreement with the behaviour of all the low-$\sigma$ 2.5D set-ups. Along the z axis, we notice the formation

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Proper density (log) for Case 8 (maximum $\sigma = 0.001$). In this case, the initial value of the rotation of the jet is close to zero. Snapshots after 0.5, one, two and three rotations of the inner jet. The minimal centrifugal force leads to a stable final state compared to Case 1.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Proper density (log) for Case 9 (maximum $\sigma = 0.1$). In this case, the initial value of the rotation of the jet is close to zero. Snapshots after 0.5 and three rotations of the inner jet. The increased magnetic tension further stabilizes the outflow, compared with Case 8.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Average Lorentz factor of the inner jet (top) and effective radius of the outer jet (bottom) after three rotations assuming a realistic Lorentz factor of 10 for the inner jet.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 16.** Left-hand column: Proper density (log) for Cases 1b (top) and 7b (bottom), with maximum $\sigma = 0.001$ and $\sigma = 0.1$, respectively. Right-hand column: same for the Lorentz factor. Snapshot after 0.5 rotations of the inner jet, view of the x, y plane, at $z = 0.5$ pc.</figcaption></figure>

of Kelvin–Helmholtz instabilities around ~0.3 rotations, mainly between the outer jet and the external medium.

A three-slice representation of the computational domain is given in the left-hand column of Fig. 16. A view of the x, y plane at $z = 0.5$ pc is given in the top left-hand (density) and top right-hand (Lorentz factor) images in Fig. 15. We highlight the rapid expansion of the outer jet due to the Kelvin–Helmholtz effects.

Overall, the set-up is more unstable than the 2.5D case, in terms of expansion of the outer jet mainly due to the development of Kelvin–Helmholtz instabilities along the z-axis (but we note that causality can still play a role in decelerating the jet). Up to half rotation time, the two components can still be distinguished in terms of density and Lorentz factor, a trend which persists roughly up to a full rotation of the inner jet. Up to one rotation, there is no significant scaling of the density or the Lorentz factor with z, for the inner jet and a part of the outer jet up to ~0.11 pc. This is shown in the bottom image of Fig. 16, where we plot the density variation with z in different regions (i.e. for different values of R). We use $R = R_{\text{in}}/2$ for the inner jet and $R = 0.25$ pc for the interface between the outer jet and the external medium (where we notice the formation of vortices due to Kelvin-Helmholtz instabilities).

## 7.2 Case 7b: $\sigma = 0.1$, $\gamma_{\text{in}} = 10$

We examine one of the most stable cases, the 3D variant of Case 7, with $\sigma = 0.1$ and $\gamma = 10$ for the inner jet. Compared to the most unstable cases, the magnetic tension is now significant and the inertia reduced, thus we expect no significant Rayleigh–Taylor-type instabilities. Kelvin–Helmholtz effects though can still be anticipated due to the density and velocity shear between the outer jet and the environment.

In the early stages of the simulation, up to ~0.3 rotations there is no significant change in the jet. After this time, Kelvin–Helmholtz instabilities begin to develop along the z axis, between the outer jet and the external medium. In this case, the instabilities are somewhat less prominent due to the strong toroidal field compared to the
<figure><img src="image_19.png" /><figcaption>**Figure 17.** Proper density for Cases 1b (left-hand column, maximum $\sigma = 0.001$, $\gamma = 30$) and 7b (right-hand column, maximum $\sigma = 0.1$, $\gamma = 10$). Snapshots after 0.3 and 0.5 rotations of the inner jet. The flow is moving upwards along the z axis. The bottom image shows the density variation along the z direction, for different radii, after 0.5 rotations (Case 1b). The dashed line refers to $R = R_{in}/2$ and the solid line to the interface between the outer jet and the surrounding medium. The periodic behaviour of the density distribution with z at the interface is an indication of Kelvin–Helmholtz instabilities. A similar behaviour is noticed in Case 7b.</figcaption></figure>

previous case. We emphasize that the velocity shear between these media is the same in every case, since we keep the Lorentz factor of the outer jet fixed to $\gamma = 3$.

A three-slice representation can be seen in the right-hand column of Fig. 16 and the x, y plane at $z = 0.5$ pc in Fig. 15. We notice that the strong tension prevents the formation of ‘arms’ between the inner and the outer jet and, at least in terms of Lorentz factor, the configuration is stable.

Overall, the setup is similar to the relevant 2.5D case, with the two jet components remaining separate, distinguishable in terms of density and Lorentz factor, and the jet collimated, at least up to one rotation. In terms of expansion of the outer jet, we notice an increase in the effective radius of the outer part, again mainly due to the development of the Kelvin–Helmholtz instability along the z axis. The outer jet has a reduced density, compared to the initial one, due to the mixing induced by Kelvin–Helmholtz.

## 8 SEMIANALYTICAL INTERPRETATION

We will examine the above cases following a semi-analytical approach, similar to Meliani & Keppens (2009). This interpretation focused on velocity perturbations in the radial direction only and excluded terms of $B_\phi$ in the momentum equation, e.g. magnetic tension. Since in our case the magnetic field includes a toroidal component, we want to examine the importance of the extra term(s) in the overall (in)stability of each different case.

The momentum equation, using the term $\frac{dV}{dt}$ is (Kalra & Gebretsadkan 2000)

$$ \gamma^2 \left[ \rho h + B^2 - (\mathbf{V} \cdot \mathbf{B})^2 \right] \frac{d\mathbf{V}}{dt} + \nabla \left( p + \frac{B^2}{2} \right) + \mathbf{V} \frac{\partial}{\partial t} \left( p + \frac{B^2}{2} \right) - \left\{ \left( \frac{\mathbf{B}}{\gamma} \cdot \nabla \right) \frac{\mathbf{B}}{\gamma} + \frac{\mathbf{B}}{\gamma} \left( \nabla \cdot \frac{\mathbf{B}}{\gamma} \right) + (\mathbf{V} \cdot \mathbf{B}) \left[ \frac{d\mathbf{B}}{dt} + \mathbf{B} \left( \nabla \cdot \mathbf{V} \right) + \left( \mathbf{B} \cdot \nabla \right) \mathbf{V} \right] + \mathbf{B} \frac{d}{dt} \left( \mathbf{V} \cdot \mathbf{B} \right) \right\} = 0, \quad (10) $$

whereas the HD limit can be found in (Bodo et al. 2004).

In Meliani & Keppens (2009), the approximate momentum equation, ignoring charge separation, was

$$ (\gamma^2 \rho h + B_z^2) \left[ \frac{\partial}{\partial t} + \mathbf{V} \cdot \nabla \right] \mathbf{V} + \nabla p_{\text{tot}} + \mathbf{V} \frac{\partial p_{\text{tot}}}{\partial t} \approx 0, \quad (11) $$

where the total pressure is $p_{\text{tot}} = p + \frac{B_z^2}{2}$. As in Meliani & Keppens (2009), we evaluate the assumption of the minimal charge separation contribution (lab frame) in the momentum equation by examining the maximum radial and toroidal velocity components. In every case, $v_\phi < 2 \times 10^{-2}$ and $v_r < 10^{-3}$ and thus the contribution of the $\rho_e E$ term is negligible.

Starting from equation (10), we extend this expression by including the toroidal component of the magnetic field, while we maintain the assumption that the velocity is mainly poloidal ($v_z \gg v_\phi$) and that $\gamma \approx \text{const.}$ initially in both components. The momentum equation will then be

$$ (\gamma^2 \rho h + B_z^2 + \gamma^2 B_\phi^2) \left[ \frac{\partial}{\partial t} + \mathbf{V} \cdot \nabla \right] \mathbf{V} + \nabla p_{\text{tot}} + \mathbf{V} \frac{\partial p_{\text{tot}}}{\partial t} - \frac{1}{\gamma^2} (\mathbf{B} \cdot \nabla) \mathbf{B} - 2(\mathbf{V} \cdot \mathbf{B})(\mathbf{B} \cdot \nabla) \mathbf{V} - \mathbf{B} \frac{d(\mathbf{V} \cdot \mathbf{B})}{dt} \approx 0. \quad (12) $$

The importance of the toroidal magnetic field component is expressed initially in the first parenthesis and via the magnetic tension term of the Lorentz force, i.e. $(\mathbf{B} \cdot \nabla)\mathbf{B}$. In the first case, the ratio $\frac{\gamma^2 B_\phi^2}{\gamma^2 \rho h}$ is of the order $\frac{\gamma^2 \sigma}{h}$, which is less than 0.4 for $\sigma$ up to 0.1.

If we further approximate equation (12) neglecting the time derivative of $p_{\text{tot}}$ and assume that the velocity perturbation does not depend on z or $\phi$, we obtain

$$ (\gamma^2 \rho h + B_z^2 + \gamma^2 B_\phi^2) \left[ \frac{\partial v_R}{\partial t} - \frac{V_\phi^2}{R} \right] + \frac{\partial p_{\text{tot}}}{\partial R} + \frac{B_\phi^2}{\gamma^2 R} \approx 0. \quad (13) $$

From this expression, we notice that the jet inertia increases as $\gamma^2$. Moreover, the increase of $B_\phi$ in the outer jet further increases the jet stability.
We expect the magnetic tension to balance the centrifugal force for sufficiently high values of magnetization. Indeed, the ratio of the centrifugal force to magnetic tension is $\frac{h}{a}\gamma^2V_\phi^2 \sim \frac{h}{a}\sigma$. The inner part of the jet has a significantly high enthalpy, whereas the outer part is cold and the tension is dominant for higher values of $\sigma$. Moreover, simulations with lower values of $\gamma$ for the inner jet are also more stable, since the ratio is proportional to $\gamma^2$.

## 9 CONCLUSIONS

We performed 2.5D and 3D simulations of relativistic, two-component jets, assuming differential rotation and a toroidal magnetic field component of different magnitude in each case. The chosen parameters (e.g. the Lorentz factor $\gamma$) correspond to AGN jets and each (2.5D) simulation is restricted to a plane $(x, y)$ perpendicular to the jet axis, far from the central source. In every case, we assumed that the jet has already been accelerated and collimated.

We keep the same density ratio and toroidal velocity profile (the latter with different coefficients) in each case, as described in Section 3. Depending on the choice of the maximum magnetization, we result in a different final state, in terms of mixing between the components, the average Lorentz factor and the effective radius of the jet.

The key difference in this work, in comparison with previous studies (Meliani & Keppens 2007, 2009), is the addition of a toroidal component in the magnetic field, in agreement with observations of real astrophysical jets. This extra component results in an additional term in the momentum equation, which represents the magnetic tension. Although we also gradually (from lower to higher $\sigma$) arbitrarily decrease the contribution of the $B_{\perp}$ component of the magnetic field in the magnetic pressure, and thus the effective inertia of the inner jet, the big picture depends mainly on the choice of the magnetization. We also note that the assumption of a continuous toroidal velocity field leads to a stronger centrifugal force in the outer jet, in contrast with Meliani & Keppens (2009), where the toroidal velocity component was discontinuous.

In all cases of Section 4, a shear region between the inner and the outer jet is formed, with mixing between the components also taking place. The shear region and the mixing are different in each setup, with cases of very low magnetization (below $\sigma = 0.01$) resulting in stronger mixing. Case 1 ($\sigma = 0.001$) is the most unstable, with the average Lorentz factor dropping to $\sim 11$ and the jet expanding to an effective radius of $\simeq 0.09$ pc. The overall image of the unstable cases, resembles that of Meliani & Keppens (2009), where the jets were decelerated and expanding.

For higher values of magnetization, we can still distinguish an inner region of high Lorentz factor, although the maximum value is reduced when compared to the initial one. This effect is prominent in Case 5 ($\sigma = 0.1$). Even though mixing is still present, its effect is not that strong, but still capable to decelerate the inner jet to $\gamma \sim 20$. In the end, we find that a toroidal magnetic field, of sufficiently high $\sigma$, damps the Rayleigh-Taylor-type instability and makes the interface between the inner and the outer jet stable.

Furthermore, we examine cases with slower inner jets to achieve a more realistic value of Lorentz factor $\gamma \sim 10$, which tend to be more stable even at low magnetizations. This is attributed to the lower effective inertia of the inner jet. If we keep the rotation of the jet to a minimum, the centrifugal force is negligible, unable to trigger any Rayleigh-Taylor-type instabilities and thus the outflow still remains collimated after three rotations. These cases of course can still be

unstable to other types of instabilities, e.g. Kelvin-Helmholtz or kink, which require a 3D study and a relaxation of the translational symmetry along the $z$ axis.

As regards the 3D simulations, at least for the asymptotic cases we examined, jets with weak toroidal magnetic field tend to remain unstable to the Rayleigh-Taylor-type instability. Increasing the magnetization, we achieve stability but we note that both cases are still unstable to Kelvin-Helmholtz instabilities, which is dominant now up to one rotation. In this study, we are modelling but a small region of the jet in 3D, with the upper $z$ boundary being equal to 10 jet radii. Using periodic boundary conditions though in the $z$ direction, we can represent a larger region (of sub-parsec to parsec scale) of the jet.

Overall, our results are complementary to the findings of Meliani & Keppens (2009), where the proposed models were used as an analogy to explain the FRI/FRII classification. Furthermore, the (in)stability of such an outflow has impact on its structure (e.g. in terms of velocity, mixing etc.), which in turn can affect the emission we observe from these objects. This work will be continued in a follow-up study, synthetic observations will be created based on our 2.5D and 3D results. With the translational symmetry relaxed, it is also possible to see if different kinds of instabilities (e.g. kink) develop. Furthermore, we plan to simulate different kinds of outflows, e.g. GRB jets and create synthetic radiation maps, where polarized emission can also be considered.

## ACKNOWLEDGEMENTS

This research was supported by projects GOA/2015-014 (2014-2018 KU Leuven) and the Interuniversity Attraction Poles Programme by the Belgian Science Policy Office (IAP P7/08 CHARM). The computational resources and services used in this work were provided by the Vlaams Supercomputer Centrum (Flemish Supercomputer Center), funded by the Research Foundation Flinders (FWO) and the Flemish Government - department EWI. Visualization was performed using PARAVIEW (more information on www.paraview.org) and GNUPLOT (http://www.gnuplot.info/index.html). DM thanks Bart Ripperda, Norbert Magyar, Fabio Bacchini and Jo Raes for their comments on various stages of this work and Chun Xia and Oliver Porth for their suggestions regarding the use of the code. Finally, we thank the anonymous referee for the useful comments on the manuscript.

## REFERENCES

* Baty H., Keppens R., 2002, ApJ, 580, 800

* Biretta J. A., Junor W., Livio M., 2002, New Astron Rev., 46, 239

* Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883

* Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433

* Boccardi B., Krichbaum T. P., Bach U., Mertens F., Ros E., Alef W., Zensus J. A., 2016, A&A, 585, A33

* Bodo G., Mignone A., Rosner R., 2004, Phys. Rev. E, 70, 036304

* Bodo G., Mamatsashvili G., Rossi P., Mignone A., 2013, MNRAS, 434, 3030

* Bogovalov S. V., 1995, Astron. Lett., 21, 565

* Bogovalov S., Tsinganos K., 2001, MNRAS, 325, 249

* Čada M., Torrilhon M., 2009, J. Comput. Phys., 228, 4118

* Casse F., Keppens R., 2004, ApJ, 601, 90

* Fanaroff B. L., Riley J. M., 1974, MNRAS, 167, 31P

* Gabuzda D. C., Reichstein A. R., O'Neill E. L., 2014, MNRAS, 444, 172

* Gabuzda D. C., Knuettel S., Reardon B., 2015, MNRAS, 450, 2441

* Ghisellini G., Tavecchio F., Chiaberge M., 2005, A&A, 432, 401
Giroletti M. et al., 2004, ApJ, 600, 127

Globus N., Levinson A., 2016, MNRAS, 461, 2605

Heyvaerts J., Norman C., 2003, ApJ, 596, 1270

Kalra G. L., Gebretsadkan W. B., 2000, Phys. Plasmas, 7, 3531

Keppens R., Meliani Z., van Marle A. J., Delmont P., Vlasis A., van der Holst B., 2012, J. Comput. Phys., 231, 718

Komissarov S. S., Barkov M. V., Vlahakis N., Königl A., 2007, MNRAS, 380, 51

Matsakos T., Tsinganos K., Vlahakis N., Massaglia S., Mignone A., Trussoni E., 2008, A&A, 477, 521

Matsakos T., Massaglia S., Trussoni E., Tsinganos K., Vlahakis N., Sauty C., Mignone A., 2009, A&A, 502, 217

Matsumoto J., Masada Y., 2013, ApJ, 772, L1

Meliani Z., Keppens R., 2007, A&A, 475, 785

Meliani Z., Keppens R., 2009, ApJ, 705, 1594

Meliani Z., Keppens R., Giacomazzo B., 2008, A&A, 491, 321

Mertens F., Lobanov A. P., 2016, A&A, 587, A52

Mertens F., Lobanov A. P., Walker R. C., Hardee P. E., 2016, A&A, 595, A54

Mizuno Y., Lyubarsky Y., Nishikawa K.-I., Hardee P. E., 2012, ApJ, 757, 16

Mizuno Y., Hardee P. E., Nishikawa K.-I., 2014, ApJ, 784, 167

Mobarry C. M., Lovelace R. V. E., 1986, ApJ, 309, 455

Piner B. G., Unwin S. C., Wehrle A. E., Zook A. C., Urry C. M., Gilmore D. M., 2003, ApJ, 588, 716

Porth O., Komissarov S. S., 2015, MNRAS, 452, 1089

Porth O., Xia C., Hendrix T., Moschou S. P., Keppens R., 2014, ApJS, 214, 4

Rossi P., Mignone A., Bodo G., Massaglia S., Ferrari A., 2008, A&A, 488, 795

Singh C. B., Mizuno Y., de Gouveia Dal Pino E. M., 2016, ApJ, 824, 48

Sądowski A., Narayan R., Tchekhovskoy A., Zhu Y., 2013, MNRAS, 429, 3533

Tchekhovskoy A., Narayan R., McKinney J. C., 2011, MNRAS, 418, L79

Teşileanu O. et al., 2014, A&A, 562, A117

Vlahakis N., 2004a, Ap&SS, 293, 67

Vlahakis N., 2004b, ApJ, 600, 324

Woitas J., Ray T. P., Bacciotti F., Davis C. J., Eislöffel J., 2002, ApJ, 580, 336

## APPENDIX: INTEGRATION CONSTANTS IN PRESSURE EQUATION

The full expressions for $\zeta$ in equation (9) are (Meliani & Keppens 2009):

$$ \zeta_{\text{in}} = p_o + \frac{\gamma_{\text{in}} - 1}{\gamma_{\text{in}}} \rho_{\text{in}} + \frac{(\Gamma - 1)(\tilde{a} + 2)}{2\tilde{a}(\alpha(\Gamma - 1) + \Gamma)} \times (-b_{\phi}^2 + (b_{\phi}V_z - v_{\phi}B_z)^2) \quad (\text{A1}) $$

for the inner jet and

$$ \zeta_{\text{out}} = \left[ p_{\text{tot}}(R_{\text{in}}) + \frac{\Gamma_{\text{out}} - 1}{\Gamma_{\text{out}}} \rho_{\text{out}} - \frac{B_{z,\text{out}}^2}{2} \right] (1 - \tilde{a}_{\text{out}})^{\frac{\Gamma_{\text{out}}}{\alpha_{\text{out}}(\Gamma_{\text{out}} - 1)}} + [-b_{\phi,\text{out}}^2 + (b_{\phi,\text{out}}V_{z,\text{out}} - v_{\phi,\text{out}}B_{z,\text{out}})^2] \\ \times \left(1 - \tilde{a}_{\text{out}}\right)^{\frac{\Gamma_{\text{out}}}{\alpha_{\text{out}}(\Gamma_{\text{out}} - 1)}} \\ \times \left(\frac{1}{2} + \frac{(1 - \tilde{a}_{\text{out}})(\Gamma_{\text{out}} - 1)(\alpha_{\text{out}} + 2)}{2\tilde{a}_{\text{out}}[\alpha_{\text{out}}(\Gamma_{\text{out}} - 1) + \Gamma_{\text{out}}]}\right) \quad (\text{A2}) $$

for the outer jet. We remind that $\gamma_{\text{in}} = 4/3$ and $\Gamma_{\text{out}} = 5/3$.

This paper has been typeset from a TeX/LA TeX file prepared by the author.