© The Authors 2025

1. Introduction

Relativistic jets from AGNs are among the most powerful and persistent astrophysical outflows (for a review, see Blandford et al. 2019). These jets extend from subparsec scales near the supermassive black hole to megaparsec distances, emitting nonthermal radiation across the electromagnetic spectrum and potentially producing high-energy neutrinos and ultra-high-energy cosmic rays (e.g., O’Sullivan et al. 2009; Petropoulou et al. 2020; Boula & Mastichiadis 2022). Despite their remarkable large-scale stability (Willis et al. 1974), jets are susceptible to a broad range of hydrodynamic and magnetohydrodynamic instabilities, including Kelvin–Helmholtz, Rayleigh–Taylor, centrifugal, current-driven and pressure-driven modes (e.g., Bodo et al. 1989, 2004, 2013; Bodo et al. 2019; Millas et al. 2017; Das & Begelman 2019; Wang et al. 2023; Bodo et al. 2022; Musso et al. 2024). High-resolution observations of nearby radio galaxies, such as M87 and BL Lac, reveal jets that undergo rapid lateral expansion and acceleration near the core, followed by gradual collimation at larger distances (Cohen et al. 2014; Walker et al. 2018). This structural evolution is often associated with recollimation shocks, which arise from pressure mismatches between the jet and the surrounding medium (Komissarov & Falle 1997; Bodo & Tavecchio 2018; Gourgouliatos & Komissarov 2018; Costa et al. 2024, 2025). These shocks not only influence the jet’s morphology but are also considered prime sites for localized particle acceleration via oblique shocks and stochastic processes in the downstream flow (e.g., Sciaccaluga et al. 2025). Observationally, recollimation shocks may also be linked to high-energy variability; for instance, Hervet et al. (2019) proposed that multiple stationary shocks in Mrk 421 could explain complex X-ray flaring patterns, suggesting that recollimation structures actively modulate nonthermal emission. Collectively, these findings underscore the importance of parsec-scale recollimation zones for both the dynamical and radiative properties of AGNs jets. While early models and axisymmetric simulations predicted a stable sequence of recollimation and reflection shocks (Komissarov & Falle 1997; Mizuno et al. 2015; Porth & Komissarov 2015; Fromm et al. 2016), 3D studies have shown that these structures become unstable beyond the first shock. Simulations reveal the rapid growth of instabilities, the onset of turbulence, and the suppression of multiple shock structures that are commonly observed in 2D (Matsumoto & Masada 2013; Mizuno et al. 2009; Gourgouliatos & Komissarov 2018; Matsumoto et al. 2021). A key mechanism in this context is the CFI, triggered by the curvature of streamlines in the recollimation flow. The CFI is intrinsically 3D and can develop even in the absence of a jet–ambient density contrast (Gourgouliatos & Komissarov 2018). Its nonlinear evolution and dependence on jet parameters, however, remain incompletely understood (Costa et al. 2025). Magnetic fields further complicate the stability picture (Gourgouliatos & Komissarov 2018; Komissarov et al. 2019; Gottlieb et al. 2020; Matsumoto et al. 2021; Hu et al. 2025). In magnetized jets, current and pressure-driven instabilities can emerge and compete with hydrodynamic modes. Recent studies indicate that the kink instability may dominate over conical shock formation in highly magnetized flows, significantly altering jet propagation and energy dissipation patterns (Barniol Duran et al. 2017). Linear analyses have provided insights into the growth of these modes under various configurations (Bodo et al. 2013, 2019; Kim et al. 2018; Vlahakis 2023, 2024; Musso et al. 2024). Overall, the interplay between jet magnetization, velocity shear, density, and the external pressure gradient critically determines jet stability. In this work we investigate the onset of instabilities in magnetized recollimating flows, under conditions representative of low-power AGNs jets, using 3D relativistic hydrodynamic simulations performed with the PLUTO code (Mignone et al. 2007) complemented by a linear stability analysis of the CFI. While previous studies have focused mainly on 2D setups or isolated instabilities, our approach captures the intrinsic 3D character of the CFI, its interaction with magnetic fields, and the transition to turbulence downstream of the first recollimation shock. The linear analysis quantifies the growth rates and parameter dependence of the instability, while the 3D simulations reveal its nonlinear evolution, energy dissipation, and impact on jet structure. By combining these two approaches, we provide new insights into jet stability, the role of magnetic stresses, and the potential connection to observable signatures such as stationary features and high-energy variability.

The paper is organized as follows. In Sect. 2, we present the model description and the numerical setup. In Sect. 3.1, we show the results of our 2D simulations. In Sect. 3.2 we present the linear analysis of the relativistic CFI. In Sect. 3.4 we present the results of the 3D simulations. In Sect. 4, we discuss the implications of the results in the context of jet stability and the development of the CFI and we summarize our conclusions and outline prospects for future work.

2. Problem description

Observations suggest the presence of quasi-stationary features in AGNs jets on parsec scales, likely associated with recollimation shocks (see for example Paraschos 2025). These regions have been argued to be the expected location of the high-energy emission zone in blazars. At these distances, variability is typically moderate, occurring on timescales of about a day, which suggests that the emission is produced by a relatively stable dissipation mechanism (Tavecchio et al. 2010). This challenges scenarios involving rapidly propagating shocks and points instead to localized, possibly stationary structures as acceleration sites.

In this work following the approach adopted in Costa et al. (2024) for a purely hydrodynamical jet, we investigate the recollimation of a relativistic jet at parsec-scale distances from the central engine (the influence of general relativity is negligible in this context). We consider a stage in which the bow shock at the jet head has dissipated and the jet propagates without the protection of a surrounding cocoon. Instead, it interacts directly with the ambient medium, which we interpret as a residual wind or warm interstellar gas still influenced by the gravitational potential of the central source. This environment exerts a confining pressure, for which we assume a decreasing power-law profile, with index η = 0.5. A free expanding jet is eventually reconfined by the external pressure, with the formation of a chain of recollimation shocks (Komissarov & Falle 1997; Porth & Komissarov 2015).

Our methodology is similar in spirit to the works of Gourgouliatos & Komissarov (2018) and Matsumoto et al. (2021), who also studied the stability of recollimated jets using a two-step simulation process. In particular, we followed the two-step approach presented in Costa et al. (2024, 2025). We note that the equilibrium configuration adopted as the initial state represents an idealized setup that may not be fully realized in nature. This choice, however, allows us to isolate and analyze the growth of specific instabilities under controlled conditions, following the approach of previous works. First, we performed 2D axisymmetric simulations to identify steady-state solutions featuring a clear recollimation structure. Specifically, a relativistic conical jet with an opening angle θ_j is launched into the computational domain at a distance z₀ from the cone’s apex, embedded within a confining ambient medium. We extended the setup by incorporating magnetic fields into the simulations. Then, we extended the analysis to 3D, relaxing the symmetry constraints to assess the stability of these solutions and capture the development of instabilities. This strategy allows us to isolate instability growth from jet propagation effects. We note that different procedures for constructing the steady initial configuration, such as starting from an analytical equilibrium or from a numerically relaxed axisymmetric state, may lead to small variations in the early 3D evolution. However, these differences do not affect the qualitative behavior of the instability during its nonlinear development.

The system of equations in conservation form reads

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}\left(\begin{array}{c}D\\ \mathbf{m}\\ E_t\\ \mathbf{B}\end{array}\right)+\mathrm{\nabla }·\left(\begin{array}{c}D\mathbf{u}\\ w_t{\mathrm{\Gamma }}^2\mathbf{u}u-\mathbf{b}b+\mathbf{I}{p}_t\\ \mathbf{m}\\ \mathbf{u}B-\mathbf{B}u\end{array}\right)=\left(\begin{array}{c}0\\ {\mathbf{f}}_g\\ \mathbf{u}·f_g\\ 0\end{array}\right).\end{array}$$(1)

We adopt units in which the speed of light c = 1, and the magnetic field B is rescaled to absorb the factor $\sqrt{4\pi }$. The set of primitive variables is given by the density in the rest frame ρ, the thermal pressure p, the three-velocity u in the laboratory frame and the three-vector magnetic field B in the laboratory frame. Here, Γ denotes the Lorentz factor, f_g is an external force density in the lab frame, and I represents the 3 × 3 identity tensor. The conserved variables include the magnetic field B, as well as the lab-frame mass, momentum, and energy densities, defined as:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill D& =\mathrm{\Gamma }\rho ,\hfill \\ \hfill \mathbf{m}& =w_t{\mathrm{\Gamma }}^2\mathbf{u}-b_0\mathbf{b},\hfill \\ \hfill E_t& =w_t{\mathrm{\Gamma }}^2-b_0{b}_0-p_t,\hfill \end{array}\text{with}\begin{array}{cc}\hfill b_0& =\mathrm{\Gamma }\mathbf{u}·\mathbf{B},\hfill \\ \hfill \mathbf{b}& =\frac{\mathbf{B}}{\mathrm{\Gamma }}+\mathrm{\Gamma }(\mathbf{u}·\mathbf{B})\mathbf{u},\hfill \\ \hfill w_t& =\rho h+\frac{B^2}{{\mathrm{\Gamma }}^2}+{(\mathbf{u}·\mathbf{B})}^2,\hfill \\ \hfill p_t& =p+\frac{{\mathbf{B}}^{\mathbf{2}}}{2{\mathrm{\Gamma }}^2}+\frac{{(\mathbf{u}·\mathbf{B})}^2}{2}.\hfill \end{array}\end{array}$$(2)

We close the set of equations with the Taub-Matthews equation of state, which approximates the Synge EoS of a single species relativistic perfect fluid (Mignone et al. 2005):

$$\begin{array}{c}\hfill h=\frac{5}{2}\mathcal{T}+\sqrt{\frac{9}{4}{\mathcal{T}}^2+1},\end{array}$$(3)

where h is the specific enthalpy and 𝒯 = p/ρ is the temperature.

As explained above at t = 0 we have a conical jet, in which the velocity is constant and the density and pressure profiles are, following Komissarov & Falle (1997):

$$\begin{array}{cc}& {\rho }_j(r,z,t=0)={\rho }_j(0,z_0,0){\left(\frac{R}{R_0}\right)}^{-2},\hfill \end{array}$$(4)

$$\begin{array}{cc}& P_j(r,z,t=0)=P_j(0,z_0,0){\left(\frac{R}{R_0}\right)}^{-2\gamma },\hfill \end{array}$$(5)

where r is the cylindrical radius (measured from the axis of the cone), $R=\sqrt{r^2+z^2}$ is the spherical radius and R₀ = z₀. Furthermore, the transverse transition of all the variables from their values in the jet to the values in the external medium are smoothed to avoid numerical noise at the contact discontinuity, as better specified in Appendix A. The external medium is characterized by density and pressure that decay with distance from the central object as a power law. The external density profile along z is ρ_ext(z) = ρ_ext, 0(z/z₀)^−η and the external pressure profile is P_ext(z) = P_ext, 0(z/z₀)^−η.

The magnetic field is helical and its components, in the laboratory frame, are (see Sciaccaluga et al. 2025, for details):

$$\begin{array}{cc}& B_r=\frac{\alpha B_0}{R^2}{\text{e}}^{-{\mathcal{X}}^2}\frac{r}{R},\hfill \end{array}$$(6)

$$\begin{array}{cc}& B_z=\frac{\alpha B_0}{R^2}{\text{e}}^{-{\mathcal{X}}^2}\frac{z}{R},\hfill \end{array}$$(7)

$$\begin{array}{cc}& B_{\varphi }=\mathrm{\Gamma }\frac{B_0}{R}\sqrt{{\text{e}}^{-2{\mathcal{X}}^2}-\frac{{\psi }_{\chi }}{2{sin}^2\theta }[{\psi }_{\chi }-{\psi }_{\chi }{\text{e}}^{-2{\mathcal{X}}^2}+\sqrt{2\pi }\mathrm{erf}\left(\sqrt{2}\mathcal{X}\right)]},\hfill \end{array}$$(8)

where Γ is the Lorentz factor of the flow, 𝒳 = (cos θ−1)/ψ_χ and θ = cos⁻¹(z/R). The parameter α is related to the pitch and the amplitude B₀ is defined in the rest frame, it defines the field strength and is related to the rest frame hot magnetization, at injection, σ by the relation

$$\begin{array}{c}\hfill \sigma =\frac{B_0^2}{\rho (z_0)h}.\end{array}$$(9)

Finally, ψ_χ is the scale length for the decay of the poloidal field, which is assumed to have a Gaussian profile. We set α = 1 meaning that the maximum strengths of the poloidal and toroidal components are comparable in the fluid rest frame, while in the lab frame the toroidal component dominates due to the Lorentz factor.

In our simulations, Table 1, we consider a jet with an initial opening angle θ_j = 0.1 and Lorentz factor Γ = 10. We set ψ_χ = 10⁻²θ_j, to make the field decay to 0 outside the jet. The density ratio of the jet to the ambient medium is ν = ρ_j(r = 0, z = z₀)/ρ_ext, 0 = 10⁻⁵, 10⁻⁴. The pressure ratio is P_ratio = P_j(r = 0, z = z₀))/P_ext, 0 = 10⁻³, 10⁰, 10¹ and the magnetization ranges from 0 (unmagnetized) to 1.

The total jet power L_jet is the sum of the kinetic (matter) power and the magnetic (Poynting) power L_jet = L_HD + L_B where: L_HD = πz₀²θ_j²ρ_jh_jΓ_j²v_j and $L_B=\pi z_0^2{\theta }_j^2\frac{B_o^2}{4\pi }{\mathrm{\Gamma }}_j^2$. Putting it together, we get:

$$\begin{array}{c}\hfill L_{\text{jet}}=\pi z_0^2{\theta }_j^2{\mathrm{\Gamma }}_j^2({\rho }_j{h}_j{v}_j+B_0^2),\end{array}$$(10)

which, for the given parameters, corresponds to low power jets (with a value of ∼ 10⁴³ erg/sec, for z₀ = 0.1 pc and ρ_ext, 0 = 10⁵ m_p cm⁻³).

The numerical simulations are carried out by solving the equations of relativistic ideal magnetohydrodynamics (MHD) first in 2D and then in 3D, by using the ideal RMHD module of the PLUTO code (Mignone et al. 2007). The equations were evolved with second-order Runge-Kutta time stepping, linear reconstruction, and the HLLD Riemann solver (Mignone et al. 2009). Beyond these methods, magnetic field evolution required additional care. Constrained transport (Balsara & Spicer 1999; Londrillo & del Zanna 2004) was used to maintain the condition ∇ ⋅ B = 0.

The setup of each 2D simulation uses cylindrical coordinates (r, z), with a domain [0,L_r,2D] × [1,L_z,2D]. All lengths are expressed in units of z₀, defined as the distance from the cone (Appendix B for details). At the final time, the 2D simulations provide a nearly steady state which is used as the initial state for the 3D simulations. 3D simulations were carried out on a Cartesian domain with coordinates in the range x ∈ [−L_x,3D/2,L_x,3D/2], y ∈ [−L_y,3D/2,L_y,3D/2], and z ∈ [1,L_z,3D] (lengths are expressed in units of the z₀). Here L_x = L_y, while we chose z as the jet propagation direction. We employed a grid of N_x × N_y × N_z zones with a uniform grid spacing along the z (axial) direction (Appendix B for details).

3. Results

To build a coherent picture of jet stability, we started from the simplest setting and add complexity step by step. Our analysis begins with 2D steady-state solutions, which show that even jets with the same Lorentz factor, opening angle, and magnetization can develop different internal structures depending on the density and pressure ratios. Therefore, they may also have different stability properties with respect to the CFI. We then derive a local onset criterion and perform a linear stability analysis to determine how the growth rate of the CFI depends on jet parameters. With this insight, we then select two contrasting 3D simulations–one prone to disruption and one remaining stable–allowing us to investigate how magnetization, curvature, and external confinement jointly govern jet stability.

3.1. Insights from the 2D simulations

In this section we present the main findings from our 2D relativistic MHD simulations, which are used to construct steady jet solutions under different physical conditions. By varying the density ratio, pressure, and magnetization, we identify equilibria that serve as a baseline for assessing stability in 3D and for comparison with linear analysis. At the same time, the 2D runs already provide valuable insight into how the CFI can be triggered or suppressed, highlighting parameter combinations that lead either to long-lived stable flows or to configurations prone to disruption. Although 2D simulations cannot capture the full spectrum of 3D modes, they reproduce many of the essential dynamical features of relativistic jets. These include the formation of recollimation shocks, oscillatory jet profiles, and the redistribution of magnetic and thermal pressure. In particular, the 2D solutions offer a diagnostic of where strong curvature and flow gradients develop – precisely the conditions that are expected to enhance the growth of centrifugal and related instabilities.

Figure 1 shows the normalized Lorentz factor profiles along the z-axis at x = 0, y = 0 for simulations A, B, C, and D. The normalization is taken at the injection point (0, 0). Despite identical initial Lorentz factors and jet opening angles within each set, the location of the first recollimation shock shifts noticeably between cases. This shift is solely due to changes in the density ratio, temperature and magnetization, indicating that these parameters influence the effective curvature radius and radial expansion rate of the jet.

Fig. 1. Profiles along the z-axis at position x = 0 of the normalized Lorentz factor for the 2D simulations presented in Table 1. The normalization is taken at the position (0, 0) for simulations A, B, C, and D with magnetization σ = 0, 0.1, and 1.

In cases A (cold, light case) the first recollimation point is the closest to the origin. We also observe that the Lorentz factor shows a steady decrease on top of the oscillations due to recollimation shocks. A similar evolution occurs in cases B: the first recollimation point is somewhat shifted while the larger inertia moderates the global deceleration relative to cases A. Cases C is warm (kT ∼ m_ec²) and injected close to pressure equilibrium with the ambient and we observe a slight radial expansion that leads to a conversion of thermal to kinetic energy and to an increase of the Lorentz factor.

Once the jet recollimates, however, a fraction of this kinetic energy is reconverted into internal energy and the Lorentz factor shows a global decrease that is even larger than in cases A. In cases D the jet is also warm but injected with a pronounced overpressure, that drives a significant radial expansion and efficient acceleration: the Lorentz factor can reach Γ ≃ 75 before the recollimation point. In this case the thermal pressure plays a crucial role in converting internal energy into bulk motion.

Across all cases, a common trend emerges: despite the different injection conditions, the flow is not able to reconvert all the thermal energy into kinetic form. Increasing the magnetization strengthens the role of magnetic tension, which provides a more effective confinement of the jet and shifts the recollimation shocks to smaller axial distances; this behaviour is evident in all the cases shown in Fig. 1. A backward displacement of the shock is visible in each case, although the effect is very small in cases A. Since it appears systematically, it is reasonable to attribute it to the same underlying mechanism, likely related to the tension of the toroidal magnetic component. The detailed dependence of the shock location on magnetization will be presented in a forthcoming paper (Boula et al., in prep.).

The internal structure of a representative case is shown in Fig. 2 for simulation A_0.1 with σ = 0.1. We adopt σ = 0.1, as this choice leads to strong recollimation shocks, unlike the case with σ = 1. The logarithmic maps of total, thermal, and magnetic pressure, along with the total magnetization, reveal how different pressure components dominate in different regions. Upstream of the recollimation shock, thermal pressure contributes to acceleration as the flow expand, while magnetic pressure becomes more significant downstream of the shock, particularly along the jet edges. At the recollimation shock, magnetic pressure enhancements appear at the boundary layers, plausibly related to compression and flow deflection. A contribution from magnetic tension cannot be excluded, and this aspect will also be investigated in more detail in the forthcoming paper (Boula et al., in prep.).

Fig. 2. Logarithmic maps of total pressure, thermal pressure, magnetic pressure, and total magnetization for the 2D simulation A_0.1 with magnetization σ = 0.1. The initial parameters are: density ratio ν = 10−5, Lorentz factor Γ = 10, jet opening angle θjet = 0.1, and pressure ratio Pratio = 10−3. The maps correspond to the final time step t = 3000 tcr in units of z0/c.

3.2. Centrifugal instability

This study focuses on the CFI, a local instability that develops in curved flows. The presence of a toroidal magnetic field may, however, lead to stabilization due to the competing effect of magnetic tension.

To quantify the local conditions for CFI, we performed an order-of-magnitude comparison between the restoring magnetic tension due to the toroidal component of the magnetic field and the effective centrifugal force due to the motion along the curved streamlines. In a cylindrical or weakly perturbed geometry, the magnetic tension per unit volume scales as B_ϕ²/R_j, while the inertial term associated with curved motion scales as (Γ²ρ_jh_j + B²)/R_c, where B is the total magnetic field, B_ϕ is its toroidal component, R_j is the local jet radius, and R_c the streamline curvature radius (Komissarov et al. 2019; Matsumoto et al. 2021). Equating these contributions yields:

$$\begin{array}{c}\hfill \frac{B_{\varphi }^2}{R_j}\sim \frac{({\mathrm{\Gamma }}^2\rho h+B^2)}{R_c}\to \frac{{\sigma }_{\mathrm{tor}}}{(\sigma +1)R_j}\sim \frac{{\mathrm{\Gamma }}^2}{R_c}\to \frac{{\sigma }_{\mathrm{tor}}}{{\mathrm{\Gamma }}^2(\sigma +1)}\sim \frac{R_j}{R_c},\end{array}$$(11)

where σ and σ_tor are respectively the local values of the full magnetization and of the magnetization related to the toroidal component of the field. This shows that σ_tor/Γ² serves as a direct measure of the competition between stabilizing magnetic tension and destabilizing curvature (the total magnetization σ is less than the unity, and we can be neglected its contribution in the denominator).

A complementary perspective is provided in Fig. 3, which shows streamlines–initialized from the same starting point–overlaid on the spatial distribution of σ_tor/Γ² for each simulation. Using identical seeds ensures a direct visual comparison of streamline curvature and magnetization across cases. However, for the quantitative analysis presented below, we varied the initial streamline position to ensure that the selected trajectories actually pass through the regions favorable to the development of the CFI. For clarity, the figure highlights a representative streamline that illustrates the systematic trend across simulations.

Fig. 3. Streamlines (initialized from the same starting point) overlaid on the spatial distribution of σtor/Γ2 for all simulation runs. This ratio serves as a local diagnostic for the onset of centrifugal instability, with low values marking regions where the jet is most susceptible to curvature-driven perturbations. In the widest jets (e.g., Case D_0.1), outer streamlines would likely reach regions satisfying the CFI criterion; here we show a representative streamline to illustrate the overall trend.

The quantity σ_tor/Γ² quantifies the relative importance of toroidal magnetic tension to inertia and serves as a local diagnostic for the onset of the CFI. Regions with large σ_tor/Γ² are more resistant to bending and tend to suppress CFI growth. In contrast, regions of strong curvature and velocity shear often correspond to lower values of σ_tor/Γ², which makes them more favorable for the onset of the instability. In several runs, such regions appear already upstream, close to the local minimum of the jet radius, where perturbations can be seeded. In fully 3D flows, these perturbations may be further amplified by Richtmyer–Meshkov–type interactions as the jet passes through the reflection shock close to the minimum radius (Matsumoto & Masada 2013).

In Table 2 we report, for each run, the normalized magnetization σ_tor/Γ², the maximum transverse displacement R_j, max, the local curvature radius R_c (computed via circle fitting using ten points centered on R_j, max), and their ratio R_j, max/R_c. The last column lists the dimensionless combination (σ_tor/Γ²)/(R_j, max/R_c), which we use as a proxy for the balance between magnetic tension (stabilizing) and streamline curvature (destabilizing).

All quantities are evaluated at the axial location of maximum jet displacement (marked by the vertical lines in Fig. 4). This choice ensures consistency across runs, and test calculations indicate that this region is the most favorable for the development of the instability. The magnetization σ_tor is measured at the jet boundary along this axial position, where curvature effects are strongest.

Fig. 4. Left: Streamlines for cases A_0.1–D_0.1 (as in Fig. 3), starting at x0 = 0.082, 0.098, 0.091, and 0.124 respectively, with curvature and the ratio Rj, max/Rc. Right: Profiles of σtor/Γ2 versus x at z = zxmax, with vertical lines marking xmax = Rj, max.

Table 2 indicates that all four cases are unstable, with the instability proxy reflecting the combined effects of σ_tor/Γ² and curvature radius. Figure 5 compares the radial profiles of σ_tor, Γ, ρ, and h at z = z_max for Cases A_0.1 and D_0.1. Both cases exhibit smooth, but Case D_0.1 displays a wider and more complex radial transition region outside the recollimation shock. This illustrates that even when the instability proxy is similar, the spatial structure of the jet can vary, which may influence the growth and localization of local instabilities.

Fig. 5. Comparison of radial profiles of toroidal magnetization σtor (blue), Lorentz factor Γ (orange), rest-mass density ρ (green), and enthalpy h/ρ (red) at z = zmax for Case A_0.1 (top) and Case D_0.1 (bottom). Case A_0.1 shows smooth, monotonic profiles, whereas Case D_0.1 exhibits sharp gradients and non-monotonic behavior, underscoring the complexity of its internal structure. The vertical dashed line marks the position of xmax determined from the streamline analysis. Profiles are plotted on a logarithmic scale.

According to the Komissarov et al. (2019) criterion, suppression of CFI modes with azimuthal wavenumber m ≥ 4 at the interface between reconfined jets and external gas occurs for σ_tor/Γ² ≲ θ₀²/16, where θ₀ is the initial half-opening angle of the jet. In our simulations θ₀²/16 ≃ 0.000625. Comparing this value with σ_tor/Γ² (see Table 2) all cases are stable according to the Komissarov et al. (2019) criterion. This parameterizes the effect of geometry on the instability by using θ₀, but our results highlight the importance of the local geometry (i.e., the local curvature and jet radius) when studying the stability. In turn, this is determined not only by θ₀ but also by the dynamics of the jet.

To further quantify these findings, we now turn to a linear stability analysis aimed at computing the CFI growth rate as a function of jet parameters. This approach allows us to connect the local diagnostics discussed above with the expected timescale for instability development.

3.3. Linear analysis

Similarly to Komissarov et al. (2019), we examine the stability of a magnetized rotating cylindrical shell with a purely axial magnetic field. In contrast to that work, we carry out a linear stability analysis to compute the growth rate of the centrifugal instability, following the methodology of Vlahakis (2023)¹.

We adopt a cylindrical coordinate system (R, z) centered on the shell axis, distinct from the jet frame used elsewhere in this paper (Fig. 6 for details). The shell has outer radius R_c, width ΔR, and rotates with a constant azimuthal velocity V_ϕ, corresponding to a Lorentz factor Γ. Its specific enthalpy h and magnetization σ are taken to be uniform. The magnetic field is expressed as $\mathbf{b}=\sqrt{\sigma h\rho }\widehat{z}$, and the thermal pressure as P = ρ𝒯. Radial force balance then yields a power-law density profile ρ(R) = ρ_j(R/R_c)^a, with exponent a = (Γ² − 1)(1 + σ)/(𝒯/h + σ/2) We fix the shell width as ΔR = R_c/a.

Fig. 6. Cartoon illustrating the configurations used in the linear analysis, adapted from Komissarov et al. (2019). Top panel: Relativistic jet with a purely azimuthal magnetic field undergoes recollimation due to the confining action of an external pressure. The white inner region represents the freely expanding, unshocked flow, while the shaded area corresponds to the shocked outer layer. The interface between the two is the conical reconfinement shock driven into the jet. Bottom panel: < agnetised, rotating cylindrical shell with a purely axial magnetic field, representing the idealised configuration analysed in this work.

The shell is embedded in an external medium characterized by constant density ρ_ext, specific enthalpy h_ext, and a uniform axial magnetic field with strength $\sqrt{{\sigma }_{\text{ext}}{h}_{\text{ext}}{\rho }_{\text{ext}}}$. Pressure equilibrium at the shell boundary gives the density ratio: ρ_ext/ρ_Rc = (𝒯_Rc + σ_Rch_Rc/2)/(𝒯_ext + σ_exth_ext/2). This ratio, and in particular the quantity h_e ≡ h_ext/h, controls the density contrast between the shell and its surroundings. For minor deviations from h_e ≈ 1, the environment has little effect; more substantial effects appear only when h_e ≫ 1.

Since the CFI eigenfunction is localized near the shell’s outer edge, we model the inner boundary as a rigid wall. Perturbations of the form f(R)e^{i(kz − ωt)} are introduced in the linearized equations, which are then solved across the shell and the external region using appropriate matching conditions, as described in Vlahakis (2023). The resulting eigenvalue problem yields the (purely imaginary) instability growth rate ℑω as a function of the wavenumber k.

The results are presented in Fig. 7. The top panel displays the growth rate for various values of Γ and σ, with the density contrast held constant (h_e ≈ 1). As shown, increasing the Lorentz factor Γ or decreasing the magnetization σ significantly enhances the instability, as both the peak growth rate and the unstable wavenumber range increase. This behavior is consistent with the qualitative criterion proposed by Komissarov et al. (2019), in which the growth of the CFI is linked to the ratio (σ_tor/σ + 1)/Γ² ∼ σ_tor/Γ², since σ < 1.

Fig. 7. Results from the linear analysis of the CFI. Top panel: Growth rate as a function of kRc, for varying Lorentz factors and magnetizations, assuming fixed external conditions (he ≈ 1). Bottom panel: Same as top, but for fixed magnetization σ = 0.1, and varying external enthalpy ratio he and Lorentz factor Γ.

The bottom panel illustrates the effect of the external medium by fixing σ = 0.1 and varying both Γ and the enthalpy ratio h_e. A comparison between h_e ≈ 1 and h_e = 1.5 shows that the environment has a negligible impact unless h_e exceeds unity. For example, the growth rate curves for h_e ≈ 1 (h_e − 1 = 10⁻⁴) are indistinguishable, while the case h_e = 1.5 exhibits a marked increase in both growth rate and unstable range – especially for higher Lorentz factors. The choice of h_e = 1.5 reflects conditions in the simulations, where the external enthalpy is measured to reach values around 1–1.5; in practice, this translates to a corresponding contrast in density between jet and environment.

Overall, these results confirm that centrifugal instability is most effective in fast, weakly magnetized flows and is largely insensitive to the external medium, unless the density contrast is large. This is consistent with the parameter ranges used in our linear analysis, which were motivated by simulation values: Γ = 5–10, σ = 0.001 (almost hydrodynamic) or 0.1 (moderately magnetized), and external enthalpy h_e characteristic of 1–1.5.

3.4. 3D results

3.4.1. The case run A_0.1

As a reference, we adopt a moderate magnetization σ = 0.1, which lies between the nearly hydrodynamic regime (small σ) and strongly magnetized flows prone to kink instabilities (e.g., Barniol Duran et al. 2017). In this regime we focus on CFI instabilities, which should be stable according, for example, to Matsumoto et al. (2021). As already discussed, we use as initial condition for the 3D simulations the steady state reached at the end of the 2D simulations that we present (for case A_0.1) in Figure 2, where we display maps of the total pressure, thermal pressure, magnetic pressure, and total magnetization. These plots focus on the region near the base of the jet and the first three recollimation shocks. We observe that the recollimation shock is strong and produces a distinct radial structure in the jet. At the jet edges, the total pressure is dominated by magnetic pressure, while along the axis, thermal pressure is stronger.

The pressure maps in Figure 2 confirm that during the early free expansion phase, magnetic pressure dominates throughout most of the jet. After recollimation, thermal pressure becomes important, particularly near the axis.

We can perform a linear stability analysis for parameters extracted from the initial jet configuration. In the linear analysis, the instability is studied in the idealized configuration of a rotating cylindrical shell of height 2πR_j, so that perturbations can be decomposed into Fourier modes with axial wavenumber k, Fig. 6, bottom panel. In our simulations, however, we probe the instability through azimuthal structures developing around the jet cross-section, which requires translating the axial description of the linear model into an azimuthal one. We identify the cylinder in the linear model with a circular cross-section of radius R_j in the jet, so that the axial periodicity in the linear model corresponds to the azimuthal periodicity around the jet in the simulation, Fig. 6, upper panel. Under this identification, the axial wavenumber k maps to the azimuthal mode number m through k = m/R_j, so that. Figure 8 shows the linear growth rate for Case A_0.1 (Γ = 10, σ = 0.1, upper profiles of Fig. 5) as a function of wavenumber. Normalizing k with the radius of curvature R_c (as done in the linear analysis) yields kR_c = m(R_c/R_j).

Fig. 8. Growth rate of the CFI as a function of axial wavenumber kRc, computed using the profiles for simulation A_0.1 presented in Fig. 5. Shaded vertical bands indicate the values of kRc corresponding to azimuthal mode numbers m = 1 to m > 6, derived from the simulation jet radius Rj = 0.08 and local curvature radius Rc = 15.28. The most unstable mode occurs at kmaxRc ≈ 190, in agreement with the effective kRc from the simulation for m = 1.

Using the measured values R_c = 15.28 and R_j = 0.08, we obtain R_j/R_c ≈ 0.0053. For the m = 1 mode (which corresponds to the maximum wavelength 2πR_j) this gives kR_c ≈ 191. Modes with a lower values of kR_c have a wavelength that cannot fit around the jet circumpherence and correspond to the green stable region in Fig. 8. The linear analysis further predicts a most unstable mode at k_maxR_c ≈ 1100. Thus, Case A lies well within the unstable range, confirming that the CFI develops physically in the simulation.

The linear stability analysis yields dimensionless growth rates expressed as Im ω R_c/c, where R_c is the curvature radius of the cylindrical shell. In our simulations, R_c ≈ 15z₀, and the natural time unit is z₀/c. For the azimuthal mode number m = 7 − 8, which corresponds to the dominant unstable mode observed in the simulations, the maximum growth rate reaches Im ω R_c/c ≈ 13. Translating this into simulation time units, the growth time is t_grow = 1/Im ω = R_c/(Im ω c)≈15z₀/(13c)≈1.15 z₀/c. This timescale indicates that the centrifugal instability develops in a time slightly longer than a single code unit z₀/c. However, the geometry is more complex than a simple cylinder, as the jet is deformed and the flow is advected along the axis. Taking advection into account, the growth time can be translated into a growth length by multiplying with the flow speed (approximated by c), which is consistent with the width of the recollimation region where the instability develops in the simulations.

Further insights are provided by Fig. 9, which displays, at t = 40 t_cr, longitudinal maps of the Lorentz factor Γ, and total pressure P_tot (left panel), along with a map and transverse cuts of ρ (right panel) at the position indicated by the dashed line in the left panel. In the density plot we can observe the beginning of the development of characteristic finger-like structures, signatures of growing CFI with m = 8, that in fact corresponds to the value in which the linear analysis predicts the maximum growth rate (see Fig. 8).

Fig. 9. First column: Maps of Γ and log Ptot, with transverse cuts of Pth and log ρ at a z-position, for the 3D case A_0.1 with magnetization σ = 0.1. The horizontal line corresponds to these cut, showing the logarithmic density.

Considering that our simulations extend up to t ∼ 100 z₀/c, the observed rapid growth confirms that the instability fully develops and saturates early, in agreement with the nonlinear features and perturbations observed in the flow. Thus, the linear analysis provides insight into the onset and location of the centrifugal instability in our numerical models. The large computational domain (z ≤ 30z₀) and long integration time (t = 100 t_cr) allow us to track the whole evolution of the instability. Figure 10 shows a volume rendering of the z-component of the jet’s four-velocity (Γv_z) from the 3D evolution of case A_0.1, highlighting the emergence of turbulent structures. High-velocity regions and their interactions with surrounding, slower-moving material are clearly visible, indicating the onset of dynamic disruptions. As shown in the pressure evolution of Fig. 11, the jet remains coherent and stable up to approximately t ≈ 40, t_cr (see also Costa et al. 2025 for a detailed discussion of the hydrodynamic case). Beyond this point, the growth of instabilities progressively disrupts the flow. In particular, the CFI continues to grow and is further amplified by the Richtmyer–Meshkov instability (RMI), leading to a coupling between the two modes (Brouillette 2002; Zhou et al. 2021). In plasmas, RMI develops when a shock wave crosses a perturbed interface between fluids of different densities, and in this context it reinforces the action of the CFI.

Fig. 10. 3D volume rendering of the z-component of the jet four-velocity (Γvz), case A_0.1, σ = 0.1 at t = 80. The black-gray bar next to the color bar indicates the opacity: gray regions correspond to opaque values of Γvz, while black regions indicate transparent regions.

Fig. 11. Time evolution of the thermal pressure distribution for case A_0.1 (σ = 0.1) in the x–z plane. The sequence of snapshots (t = 35–40) shows the growth of the centrifugal instability: pressure perturbations at the jet edges gradually amplify, forming the characteristic lateral fingers and oscillatory structures along the jet axis. These features indicate the progressive destabilization of the jet due to centrifugal effects.

Figure 12 shows the plasma β and the toroidal component of the magnetization σ from the 3D simulation. These maps indicate that the jet remains coherent up to the second recollimation point (z ≈ 7z₀), beyond which it gradually decelerates and develops turbulence. The plasma β provides a direct diagnostic of dissipation, distinguishing magnetically dominated regions (low β) from kinetically dominated ones (high β). In the early phase, the jet’s β distribution is relatively uniform, while further downstream it exhibits strong fluctuations driven by nonlinear interactions. A central spine surrounded by a sheath becomes apparent, reflecting the impact of recollimation on the flow structure. The toroidal magnetization shown in the right panel emphasizes the role of the magnetic field in shaping this evolution, revealing how magnetic stresses contribute to both the stability and the eventual disruption of the jet.

Fig. 12. Same as Fig. 9, but including MHD diagnostics (plasma β on the left and the toroidal component of magnetization on the right) for the 3D case A_0.1 with σ = 0.1.

We conclude with one final technical note. In starting a 3D calculation from the configuration produced by a 2D calculation, as we do here, one must be careful. As described in Appendix C, the interpolation of the jet profile from a cylindrical grid to a Cartesian one at the injection boundary may introduce minor numerical artifacts. These arise from small mismatches at the boundary, which act as transient perturbations at the jet base. In marginally unstable cases, the combination of this initial perturbation and the underlying Cartesian symmetry can produce a weak cross-shaped pattern in the jet cross-section that derives from a preferential excitation of grid-aligned perturbations. In more unstable regimes, where physical modes grow faster, such numerical imprints are rapidly suppressed. This perturbation, which we refer to as a “wave,” rapidly propagates up the jet axis, altering the jet structure. In magnetized cases, particularly those with higher σ, this transient wave can, in principle, seed asymmetric distortions of the magnetic field, potentially acting as a trigger for current-driven instabilities such as the kink mode further downstream. However, in the parameter regime explored here –where the focus is on the early development of the CFI– the influence of this perturbation remains limited. After the passage of the wave, the system relaxes back toward its initial equilibrium, and no sustained kink-like growth is observed within the simulated timescales. A more detailed analysis of the long-term evolution and possible onset of kink instabilities will be presented in a forthcoming work.

3.4.2. A stable case

As shown in Sect. 3.1, Case A_0.1 provides insights into unstable configurations. Motivated by linear stability analysis for CFI, which suggests that jets with lower Lorentz factors are more stable, and by the results of Komissarov et al. (2019), which indicate that smaller opening angles suppress instabilities, we select parameters that should favor a stable configuration. In this simulation the adopted parameters are: density ratio ν = 10⁻⁵, Lorentz factor Γ = 5, jet opening angle θ_jet = 0.05, and pressure ratio P_ratio = 10⁻³.

The results confirm that the jet remains stable. To test the robustness of this conclusion, we extended the simulation up to t = 100 t_cr and found no transition to turbulence or structural deformation of the jet. The Lorentz factor remains steady along the axis, with no evidence of deceleration or anomalous pressure structures throughout the domain. An independent check based on the field geometry gives σ_tor/Γ² ≈ 0.0523 and (σ_tor/Γ)/(R_j/R)≈15.5 ≫ 1. This places the system well below the critical threshold for centrifugal instabilities, in agreement with the absence of unstable modes in the simulation. A posteriori, this evidence is confirmed by the linear analysis outlined in Sect. 3.4.1. The maximum growth rate occurs at k_maxR_c ≈ 350. With R_c ≈ 16.86 and R_j ≈ 0.0567, the resulting ratio R_j/R_c ≈ 3.4 × 10⁻³ falls below the numerical threshold (∼0.004) required for the m = 4 mode to develop, confirming that the configuration is linearly stable, as the linear analysis shows in Fig. 13. As shown in Figs. 14 and 15, the jet evolves smoothly without signs of disruption or turbulence.

Fig. 13. Top panel: Comparison of radial profiles of toroidal magnetization σtor (blue), Lorentz factor Γ (orange), rest-mass density ρ (green), and enthalpy h/ρ (red) at z = zmaxComparison of σtor/Γ2 for jets with σ = 0.1 with an opening angle θ = 0.05, Lorentz factor Γ = 5. Bottom panel: Growth rate of the CFI as a function of axial wavenumber kRc. Shaded vertical bands indicate the values of kRc corresponding to azimuthal mode numbers m = 1 and m = 2 is, derived from the simulation jet radius Rj = 0.0567 and local curvature radius Rc = 16.86. The most unstable mode occurs at kmaxRc ≈ 300.

Fig. 14. Comparison of σtor/Γ2 for jets with σ = 0.1 with a. opening angle θ = 0.05, Lorentz factor Γ = 5 (left panel) and b. θ = 0.1, Lorentz factor Γ = 10 (left panel) at t = 0.

These parameters are in the stable regime and we further verified stability by increasing the jet density (i.e., varying ν to higher values; see Appendix D). The jet remains collimated and stable across these different setups. This indicates that suppression of CFI in this parameter regime is robust and not sensitive to moderate changes in ν.

The stabilizing effect is evident in Figs. 14 and 15 (left panels), which show the distribution of σ_tor/Γ², streamlines, and plasma β. These diagnostics demonstrate that magnetic tension dominates over destabilizing forces, preventing the growth of unstable modes. A direct comparison with the unstable case underscores the crucial role of geometry and magnetization in controlling jet stability. Together, these results show that narrower opening angles and stronger magnetic support shift the system away from the CFI threshold, enabling the jet to maintain coherence over large distances.

Fig. 15. Comparison β for jets with σ = 0.1 with a. opening angle θ = 0.05, Lorentz factor Γ = 5 (left panel) and b. θ = 0.1, Lorentz factor Γ = 10 (left panel) at t = 80.

4. Discussion and conclusions

In this work we examine how the addition of a magnetic field influences the stability of relativistic jets launched into an external medium with a vertical pressure profile that causes the jet to recollimate. In 2D simulations, where axial symmetry is enforced, a stable sequence of recollimation shocks develops and the jet flow remains ordered. In a more realistic 3D calculation, non-axisymmetric modes and instabilities are now allowed to develop and may disrupt the jet. Indeed, earlier hydrodynamic calculations, for example, Costa et al. (2025), show that the curvature of the jet fluid flow lines that near a recollimation shock often enables the development and significant growth of the CFI. This in turn can lead to very efficient disruption of the jet flow after just one or two recollimation shocks, which may be relevant to the structure of jets observed in FRI and FR0 radio sources.

The addition of an internal jet magnetic field to the problem of a recollimating jet can have three important effects on the subsequent evolution of the jet. First, as we demonstrate in Sect. 2, the addition of a magnetic field changes the vertical pressure profile of the jet, which in turn changes the position of the recollimation shocks and the curvature of the flow lines near the shocks. This quantitatively changes the values of the jet injection parameters that lead to a CFI-unstable flow. Second, to the extent that fluid flow lines and magnetic field lines are tied together in the MHD approximation, the development of the CFI “corrugates” and stretches the magnetic field lines, increasing the effective tension of the field lines and thus resisting the further development of the CFI. A sufficiently strong field may thus quench the CFI and allow the jet propagate in a stable manner. Third, for sufficiently strong fields, other instabilities such as the current-driven instability (CDI) eventually become important and may again disrupt the jet. In this study we restrict ourselves to the case of low to moderate jet magnetizations (σ ∼ 0 − 1), where the CFI can still operate. For such magnetizations, our 3D results indicate that the CFI still seems to be the fastest growing initial instability, and thus the one relevant for determining where the jet flow eventually becomes unstable. Furthermore, we present a case that remains stable without triggering the CFI, indicating that the underlying stability problem is more complex than it may initially appear.

This conclusion is consistent with Matsumoto et al. (2021), who also emphasized the relevance of the CFI in magnetized flows. However, we find that the problem of determining exactly which initial jet parameters will lead to unstable jet flow is not as straightforward as implied in Komissarov et al. (2019). For instance, a 3D case in Matsumoto et al. (2021) considered stable by the Komissarov et al. (2019) criterion becomes unstable when extended to a larger domain and longer runtime. A further complication is that the nonlinear evolution of the jet after the CFI begins appears nontrivial. In our 3D calculations of the jet flow evolution, the CFI appears to trigger other instabilities, and as in the hydrodynamic calculation of Costa et al. (2025), different regions of the jet can eventually influence each other, possibly leading to mixing of jet material with the external medium and significantly changing the pressure profile of the medium near jet (as in the right panel of Fig. 15). Again, one will not see this long-term behavior if the simulation is truncated in space or time. Conversely, if one sees an intrinsically young jet system in nature that does not show evidence (yet) of disruption, one should not conclude that the jet parameters are such that the jet is stable in the long-term.

Our 2D survey shows that the location of the first recollimation shock is sensitive to the choice of initial jet parameters, namely the density contrast ν, the pressure ratio P_ratio, the jet magnetization and the relative strength of the toroidal and poloidal field components, the Lorentz factor Γ, and the jet opening angle θ_jet. Increasing the value of the nonmagnetic field parameters tends to push the recollimation point to larger distances, modifying the jet’s geometry. Physically, this reflects the enhanced inertial and pressure content of the flow, which requires more expansion and interaction with the external medium before reconfinement occurs. By contrast, higher magnetization moves the recollimation shock closer to the source. These trends may offer a valuable guide for interpreting observed structures in AGNs jets, where knot positions are often correlated with recollimation features.

In summary, magnetization is a key factor in controlling jet collimation and the development of instabilities in the jet. While 2D studies and linear theory offer important guidance, only 3D simulations capture the full nonlinear evolution. These are currently computationally expensive to run for long enough times and with sufficient spatial resolution in a sufficiently large domain. We therefore presented only a few illustrative 3D calculations, and further exploration of the full 3D problem is warranted.

Our results have direct implications for understanding the structure of AGNs jets, the origin of knots and variability, and their polarized emission signatures. For example, the development of the CFI-driven turbulence and magnetic field disordering could strongly influence the observed polarized emission (Sciaccaluga et al. 2025). In particular, the transition from a magnetically dominated spine to a turbulent sheath, or the radial mixing induced by recollimation and instability growth, can produce variable polarization degrees and swings. This is especially relevant in light of recent polarimetric observations from instruments like IXPE. Future work should focus on connecting the internal instability-driven structure of the jet with synthetic polarization maps to explore these effects in more detail.

Acknowledgments

We thank the anonymous referee for their useful comments. We thank Brian Reville, Konstantinos Nektarios Gourgouliatos, Dimitrios Millas, Paola Rossi and Om Sharan Salafia for fruitful discussions – special thanks to Om for Fig. 6 – and SB acknowledges helpful manuscript feedback from Brian and Dimitrios. We gratefully acknowledge financial support by INAF Theory Grant 2022 (PI F. Tavecchio) and NASA IXPE Theory Grant80NSSC24K1177 (P.I. P. Coppi). This work was also funded by the European Union-Next Generation EU, PRIN 2022 RFF M4C21.1 (2022C9TNNX). We acknowledge support by CINECA, through ISCRA and Accordo Quadro INAF-CINECA, and by PLEIADI, INAF – USC VIII, for the availability of HPC resources (PI S. Boula). G.B. and A.C. acknowledge the support by the Spoke-1 “FutureHPC & BigData of the ICSC – Centro Nazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing – funded by European Union – NextGenerationEU. We used the following Python libraries: Numpy (Harris et al. 2020), Matplotlib (Hunter 2007), Scipy (Virtanen et al. 2020), and PyPluto Mattia et al. (2025).

References

Appendix A: Smoothed profiles and initial conditions

The simulations adopt smoothed radial profiles to ensure a continuous transition between the jet and ambient medium, thereby minimizing numerical artifacts at the jet boundary. These profiles are implemented using the hyperbolic secant (sech) function, which allows for smooth yet steep gradients.

The general expression for a smoothed quantity q is given by:

$$\begin{array}{c}\hfill q(r,z)=q_{\text{ext}}+(q_{\text{j}}-q_{\text{ext}})\text{sech}\left[{\left(\frac{r}{z{\theta }_q}\right)}^{{\alpha }_q}\right],\end{array}$$(A.1)

where q_ext and q_j are the values of the quantity in the ambient medium and jet, respectively; θ_q sets the angular width of the transition layer, and α_q controls the steepness of the profile.

The cylindrical radial distance r and spherical radius r′ are defined as:

$$\begin{array}{cc}\hfill r& =\sqrt{x^2+y^2},\hfill \end{array}$$(A.2)

$$\begin{array}{cc}\hfill r^{\prime }& =\sqrt{x^2+y^2+z^2}.\hfill \end{array}$$(A.3)

The Lorentz factor profile is initialized as:

$$\begin{array}{c}\hfill \mathrm{\Gamma }(r,z)=1+({\mathrm{\Gamma }}_j-1)\text{sech}\left[{\left(\frac{r}{z{\theta }_{\mathrm{\Gamma }}}\right)}^{{\alpha }_{\mathrm{\Gamma }}}\right],\end{array}$$(A.4)

from which the three-velocity magnitude is calculated using:

$$\begin{array}{c}\hfill v_0=\sqrt{1-\frac{1}{{\mathrm{\Gamma }}^2}}.\end{array}$$(A.5)

The density profiles are given by:

$$\begin{array}{cc}\hfill {\rho }_{\text{ext}}(z)& =z^{-\eta },\hfill \end{array}$$(A.6)

$$\begin{array}{cc}\hfill {\rho }_j(r^{\prime })& =\nu r^{\prime -2},\hfill \end{array}$$(A.7)

$$\begin{array}{cc}\hfill {\rho }_f(r,z)& ={\rho }_{\text{ext}}+({\rho }_j-{\rho }_{\text{ext}})\text{sech}\left[{\left(\frac{r}{z{\theta }_{\rho }}\right)}^{{\alpha }_{\rho }}\right],\hfill \end{array}$$(A.8)

where ν is the density contrast between the jet and the ambient medium, and η controls the stratification of the external atmosphere.

Similarly, the pressure profiles follow:

$$\begin{array}{cc}\hfill p_{\text{ext}}(z)& =P_0{z}^{-\eta },\hfill \end{array}$$(A.9)

$$\begin{array}{cc}\hfill p_j(r^{\prime })& =P_0{r}^{\prime -2\gamma }{P}_{\text{ratio}},\hfill \end{array}$$(A.10)

$$\begin{array}{cc}\hfill p_f(r,z)& =p_{\text{ext}}+(p_j-p_{\text{ext}})\text{sech}\left[{\left(\frac{r}{z{\theta }_p}\right)}^{{\alpha }_p}\right],\hfill \end{array}$$(A.11)

with P_ratio denoting the jet-to-ambient pressure ratio.

For points inside the jet cone, defined by r/z < θ_j, the velocity components are initialized as:

$$\begin{array}{c}\hfill u_x=v_0\frac{x}{r^{\prime }},u_y=v_0\frac{y}{r^{\prime }},u_z=v_0\frac{z}{r^{\prime }}.\end{array}$$(A.12)

The parameters used to define the smoothed profiles and physical conditions in the simulations are: θ_j = 0.1, Γ_j = 10, P₀ = 3 × 10⁻⁶, θ_Γ = 0.07, α_Γ = 13, θ_p = 0.13, α_p = 15, θ_ρ = 0.1, α_ρ = 15, θ_β = 0.004, b_z = 10⁻⁵, and α_z = 10⁻¹⁰.

Appendix B: Resolution of the simulations

The 2D simulations use a grid of 1400 × 2200 zones. In the central region [0, 2]×[1, 20], the grid is uniform, providing a resolution of 50 points per initial jet radius, while the outer regions are geometrically stretched. The X1-grid comprises two regions: an inner uniform region r ∈ [0, 1000] with 400 zones, and an outer stretched region with up to 6 zones. Similarly, the X2-grid has two regions: an inner uniform region z ∈ [1, 1500] with 700 zones, and an outer stretched region up to 30. At the final time, the 2D simulations reach a nearly steady state, which is used as the initial condition for the 3D runs.

The 3D simulations are performed on a Cartesian domain with coordinates

$$x\in [-\frac{L_{x,3D}}{2},\frac{L_{x,3D}}{2}],y\in [-\frac{L_{y,3D}}{2},\frac{L_{y,3D}}{2}],z\in [1,L_{z,3D}],$$

where lengths are expressed in units of z₀, L_x = L_y, and z is the jet propagation direction. A uniform grid spacing is adopted along z in the central jet region, while stretched grids are used in the transverse directions and in the far axial domain. Table B.1 summarizes the resolution, grid structure, and number of cells per initial jet radius for the cases discussed in Sec. 3.4 (A1 and B1) and in Appendix D (A2 and B2).

Appendix C: Comparison 2D and 3D

To initialize the 3D simulations from 2D axisymmetric profiles, we employed the InputDataInterpolate and StaggeredRemap utility functions provided by PLUTO. These tools ensure consistent remapping of primitive and staggered quantities from 2D cylindrical coordinates to the 3D Cartesian grid. Figure C.1 shows a comparison of the pressure profiles–thermal, magnetic, and total–evaluated at the same axial position z = z₀. The agreement between the 2D and 3D profiles confirms the accuracy of the interpolation procedure. Minor discrepancies appear in the thermal pressure, particularly for x < 0.1, where the 3D profile is slightly smoother. This is due to the interpolation onto a lower-resolution 3D grid, which tends to average out sharp features. Nonetheless, the essential physical structure is preserved, validating the initialization approach.

Figure C.2 shows the Lorentz factor at t = 0 and t = 40, illustrating that the first recollimation point remains unaffected by the initial numerical wave, which acts as a perturbation. The observed differences at the second recollimation point are instead due to the physical evolution of the system, as this is the time when the instability begins to develop.

Fig. C.1. Comparison of thermal, magnetic, and total pressure profiles along the x-axis at z = z0, for the 2D simulation at t = 3000 and the 3D simulation at t = 0. The 3D data have been interpolated onto the 2D spatial grid for a consistent comparison.

Fig. C.2. Comparison of the Lorentz factor at t = 0 and t = 40 for case A_0.1 with σ = 0.1, illustrating the jet structure before and after the impact of the imposed wave.

Appendix D: 3D simulations: other parameters

The appendix provides additional diagnostics of the jet evolution. We show logarithmic normalized profiles along the z-axis at x = 0 for the density ρ, Lorentz factor Γ, and total pressure P_tot, together with maps and transverse cuts of Γ and log P_tot at several positions along the jet. These allow us to compare the longitudinal evolution of the main dynamical quantities with the corresponding cross-sectional structure, highlighting differences between unstable and stable configurations. Figures D.1 and D.2 show the evolution of two jets with different initial parameters. In Fig. D.1, the jet is initialized with magnetization σ = 10⁻³, density ratio ν = 10⁻⁵, Lorentz factor Γ = 10, opening angle θ_jet = 0.1, and pressure ratio P_ratio = 10⁻³. This configuration leads to a rapid destabilization, closely resembling the behavior seen in hydrodynamic jets. The structure breaks down early in the simulation, with no clear recollimation features, indicating that the magnetic support is insufficient to suppress growing instabilities.

By contrast, the jet in Fig. D.2, though having σ = 10⁻¹, is initialized with a higher density ratio ν = 10⁻⁴, a lower Lorentz factor Γ = 5, and a smaller opening angle θ_jet = 0.05. This jet remains stable over a longer distance and preserves its structure. Longitudinal profiles show smoother trends in ρ, Γ, and P_tot, while transverse cuts of Γ and log P_tot reveal only mild deformations. These distortions are not associated with CFI, but rather suggest the presence of weaker or subdominant modes. In the limit where the jet remains stable, the cuts reveal a weak but coherent m = 4 signature, likely related to the underlying grid symmetry. However, this grid-induced mode is naturally enhanced due to symmetry and grows faster than physical instabilities, such as the CFI. Although present, it does not affect the overall stability or evolution of the jet. The increased stability is attributed not only to the slightly higher effective magnetization but also to the narrower opening angle, which reduces the jet’s initial expansion and limits the triggering of instabilities. This highlights how geometric and dynamical parameters jointly influence the onset and growth of jet instabilities.

The effect of the imposed wave perturbation also differs significantly between these two cases. In the stable jet configuration (Fig. D.2), the wave acts as a mild perturbation rather than a destabilizing driver. Here, the jet’s conditions effectively suppress the CFI, and the wave does not significantly amplify instability growth. However, grid-related numerical instabilities can still manifest, which explains the appearance of four discrete modes observed in the transverse cuts. These modes correspond to perturbations seeded by the numerical grid rather than physical instabilities, and crucially, they do not disrupt the jet’s overall stability. Additionally, the presence of jet rotation induces a kink-like deformation, visible in the cuts; however, this kink instability remains benign and does not interrupt the jet flow. Conversely, in the unstable jet case (Fig. D.1), the jet is inherently prone to CFI due to its initial conditions. Here, the imposed wave perturbation has little to no influence on the jet’s evolution; the CFI dominates and drives the rapid destabilization observed. The intrinsic CFI dominates the jet dynamics, rendering the wave perturbation negligible in driving the disruption of the jet structure.

These results highlight the intricate interplay between physical instabilities and perturbations. While waves can seed or amplify instabilities in marginally stable jets, they are insufficient to destabilize strongly unstable configurations dominated by current-driven modes.

Fig. D.1. First column: logarithmic normalized profiles along the z-axis at position x = 0 for density ρ, Lorentz factor Γ, and total pressure Ptot for the 3D case A with magnetization σ = 0.001. The other columns include the map of Γ and log Ptot, and cuts at three different z-positions, where the lines in the profiles and maps correspond to these cuts, showing the Lorentz factor and logarithmic density.

Fig. D.2. Same as in Fig. D.1. For the 3D case with magnetization σ = 0.1, θjet = 0.05, ν = 10−4 and Γ = 5.

All Tables

All Figures

	Fig. 1. Profiles along the z-axis at position x = 0 of the normalized Lorentz factor for the 2D simulations presented in Table 1. The normalization is taken at the position (0, 0) for simulations A, B, C, and D with magnetization σ = 0, 0.1, and 1.
In the text

	Fig. 2. Logarithmic maps of total pressure, thermal pressure, magnetic pressure, and total magnetization for the 2D simulation A_0.1 with magnetization σ = 0.1. The initial parameters are: density ratio ν = 10−5, Lorentz factor Γ = 10, jet opening angle θjet = 0.1, and pressure ratio Pratio = 10−3. The maps correspond to the final time step t = 3000 tcr in units of z0/c.
In the text

Fig. 3. Streamlines (initialized from the same starting point) overlaid on the spatial distribution of σtor/Γ2 for all simulation runs. This ratio serves as a local diagnostic for the onset of centrifugal instability, with low values marking regions where the jet is most susceptible to curvature-driven perturbations. In the widest jets (e.g., Case D_0.1), outer streamlines would likely reach regions satisfying the CFI criterion; here we show a representative streamline to illustrate the overall trend.

In the text

	Fig. 4. Left: Streamlines for cases A_0.1–D_0.1 (as in Fig. 3), starting at x0 = 0.082, 0.098, 0.091, and 0.124 respectively, with curvature and the ratio Rj, max/Rc. Right: Profiles of σtor/Γ2 versus x at z = zxmax, with vertical lines marking xmax = Rj, max.
In the text

Fig. 5. Comparison of radial profiles of toroidal magnetization σtor (blue), Lorentz factor Γ (orange), rest-mass density ρ (green), and enthalpy h/ρ (red) at z = zmax for Case A_0.1 (top) and Case D_0.1 (bottom). Case A_0.1 shows smooth, monotonic profiles, whereas Case D_0.1 exhibits sharp gradients and non-monotonic behavior, underscoring the complexity of its internal structure. The vertical dashed line marks the position of xmax determined from the streamline analysis. Profiles are plotted on a logarithmic scale.

In the text

Fig. 6. Cartoon illustrating the configurations used in the linear analysis, adapted from Komissarov et al. (2019). Top panel: Relativistic jet with a purely azimuthal magnetic field undergoes recollimation due to the confining action of an external pressure. The white inner region represents the freely expanding, unshocked flow, while the shaded area corresponds to the shocked outer layer. The interface between the two is the conical reconfinement shock driven into the jet. Bottom panel: < agnetised, rotating cylindrical shell with a purely axial magnetic field, representing the idealised configuration analysed in this work.

In the text

	Fig. 7. Results from the linear analysis of the CFI. Top panel: Growth rate as a function of kRc, for varying Lorentz factors and magnetizations, assuming fixed external conditions (he ≈ 1). Bottom panel: Same as top, but for fixed magnetization σ = 0.1, and varying external enthalpy ratio he and Lorentz factor Γ.
In the text

Fig. 8. Growth rate of the CFI as a function of axial wavenumber kRc, computed using the profiles for simulation A_0.1 presented in Fig. 5. Shaded vertical bands indicate the values of kRc corresponding to azimuthal mode numbers m = 1 to m > 6, derived from the simulation jet radius Rj = 0.08 and local curvature radius Rc = 15.28. The most unstable mode occurs at kmaxRc ≈ 190, in agreement with the effective kRc from the simulation for m = 1.

In the text

	Fig. 9. First column: Maps of Γ and log Ptot, with transverse cuts of Pth and log ρ at a z-position, for the 3D case A_0.1 with magnetization σ = 0.1. The horizontal line corresponds to these cut, showing the logarithmic density.
In the text

	Fig. 10. 3D volume rendering of the z-component of the jet four-velocity (Γvz), case A_0.1, σ = 0.1 at t = 80. The black-gray bar next to the color bar indicates the opacity: gray regions correspond to opaque values of Γvz, while black regions indicate transparent regions.
In the text

Fig. 11. Time evolution of the thermal pressure distribution for case A_0.1 (σ = 0.1) in the x–z plane. The sequence of snapshots (t = 35–40) shows the growth of the centrifugal instability: pressure perturbations at the jet edges gradually amplify, forming the characteristic lateral fingers and oscillatory structures along the jet axis. These features indicate the progressive destabilization of the jet due to centrifugal effects.

In the text

	Fig. 12. Same as Fig. 9, but including MHD diagnostics (plasma β on the left and the toroidal component of magnetization on the right) for the 3D case A_0.1 with σ = 0.1.
In the text

Fig. 13. Top panel: Comparison of radial profiles of toroidal magnetization σtor (blue), Lorentz factor Γ (orange), rest-mass density ρ (green), and enthalpy h/ρ (red) at z = zmaxComparison of σtor/Γ2 for jets with σ = 0.1 with an opening angle θ = 0.05, Lorentz factor Γ = 5. Bottom panel: Growth rate of the CFI as a function of axial wavenumber kRc. Shaded vertical bands indicate the values of kRc corresponding to azimuthal mode numbers m = 1 and m = 2 is, derived from the simulation jet radius Rj = 0.0567 and local curvature radius Rc = 16.86. The most unstable mode occurs at kmaxRc ≈ 300.

In the text

	Fig. 14. Comparison of σtor/Γ2 for jets with σ = 0.1 with a. opening angle θ = 0.05, Lorentz factor Γ = 5 (left panel) and b. θ = 0.1, Lorentz factor Γ = 10 (left panel) at t = 0.
In the text

	Fig. 15. Comparison β for jets with σ = 0.1 with a. opening angle θ = 0.05, Lorentz factor Γ = 5 (left panel) and b. θ = 0.1, Lorentz factor Γ = 10 (left panel) at t = 80.
In the text

	Fig. C.1. Comparison of thermal, magnetic, and total pressure profiles along the x-axis at z = z0, for the 2D simulation at t = 3000 and the 3D simulation at t = 0. The 3D data have been interpolated onto the 2D spatial grid for a consistent comparison.
In the text

	Fig. C.2. Comparison of the Lorentz factor at t = 0 and t = 40 for case A_0.1 with σ = 0.1, illustrating the jet structure before and after the impact of the imposed wave.
In the text

Fig. D.1. First column: logarithmic normalized profiles along the z-axis at position x = 0 for density ρ, Lorentz factor Γ, and total pressure Ptot for the 3D case A with magnetization σ = 0.001. The other columns include the map of Γ and log Ptot, and cuts at three different z-positions, where the lines in the profiles and maps correspond to these cuts, showing the Lorentz factor and logarithmic density.

In the text

	Fig. D.2. Same as in Fig. D.1. For the 3D case with magnetization σ = 0.1, θjet = 0.05, ν = 10−4 and Γ = 5.
In the text