© The Authors 2024

1. Introduction

Active galactic nuclei (AGN) power relativistic outflows that propagate up to hundreds and even thousands of kiloparsecs (Begelman et al. 1984; Blandford et al. 2019). A fundamental ingredient for their capacity to extend through nine orders of magnitude in distance is their inertia (e.g., Perucho 2019, and references therein), which is achieved by means of the investment of energy into the acceleration of the flow to relativistic velocities. The mechanisms behind the establishment of such outflows have been widely studied over recent years, and substantial progress has been made (e.g., Vlahakis & Königl 2004; Komissarov et al. 2007). However, such studies have focused on the assumption that acceleration is magnetically driven, ignoring thermal acceleration.

The crucial region in which this process can be explored in the jet is the so-called acceleration and collimation region. There, jets collimate, evolving from a parabolic geometry (i.e., r ∝ z^∼0.5, where z is the radial distance from the core and r is the jet radius) to a conical one (r ∝ z^∼1), and accelerate up to high relativistic velocities along the same distance (see e.g., Boccardi et al. 2017, and references therein). Therefore, in cold jets, the nature of the collimation and acceleration phenomena seems to be strictly correlated. The current common view of this process suggests that equipartition between magnetic field and particle energy is reached (Nokhrina et al. 2019, 2020) by the end of the acceleration region, and that external conditions, such as the crossing of the Bondi radius (Kovalev et al. 2020), may play a relevant role in shaping the jet expansion profile. It is often considered that equipartition refers to an equilibrium between magnetic and kinetic energy densities in the jet, but only the equipartition between magnetic and internal energies can be inferred, that is, from synchrotron radiative output. It is therefore plausible that even though magnetic and internal energies are in or close to equipartition at the end of the acceleration region, it is the kinetic energy that dominates jet dynamics from that point on.

Jets in Fanaroff-Riley I (FR I; Fanaroff & Riley 1974) radio galaxies, which are used as an example in the present study, are expected to reach maximum Lorentz factors of γ_max ∼ 10, as seen from statistical studies of BL Lac objects, their blazar counterparts (Hovatta et al. 2009). As mentioned above, acceleration can be related to two main possible sources of energy (see, e.g., Komissarov 2012): (i) thermal energy via the Bernoulli mechanism, or (ii) magnetic energy via the differential collimation mechanism. The former allows the jet to reach a maximum Lorentz factor γ_max = γ[1 + 4ε/(3c²)], where ε is the jet-specific internal energy at injection and γ is the initial Lorentz factor, whereas for the latter, γ_max = γ[1 + B²/(4πρc²)],  where B and ρ are the jet magnetic field and rest-mass density, respectively.

The current paradigm of magnetic jet acceleration relies on the theoretical work of Vlahakis & Königl (2003a,b, 2004). The authors suggest that thermal acceleration would be relevant only at very compact scales, and magnetic acceleration would take over and prolong the accelerating process up to parsec scales. Magnetic acceleration was further studied by Komissarov et al. (2007), who showed how the conversion of magnetic energy is a viable way to accelerate the bulk flow by means of numerical simulations. A relevant aspect revealed by these works is that magnetic acceleration does not occur when the outflow is following a conical expansion, but becomes effective when the jet expands following a parabolic law. In this process, the toroidal field acts not only as a collimating agent, but also as the main driver of jet bulk acceleration. The basic theory behind this process is summarized in Komissarov (2012).

However, this view is founded on a more fundamental assumption that the jet internal energy is relatively small in comparison to the total jet energy flux. Vlahakis & Königl (2004) used classical expressions to derive this conclusion, which implicitly means that the rest-mass energy of the particles is significantly larger than their internal energy. It is nevertheless unclear whether or not jet thermodynamics can be elucidated without considering the relativistic nature of jets (see, e.g., Perucho et al. 2017). On the one hand, jets can be thermodynamically relativistic (i.e., ε ≥ c²). On the other hand, strong jet expansion can be prevented by either the ambient pressure or the magnetic tension provided by a toroidal field, which can keep the jet plasma at relativistic internal energies, avoiding its cooling to ε ≪ c². This is, at least, observed in numerical simulations (e.g., Martí et al. 2016; Moya-Torregrosa et al. 2021; Anglés-Castillo et al. 2021). As a consequence, internal energy could still be relevant for the overall process of jet acceleration up to parsec scales. Thermal acceleration (together with magnetic acceleration) in jets is suggested in López-Miralles et al. (2022) for the case of outflows from high-mass X-ray binaries. In this paper, we relax the assumption of cold jets and explore the different roles that both accelerating mechanisms can play.

As an initial template for the physical conditions, we use the prototypical radio galaxy NGC 315. NGC 315 is a nearby (z = 0.0165, Trager et al. 2000) giant FR I (e.g., Laing et al. 2006) radio galaxy whose acceleration and collimation properties have been extensively studied by employing very long baseline interferometry (VLBI) observations. Such observations show that the jet in NGC 315 accelerates up to v ∼ (0.8 − 0.9) c, and collimates on similar scales, with the highest Lorentz factor γ ∼ 4 reached at ∼2 pc (Boccardi et al. 2021; Park et al. 2021; Ricci et al. 2022). Other physical properties, such as the magnetic field strength, the jet power, and the jet width, have been estimated (Boccardi et al. 2021; Ricci et al. 2022) and are also used as initial conditions in our models. A complete characterization of the jet at injection would require initial values for density and pressure. However, they are unknown and we use them as free parameters, allowing us to vary the jet configuration and dominant energy channel at the boundary condition of our numerical simulations at the initial state.

The goal of the present paper is twofold: while exploring the acceleration in different types of jets, namely hot and cold, and the interplay between the separate energy channels along the propagation, we aim to reproduce the final part of the acceleration in NGC 315 and to test whether or not the inferred physical properties on subparsec scales allow the jet to remain stable on parsec scales and to accelerate up to γ ∼ 4. To this end, we investigate a large number of setups using the 2D relativistic magnetohydrodinamical (RHMD) code presented in Martí (2015b) and Martí et al. (2016).

The paper is structured as follows. In Sect. 2 we describe the numerical setup; in Sects. 3 and 4 we present our results and discuss them, respectively; and in Sect. 5 we draw our conclusions.

2. Numerical simulations

2.1. Numerical methods

The simulations presented in this paper are performed using the conservative, 2D axisymmetric RMHD code presented in Martí (2015b). While for an in-depth study on the jet propagation, a 3D approach could allow us to study the phenomenon with fewer assumptions (allowing for example to explore the effects of the development of helical modes, see, e.g., López-Miralles et al. 2022), we consider the 2D axisymmetric approach well suited for our goals. Indeed, jet acceleration is connected with jet expansion, which is an axisymmetric process. Furthermore, a 2D grid allows to test a large number of models due to the lower computational resources required.

Regarding the effects of our assumption on jet stability, the instability modes that we observe are mostly associated with small-scale structures (see Sect. 3). In the long term, the oscillation of the jet outer radius could trigger the development of pinching modes, which are, in general, not disruptive for relativistic flows (e.g., Perucho et al. 2005). Moreover, the dynamics of the simulated jets at injection is completely dominated by expansion and, in such a situation, the growth of instability modes is negligible (Hardee 1986).

In our code, we employ a second-order Godunov-type scheme with the HLL Riemann solver (Harten et al. 1983) and the piecewise linear method (PLM) for cell reconstruction, where the VanLeer slope limiter preserves monotonicity (van Leer 1974; Mignone & Bodo 2006). The limiter is degraded to MinMod (Roe 1986) when strong shocks are detected in order to avoid numerical oscillations. Time integration follows the third-order TVD-preserving Runge-Kutta scheme of Shu & Osher (1989) with CFL = 0.2. The relativistic correction algorithm CA2 of Martí (2015a) is used to correct the conserved variables after each time iteration. The magnetic field divergence-free constraint is preserved with the constrained transport (CT) method (Balsara & Spicer 1999).

2.2. Model assumptions and transversal equilibrium

We study the acceleration problem by means of the ideal RMHD equations, to which we impose axisymmetry. We describe the system using cylindrical coordinates (r, ϕ, z) where r, ϕ and z are the radial, azimuthal and axial coordinates. Axisymmetry allows us to drop the dependence on the azimuthal cylindrical coordinate in the RMHD equations. We assume both the jet and ambient plasma to be perfect gas with a constant adiabatic index of Γ = 4/3. Taking into account that the ambient medium is used as a passive element in the simulations and that the jet plasma is hot, we consider this approach sufficient for our purposes.

The code evolves the rest-mass, momentum and total energy densities, and the magnetic field in a conservative way. Besides this, and in order to minimize the perturbation of the flow, the jet is injected in conditions of transversal equilibrium. As a consequence, the radial components of the magnetic field and the flow velocity, B^r(r), v^r(r), at injection are set to zero, and the equilibrium solutions are characterized by six functions, namely the density and pressure, ρ(r), p(r), and the two remaining components of the flow velocity, v^ϕ(r), v^z(r), and of the magnetic field, B^ϕ(r), B^z(r).

In what follows, all the quantities will be expressed in code units, in which the initial jet radius, r_j, the ambient density at the jet base, ρ_a, and the speed of light, c, are equal to 1 (and a factor $1/\sqrt{4\pi }$ is absorbed in the definition of the magnetic field).

The transversal equilibrium condition is given by a single ordinary differential equation (see e.g., Martí 2015b):

$$\begin{array}{c}\hfill \frac{\mathrm{d}{p}^{\ast }}{\mathrm{d}r}=\frac{\rho h^{\ast }{\gamma }^2{(v^{\varphi })}^2-{(b^{\varphi })}^2}{r},\end{array}$$(1)

in which, $\gamma =1/\sqrt{1-v^i{v}_i}$ (i = r, ϕ, z, and summation over repeated indices is assumed) is the Lorentz factor, p^* and h^* are the total (gas and magnetic) pressure and total specific enthalpy of the jet, namely p^* = p + b²/2 and h^* = 1 + ε + p/ρ + b²/ρ. In the previous expressions, b² = b^μb_μ (μ = t, r, ϕ, z), where b^μ is the magnetic field four-vector in the fluid rest frame. Its components are defined as b⁰ = γBⁱv_i and bⁱ = Bⁱ/γ + b⁰vⁱ which leads to

$$\begin{array}{c}\hfill b^2=\frac{B^2}{{\gamma }^2}+{(B^i{v}_i)}^2.\end{array}$$(2)

Moreover, we assume constant initial density and axial velocity across the jet, and a nonrotating jet model

$$\begin{array}{cc}& \rho (r)=\{\begin{array}{c}{\rho }_j,0\le r\le 1\hfill \\ 1,r>1,\hfill \end{array}\hfill \end{array}$$(3)

$$\begin{array}{cc}& v^z(r)=\{\begin{array}{c}{v}_j^z,0\le r\le 1\hfill \\ 0,r>1,\hfill \end{array}\hfill \end{array}$$(4)

and v^ϕ(r) = 0. With these definitions and the corresponding profiles for the azimuthal and axial components of the magnetic field, B^ϕ(r), B^z(r), Eq. (1) gives the equilibrium pressure profile of the jet at injection.

2.3. Physical model

In code units, the total jet power is given by

$$\begin{array}{c}\hfill L_j=\pi [{\rho }_j{h}_j{\gamma }_j^2+{(B_j^{\varphi })}^2]v_j^z.\end{array}$$(5)

The specific enthalpy is defined as:

$$\begin{array}{c}\hfill h_j=1+{\epsilon }_j+p_j/{\rho }_j\end{array}$$(6)

with pressure, density and specific internal energy related through an ideal gas equation of state, ε_j = p_j/[ρ_j(Γ − 1)], where Γ is the corresponding adiabatic index.

The different energy flux channels that account for the total jet power (Eq. (5)) are computed as follows:

• Magnetic: $F_{\text{M}}={(B_j^{\varphi })}^2{v}_j^z$;

• Internal: $F_{\text{I}}={\rho }_j{\gamma }_j^2{v}_j^z(h_j-1)$;

• Kinetic: $F_{\text{K}}={\rho }_j{\gamma }_j^2{v}_j^z$.

In addition, we compute the total rest mass energy as $F_{\text{R}}={\rho }_j{\gamma }_j{v}_j^z$. As far as the magnetic field configuration is concerned, we use two different approaches: the nonforce-free configuration, as defined in Martí (2015b), and the force-free configuration as defined in Moya-Torregrosa et al. (2021). The two models are described in Sect. 2.3.1 and in Sect. 2.3.2, respectively.

We note that the transition of the initial magnetic field and pressure profiles at the jet/ambient medium interface is smoothed out by convolving the original profiles with a continuous function mimicking a shear layer. In this way, we avoid discontinuities and abrupt changes in our initial conditions, making them more stable. As shear layer, we consider the function f(r) = 1/cosh(r^m) (Bodo et al. 1994; Perucho et al. 2004), in which m = 2, 4, and 8 with m = 8 being the value used the most in our models.

Since our simulations are axisymmetric, we assume reflecting boundary conditions (b.c.) on the axis (r = 0), while we assume zero gradient b.c. at large values of the radial coordinate. Along the z direction, we use inflow/outflow b.c. at the inner boundary, for the jet and the ambient medium, respectively, and outflow b.c. through the outer one.

In this study, we do not take into account the effects associated with resistivity and, consequently, magnetic reconnection. Although the length scales at which reconnection occurs are expected to be significantly smaller than those addressed in this work, its impact might still be relevant. In particular, the conversion of magnetic energy into thermal energy may elevate the pressure within the jet, increasing its internal energy budget. This effect could potentially manifest on larger scales when an important fraction of magnetic energy across the jet is involved, potentially altering the physical characteristics of the bulk flow and enhancing the role of thermal acceleration, due to the increased internal energy. Although investigating the role of resistivity and magnetic reconnection would be of great interest (see, e.g., Medina-Torrejón 2021; Mattia et al. 2023), it is beyond the scope of this work.

2.3.1. Nonforce-free configuration

In the nonforce-free configuration, the azimuthal magnetic field component is expressed as (Martí 2015b):

$$\begin{array}{c}\hfill B^{\varphi }(r)=\{\begin{array}{c}{\displaystyle \frac{2B_{j,\mathrm{m}}^{\varphi }(r/R_{\mathrm{m}})}{1+{(r/R_{\mathrm{m}})}^2},0\le r\le 1}\hfill \\ 0,r>1,\hfill \end{array}\end{array}$$(7)

in which $B_{j,\text{m}}^{\varphi }$ is the maximum toroidal field strength reached at the distance R_m = 0.37. In this configuration, the azimuthal magnetic field grows linearly in the region r ≪ R_m, reaches its maximum at r = R_m, and then decreases as ∝r⁻² when r ≫ R_m as shown in the upper left panel of Fig. 1. The axial component is assumed to be constant across the jet radius,

$$\begin{array}{c}\hfill B^z(r)=\{\begin{array}{c}{B}_j^z,0\le r\le 1\hfill \\ 0,r>1.\hfill \end{array}\end{array}$$(8)

Fig. 1. Initial profiles along the radial distance in the different magnetic field configurations. We use the representative values of Bϕ = 0.24, Rm = 0.37, vj = 0.8, pa = 0.01, and Bz = 0.05 (for the nonforce-free configuration) in code units. Upper panels: nonforce-free configuration. Lower panel: force-free configurations. Left panels: Magnetic field distribution of the azimuthal (blue line) and axial (green line) components. Right panels: profile distribution of the gas pressure (blue line), magnetic pressure (green line), and total pressure (black line). The dotted black line indicates the shear layer with m = 8 normalized at the higher value in the plot for an easier representation.

Therefore, the mean¹ azimuthal magnetic field component within the jet is:

$$\begin{array}{c}\hfill \overline{B_j^{\varphi }}=4B_{j,\mathrm{m}}^{\varphi }{R}_{\mathrm{m}}[1-R_{\mathrm{m}}\mathrm{atan}(\frac{1}{R_{\mathrm{m}}})]\end{array}$$(9)

while the axial one is simply $\overline{B_j^z}=B_j^z$.

By solving Eq. (1) with these magnetic field components, we obtain a gas pressure profile expressed as

$$\begin{array}{c}\hfill p(r)=\{\begin{array}{c}{\displaystyle 2[\frac{B_{j,\mathrm{m}}^{\varphi }}{{\gamma }_j(1+{(r/R_{\mathrm{m}})}^2)}{]}^2+C,0\le r\le 1}\hfill \\ \hfill \\ p_a,r>1.\hfill \end{array}\end{array}$$(10)

in which the integration constant can be computed using the boundary condition at the jet/ambient medium interface p^*(r = 1) = p_a, leading to

$$\begin{array}{c}\hfill C=p_a-\frac{{(B_j^z)}^2}{2}-\frac{{(B_1^{\varphi })}^2}{2{\gamma }_j^2}[1+{(R_{\mathrm{m}})}^2]\end{array}$$(11)

(in this expression, $B_1^{\varphi }=B^{\varphi }(r=1)$). The mean gas pressure is $\overline{p_j}=p_a-{(B_j^z)}^2/2$ and the magnetic pressure $\overline{p_{\mathrm{m}}}={\overline{b}}_j^2/2$. Fiducial magnetic and thermal pressure profiles are shown in the lower left panel of Fig. 1.

The magnetic tension (τ_m = −(b^ϕ)²/r, see, e.g., Martí et al. 2016; Moya-Torregrosa et al. 2021) is:

$$\begin{array}{c}\hfill {\tau }_{\mathrm{m}}(r)=-\frac{4{(B_{j,\mathrm{m}}^{\varphi })}^2}{{\gamma }^2{R}_{\mathrm{m}}^2}\frac{r}{{(1+{(r/R_{\mathrm{m}})}^2)}^2}.\end{array}$$(12)

2.3.2. Force-free configuration

In the force-free configuration, the azimuthal magnetic field component has the same profile as in the nonforce-free configuration, while the axial component becomes:

$$\begin{array}{c}\hfill B^z(r)=\{\begin{array}{c}{\displaystyle \frac{2B_{j,\mathrm{m}}^{\varphi }/{\gamma }_j}{1+{(r/R_{\mathrm{m}})}^2},0\le r\le 1}\hfill \\ 0,r>1.\hfill \end{array}\end{array}$$(13)

The axial magnetic field has its maximum at the jet center and decreases towards the outer radius as ∝r⁻². Its mean value is:

$$\begin{array}{c}\hfill \overline{B_j^z}=\frac{2B_{j,\mathrm{m}}^{\varphi }{R}_{\mathrm{m}}^2}{{\gamma }_j}[\ln (R_{\mathrm{m}}^2+1)-\ln (R_{\mathrm{m}}^2)].\end{array}$$(14)

In this configuration, the gas pressure is constant in the jet, with $p(r)=\overline{p_j}=p_a$, 0 ≤ r ≤ 1, and the pressure equilibrium is achieved by compensating the magnetic tension produced by the toroidal component of the field (see Eq. (12)) with the magnetic pressure produced by the axial one.

The advantage of this magnetic configuration is that the gas pressure becomes independent from the jet magnetization, allowing the field to increase without generating nonphysical, negative thermal pressures at the jet/ambient boundary (Moya-Torregrosa et al. 2021). This configuration is thus necessary to simulate models of highly magnetized jets in pressure equilibrium. The resulting initial radial magnetic field and pressure profiles are shown in the bottom panels of Fig. 1.

2.4. Jet and ambient parameters

With the aim to study jet acceleration from subrelativistic to relativistic outflow velocities, we inject the jets with the lowest possible axial velocity. The injection of low velocity, highly magnetized flows makes them submagnetosonic, which allows waves to propagate upstream and generate inconsistencies between the jet in the grid and its injection boundary condition. In particular, the axial magnetic field is changed and a jump appears between the boundary and the first cell in the grid, which represents an unphysical situation.

To overcome this problem, the initial jet velocity v^z was adjusted in the different models to result in magnetosonic Mach numbers of at least 1. While this approach prevents us from exploring the acceleration starting from subrelativistic velocities, it still allows us to achieve our main goal, which is understanding the role played by internal and magnetic energy fluxes in accelerating jets. Specifically, FR I jets can reach velocities up to γ ∼ 10 (e.g., Hovatta et al. 2009) at parsec scales, and from our injection velocities, which imply γ ∼ 1 − 2, we can explore the increase from mildly relativistic to relativistic speeds for different types of jets, namely magnetically, internally or kinetically dominated, and the velocity structures that they generate.

Using our information on NGC 315 as starting point, we design our grids to simulate outflows along the z-coordinate at z_i = 0.3 pc, since this corresponds to a distance at which the jet has already velocities in the range $v_j^z~0.3c$ to $v_j^z~0.8c$, depending on the frequencies and epochs considered (see Fig. 3, Ricci et al. 2022). The grid extends up to z_o = 3.3 pc, covering a distance that is large enough to let us explore the acceleration phenomenon from subparsec to parsec scales. In the radial direction, the grid extends from r_i = 0 up to r_o = 0.9 pc, with an initial jet radius of r_j = 0.03 pc (see Fig. 3, Boccardi et al. 2021). Taking into account that jets are built with large overpressure factors (see below), this radial size of the grid is large enough to allow for jet expansion. With our choice of the initial jet radius, our numerical domain consists of a grid of 30r_j × 100r_j with n_r × n_z = 600 × 2000 cells.

To determine the initial magnetic field we refer to Fig. 2, which highlights the possible magnetic field values, in our region of interest, extracted from the results published in Ricci et al. (2022). According to the different models and different magnetic field evolution, the initial total strength is allowed to span between B = 0.02 G and B = 0.7 G. The magnetic field structure considered in this paper ranges from helical in the force-free configurations (with the mean pitch angle defined as ${\overline{\varphi }}_{B_j}={\tan }^{-1}({\overline{B}}_j^{\varphi }/{\overline{B}}_j^z$)), to a dominating toroidal field in the nonforce-free configurations.

Fig. 2. Expected magnetic field distribution between 0.1 pc and 3 pc. The upper and lower magnetic field limits are extrapolated from the core shift measurement under the different assumption (for the details see Ricci et al. 2022). The black vertical line highlights the starting point of zi = 0.3 pc of our simulations, while the horizontal dashed lines are the expected range of magnetic field at zi.

The numerical grid is initially filled with a jet with the injected properties (see Table 1) in the region (r, z)∈[0, 1r_j]×[0, 100r_j] and an ambient medium that occupies the rest of it. The initial jet overpressure forces its expansion into this ambient medium, which we let evolve until an equilibrium, steady solution is reached. The description of the criteria by which we decide when equilibrium has been reached are described in Sect. 3.1.

By establishing a jet across the whole grid, we avoid jet injection and the generation of a bow shock that would arise from the interaction of the newly ejected jet with the external medium, while focusing on jet physics alone.

We assume jets with kinetic power in the range L_j ∼ 10⁴³ − 10⁴⁴ erg s⁻¹, in accordance with the estimate given by Morganti et al. (2009), Ricci et al. (2022) for NGC 315. To build the required jet powers, we distribute the energy among the different channels (see Sect. 2.3). Regarding the external medium, we refer to the X-ray Chandra observations for the environment surrounding NGC 315. The measurements for the galactic core region (the inner 0.3 kpc) give pressure and density of p₀ = 4.5 × 10⁻¹⁰ dyn cm⁻² and ρ₀ = 0.46 × 10⁻²⁴ g cm⁻³, respectively (Worrall et al. 2007). However, jet configuration for the aforementioned powers requires pressures that can be six to eight orders of magnitude over these values.

While this raises questions about jet equilibrium and collimation at these scales, it is well possible that the jet is indeed strongly overpressured with respect to its environment, but that this extreme overpressure has been slowly achieved with time. It is relevant to recall that when formed and injected for the first time, jets are surrounded by their own shocked plasma, which is probably overpressured with respect to the jet itself. As it expands, the pressure around the jet falls gradually (see, e.g. Perucho et al. 2014, 2022; Perucho 2019) and the jet keeps its collimation due to its own velocity, which limits opening angles to ∼1/γ. In conclusion, in the case of evolved jets (like those in radio galaxies) one could expect free expansion with large opening angles while the jet accelerates, asymptotically reaching ∼1/γ. This is in agreement with the transition from parabolic to conical jet expansion, precisely in the acceleration region.

Unfortunately, in our case, we cannot let the jet to adapt to such extreme conditions in a gradual way, and we establish the overpressure since the beginning of our runs. Furthermore, such large overpressure causes the code to crash as the jet violently expands into the ambient medium. However, since we only need the jet to freely expand, it is enough to give the ambient medium a pressure that allows this, even if the simulated ambient pressure is much higher than that estimated by Chandra. Therefore, we assume arbitrarily initial high pressure and density, with values of p_a = 1.0 × 10⁻³ dyn cm⁻² and ρ_a = 1.67 × 10⁻²² g cm⁻³, respectively. This medium is set as isothermal, and the pressure (and density) are given a gradient in the z direction to fall with distance as ∝z⁻¹. In this way, at injection the jet overpressure factor spans from 1 (i.e., jet in thermal pressure equilibrium with the environment) to 4 and reaches values up to 10² at the end of the grid. The initial jet density spans between ρ_j = 1.67 × 10⁻²⁵ − 1.67 × 10⁻²³ g cm⁻³.

Jet pressure and density, together with the magnetic field strength, are modified from model to model in order to explore a variety of initial jet conditions, from highly magnetized, cold jets (with a magnetic energy flux F_M, representing ∼75% of the total) to hot jets (with an internal energy flux, F_I, representing ∼88% of the total). The parameters that determine each of the simulated models are given in Table 1. The model names are given using the following code: the first letter represents the magnetic field configuration, with F for force-free and N for nonforce-free, the second letter stands for the type of jet, namely magnetic (M), internal (I) or kinetic (K) energy dominated, the third letter gives the code for the jet density, with S for ρ_j ≥ 0.1 ρ_a, H for 0.01 ρ_a ≤ ρ_j < 0.1 ρ_a, and L for ρ_j < 0.01 ρ_a, while the last numbers refer to the jet thermal overpressure factor at injection. When two different models are expected to have the same nomenclature, we add at the end of the name a letter in progressively alphabetical order (e.g. FMH1b, FMH1c). Models with _m4 have the shear layer exponent set to m = 4, a wider (thicker) shear layer, while for all the other models the exponent is set to m = 8.

Table 2 displays the ratio of the different flux channels over the total one (index i) as derived from the jet parameters, together with the simulated ones, which are convolved to the shear layer (index s), and final ones computed at the final row of cells at 3.3 pc (index f). In Appendix A, we discuss the reason and implications in the initial radial profile of the fluxes within the jet, considering their convolution with the shear layer.

3. Results

3.1. From initial conditions to equilibrium

The simulation setup consists of a jet that extends the injection conditions along the numerical grid with constant radius r  =  r_j, with a parameter smoothing function out towards the ambient medium, a shear layer, as described in Sect. 2.4. The jet structure is computed considering pressure equilibrium with its environment.

Jets immediately start their expansion into the under-pressured ambient, transferring energy via radial propagating waves. These waves accelerate the ambient medium in the radial direction. As a consequence, some material is pushed out of the grid, slightly changing the conditions in the ambient. We have run tests with an extended grid to check the influence that this loss may have and found that the jet is unaffected, probably due to the large jet overpressure, which does not change significantly. For proofs in support of this and for a proper description of the pressure wave and its evolution, see Appendix B.

We run the simulations for the necessary amount of time for the jets to reach equilibrium after the expansion is completed. To determine when the steady solution is reached we use the radial velocity as the indicator. Specifically, we consider equilibrium to have been reached in a certain model when |v^r|< 0.01 all over the grid, with the exception of jet expansion/recollimation associated with the achieved equilibrium configuration. An example of this is shown in Fig. 3 for model FMH1_m4. In the figure, the white contours highlight the levels of the tracer, that is the jet mass fraction. In this case, outside of the jet, the grid shows radial velocities lower than 0.01 everywhere.

Fig. 3. Absolute value of radial velocity for the model FMH1_m4. The white contours highlight the tracer values at levels of 0.2, 0.4, 0.6, and 0.8 (the tracer levels are in decreasing order, i.e., from 0.8 to 0.2, from left to right). Outside of the jet structure delineated by the shear layer (white contours), the velocity is smaller than 0.01 everywhere in the grid, meaning that this can be considered a steady solution.

Figure 4 shows an example of radial equilibrium profiles at different locations along the jet and for different resolutions, for model FMH1. We show the radial distribution of different physical quantities at t = 0 and at the final step, at three different distances, injection (cell 1), half grid (cell n_z/2), and z = z_max (cell n_z). The figure shows three sets of panels for the results at three different numerical resolutions: the upper set shows the results for a simulation with 10 cells/r_j, the middle set shows the case with 20 cells/r_j, while the bottom set shows the case with 40 cells/r_j. For the latter case we show only the initial and final cell profiles, since we used a smaller grid of [30r_j × 30r_j] for computational time reasons. The bottom right plot of the 10 cells/r_j case shows a change of the axial field at injection, which implies a discontinuity with respect to the boundary condition. This problem, caused by numerical diffusion, disappears at 20 and also at 40 cells/r_j. Having found convergence between our results at 20 cells/r_j, we adopt this resolution for all our simulations. In Appendix C we show further proofs in support of this choice.

Fig. 4. Radial profiles of the different physical properties for model FMH1 at different distances from the injection (cell 1 in red, half grid in blue, and final cell in green) and at different times (initial in dotted lines and final in continuous lines). In the top panels, we show the radial profiles with a resolution of 10 points, in the middle the ones at a resolution of 20 points, and in the bottom ones the profiles with 40 points (here, the half grid is not shown due to the smaller simulated grid). The convergence in the profiles between 20 and 40 points led us to use resolutions of 20 points for our simulations.

The plots show how the initial condition is preserved at injection, but the jet expands significantly, becoming much wider along the grid. These plots also show an increase in axial velocity, which appears to be displaced from the axis. We discuss this result in the following section.

3.2. Jet acceleration and expansion profiles

Figures 5 and 6 show the results obtained for all the different models. A clear relation between jet velocity profile and its intrinsic properties is inferred: jets developing strong internal shocks -Mach disks- present stronger acceleration in their outer, expanding layers, whereas those developing milder, conical shocks present a fast spine. Pressure maps (Fig. 6) show that the strong discontinuities observed in several models such as FMH1, FIH4, and FIH3b, are caused by the generation of a Mach disk in the jet. The Mach disk formation is driven by the strong expansion, due to the large angles formed by the shock waves with the jet axis (see, e.g., Martí et al. 2018, and references therein). In general, the main difference between the jets showing acceleration at the outer layers and those showing spine acceleration is their force-free versus nonforce-free nature, or magnetically versus internal energy dominated.

Fig. 5. Lorentz factor maps for the simulated models. The white contours are representative of the tracer at levels of 0.2, 0.4, 0.6, 0.8. The results show the clear correlation between the the jet speed profiles and the intrinsic properties, with the main differences shown between force-free and nonforce-free along with magnetically dominated or internally dominated models. For details, see the text.

Fig. 6. Pressure maps for the simulated models. The white contours are representative of the tracer at levels of 0.2, 0.4, 0.6, 0.8. As in Fig. 5 the different jet nature leads to the different pressure profiles in the spine and in the shear layer regions.

However, we observe that an increase in the width of the shear layer (FMH1_m4 vs. FMH1 models) can completely change the jet acceleration pattern obtained. In the case of FMH1 (top, left panel), the jet expands until the point at which the extension of the planar shock wave at the center of the jet reaches its surface. Beyond this point, we observe a mild expansion/recollimation process. Downstream, in the region that is unaffected by the Mach disk, the jet develops instability patterns and locally dissipates some of the gained kinetic energy, which explains why the internal energy flux at the end of the grid is larger than at ∼0.5 pc (see Sect. 4.3 for the discussion on this effect). In FMH1_m4 (top, second panel from the left), on the contrary, a Mach disk is not formed and the acceleration of the flow is concentrated towards its axis, unlike all the other force-free models. The corresponding pressure maps show that the main difference between both models is the delay in the wave propagation that the thicker shear layer causes in FMH1_m4 with respect to FMH1: the outgoing waves have to cross a thicker portion of the jet before bouncing back, which delays the bounce and moves it downstream. This delay relaxes the angle formed by the incident wave and the axis, changing the acceleration pattern and jet structure. Figure 7 shows the modulus of the gradient of the rest-mass density in the (r, z) plane², |∇ρ|, of both FMH1 and FMH1_m4 (top panels), and a zoom of the relevant region where the first recollimation shock forms (bottom panels). The modulus of the gradient of the rest-mass density highlights the location of steep sound waves and shock waves. The images in Fig. 7 show how the different widths of the shear layer changes the wave structure of the models. In model FMH1, the jet dynamics is dominated by instabilities both in the spine and at the jet/ambient medium transition on parsec scales, where we see the formation of vortexes typical of the growth of Kelvin-Helmholtz instabilities. On the contrary, in model FMH1_m4 the jet pinches, but substantially fewer instabilities arise and the jet/ambient medium transition is smoother. In the bottom panels, we can see that the opening angle of the transition between the inner jet and the shear layer is larger in FMH1. The bottom panels display contours for tracer f = 0.9 to show the transition between the inner jet and the shear layer, which, at injection, occurs at 0.91r_j in model FMH1, while it is at 0.82r_j in model FMH1_m4. We propose that this difference determines the two possible evolution patterns observed (see Sect. 4.3 for details). We have run other models with wider shear layers (with the same initial conditions as FIH2, FIH4, FMH1c, FIH3b) and we have observed the same behavior: thicker shear layers may result in a dramatic change in the acceleration pattern of force-free jets, avoiding the formation of Mach disks and favoring spine acceleration, as opposed to thin-layered models (see Appendix D).

Fig. 7. Modulus of the gradient of the rest-mass density (in logarithmic scale) for models FMH1 (left panels) and FMH1_m4 (right panels). The upper panels show the entire grid, while in the bottom ones we show a zoom-in in the 0.3–0.9 pc region. The modulus of the gradient of the rest-mass density highlights the location of steep sound waves and shock waves. The black squares highlight the triple point associated with the incident/reflected/Mach shocks in the jet-ambient transition surface (for details, see Sect. 4.3). In the bottom panels, the contours represent the tracer at the level f = 0.9, showing the transition between the inner jet and the shear layer.

In all the force-free models (except FMH1_m4), we observe the formation of the Mach disk with the resulting transversal velocity structure consisting of a slow spine surrounded by a fast sheath. In all cases, although this is better observed in model FMH1, the accelerated region slowly expands towards the jet axis, even if at the end of the grid none of them shows a completely accelerated jet cross-section. Instability modes develop at the interface between the spine and the outer accelerated layer. This triggers dissipation and we thus expect it to contribute to increasing the internal energy of the region. The fast development of instabilities is expected in the case of hot, slow flows (see, e.g., Perucho et al. 2004, 2005).

Nonforce-free models show, in contrast, an acceleration pattern concentrated around the jet axis, meaning that they produce a fast spine. The main difference with respect to force-free models is the absence of a Mach disk. This acceleration pattern corresponds to that expected from the Bernoulli process. Highly overpressured jets, such as NIH4, show strong acceleration on the axis before the flow reaches a strong recollimation shock, where the flow is decelerated. The next expansion process produces a second acceleration region, but the Lorentz factor reached is 20% smaller than the one prior to the first recollimation shock. In the case of NIH1, in thermal pressure equilibrium with the environment at injection, acceleration is milder, modulated by recollimation waves/shocks, with a progressively growing velocity from one acceleration region to the next.

Figure 8 shows the radially-averaged Lorentz factor for all the simulated models. For clarity, we divided the profiles into two separate plots: in the top panel we show the magnetically dominated models, while in the bottom panel the thermally dominated ones together with the kinetically dominated FKS1 are displayed. The average Lorentz factor value is computed by weighting the value at each cell with its volume (given by a ring). As a discriminant to determine the jet width we use the tracer parameter with a value of 0.5 as a threshold, implying that cells with tracer > 0.5 are considered to be part of the jet while cells with tracer < 0.5 are considered ambient medium. Along with the speed profiles, Fig. 8 shows the data points for NGC 315 in the frequency range 5–43 GHz and at different epochs inferred in Ricci et al. (2022). These results were obtained based on the observed jet-to-counterjet intensity ratio, which provides an average velocity in a likely stratified jet structure (hints of an edge-brightening and possible stratified velocity are given in Park et al. 2021). The comparison of the simulation profiles with the observational data is presented and discussed in Sect. 4.4. The plot shows that the simulations with different setups result in average jet acceleration from γ ∼ 1 − 2 up to γ ∼ 2 − 5 within a spatial propagation of only 3 pc. The weakest acceleration is observed for model FKS1, as expected for kinetically dominated models, whereas the strongest acceleration is obtained in models with the highest relative magnetic and/or internal fluxes at injection, such as FMH1c, FIH4 or FML1. It is also interesting to mention that force-free models accelerate larger volumes of jet plasma, as the outer regions of the jet represent larger rings in cylindrical coordinates. All models show modulations in the average Lorentz factor produced, as stated above, by expansion and recollimation episodes, in which internal energy flux is invested into acceleration and the jet flow is then decelerated by shocks, where kinetic energy is converted into internal energy. As expected, this is more evident in significantly overpressured models, such as FIH4, where the initial acceleration of the jet is followed by a drop in the average Lorentz factor. In the thin shear-layer, force-free jets, the difference in velocity between the plasma flowing along the spine, decelerated by the planar (Mach) shock, and that accelerated at outer radii triggers the development of instabilities. In models FMH1, FMH1.3, and FMH0.3 these instabilities develop to nonlinear amplitude within the grid and, on the one hand, allow the gained axial momentum at the accelerated region to be shared with the inner jet spine and, on the other hand, dissipate part of the gained kinetic energy. Altogether, this results in a drop of the maximum average Lorentz factor achieved by the jet during its initial expansion, and a slow recovery with distance. Nevertheless, these models stabilize their Lorentz factors at values between 2 and 3.

Fig. 8. Radially averaged Lorentz factor as a function of distance for all the simulated models. Top panel: magnetically dominated jet models. Bottom panel: thermally and kinetically dominated jet models. The different data points are the Lorentz factor for different epochs and frequencies of NGC 315 inferred by Ricci et al. (2022). The continuous lines represent the force-free magnetically dominated models, the dashed lines the force-free thermally dominated ones, the dotted lines are for the kinetically dominated model, and the dashed-dotted lines for the nonforce-free configurations.

Figure 9 shows the half-opening angle of the different models versus distance together with the observational points for NGC 315. We divided the plot as before, the upper panel for the magnetically dominated models, and the bottom panel for the thermally dominated (plus FKS1) ones. In all models, we observe an initial increase of the opening angle due to the expansion of the jet towards the less-pressured ambient medium, followed by a decrease that starts at different distances, but always before ∼1 pc. The fall in the opening angle profiles is associated with jet collimation, r ∝ z^α, with α < 1. None of the models shows completed collimation on the simulated scales (α ≈ 1). All models tend to final half-opening angles of 1 − 3°. This plot is very clarifying regarding the causes for jet expansion in the different models: force-free models expand faster than nonforce-free models, as can be seen by comparing cases FIH4 (green dashed line) and NIH4 (blue dashed line). This difference also allows us to group the models into two main behavior patterns: 1) the models with the largest initial opening angles develop the Mach disk, while 2) those with smaller initial opening angles develop conical shocks.

Fig. 9. Half-opening angle as a function of distance for all the simulated models. Top panel: magnetically dominated jet models. Bottom panel: thermally and kinetically dominated jet models. The force-free magnetically dominated setups are shown in continuous lines, the force-free thermally dominated in dashed lines, the kinetically dominated in dotted lines, and the nonforce-free in dashed-dotted lines. The data points represent the intrinsic half-opening angle profiles of NGC 315 obtained by Ricci et al. (2022) for multi-frequency and multi-epoch VLBI observations.

The opening angle at injection is also clearly correlated with initial jet overpressure, so force-free models with no thermal contribution to overpressure at injection (p_j = p_a) show smaller opening angles, albeit longer expansion regions (FMH1b, FMH1c and FML1). These models also show the largest average Lorentz factors (see Fig. 5), mainly driven by magnetic energy (as can be seen in Table 2 and Fig. 10). Another relevant aspect to consider is that these models in (thermal) pressure equilibrium at injection also have larger initial velocities (see Table 1), which causes smaller opening angles and therefore longer expansion lengths. Figure 9 also shows the half-opening angles derived from the same observational results (Ricci et al. 2022). For the discussion of such comparison, we refer to Sect. 4.4.

Fig. 10. Evolution of the integrated energy fluxes along the axial direction for all the models reported in Table 1. Black, red, green, blue, and orange lines represent, respectively, the total, kinetic, internal, magnetic, and rest-mass energy fluxes. Depending on the dominating energy channel, its dissipation in favor of the jet kinetic flux is visible in all the models. An exception is the kinetically dominated model FKS1, in which both the thermal and magnetical channels are dissipated to accelerate the jet, even if not dominant.

Figure 10 shows the evolution of the different energy flux channels, together with the total one, along the z direction. These plots clearly show the correlation between the increase in kinetic energy flux and the drops in either magnetic or internal energy fluxes, or both. From the plots, it becomes evident that, at injection (z  <  0.6 pc), a clear anti-correlation between internal energy and kinetic energy fluxes is only seen in internal energy dominated (FI and NI models). In the case of magnetically dominated models, it is the magnetic energy flux that shows anti-correlation with the kinetic energy flux, and we can observe regions in which the kinetic and internal energy fluxes are correlated (e.g., FMH1b). We therefore deduce that the dominant acceleration mechanism is controlled by the dominating energy flux, but that, in specific conditions, for example, in expansions for the case of internal energy, the secondary energy flux can also contribute to acceleration. FKS1 represents an interesting case (Fig. 10, third row, central panel), in which the initial expansion results in a small increase of the kinetic energy flux favored by a drop in both internal and magnetic energy fluxes (although this jet reaches the smallest terminal values of the Lorentz factor; see Fig. 5).

Finally, a slight global increase in magnetic flux is observed in the FMH1_m4 model, beyond a distance of z = 2.5 pc. The latter is in correspondence with the jet recollimating after an expansion phase (see Fig. 5). During the recollimation, the magnetic field is compressed and its value increases, but an increase in magnetic flux is typically not expected unless there is a boost in velocity due to internal energy (Bernoulli processes). However, this is not the case in the indicated region, so we propose such a small increase to be due to i) a physical process that transfers kinetic energy to magnetic energy as the reverse of magnetic acceleration, likely a consequence of local differential recollimation effects (a similar behavior is observed in Komissarov et al. 2015, Fig. 10, right panel); ii) small inaccuracies arising during the summing of the different integrated quantities, or iii) numerical diffusion at the jet boundary, which is the magnetic field smearing out towards the outer cells. The exploration of such an effect is out of the scope of this paper and will be tackled by future works. Nevertheless, we remark that the fraction of such an increase is negligible with respect to the total jet and does not affect our results.

4. Discussion

Our results show that jet acceleration is driven by the dominant energy flux at injection, which can be either magnetic or internal, but both are successful in increasing the jet Lorentz factor in expanding regions. However, there are significant differences in the acceleration patterns in both cases, as shown in Sect. 3.2.

In Sect. 4.1, we discuss the role of thermal acceleration in relativistic jets. The differences between the velocity structures due to the force-free versus nonforce-free and thin versus thick shear layers are discussed in Sects. 4.2 and 4.3, respectively. In Sects. 4.4 and 4.5, we focus on the comparison with observational data and derive conclusions as to how our results can provide insights for jet evolution. Finally, in Sect. 4.6, we comment on the limitations of our study.

4.1. The role of the thermal acceleration

Jets produced by the Blandford-Znajek mechanism (Blandford & Znajek 1977) are probably dilute, hot, and significantly magnetized at subparsec scales, where entrainment has still not affected jet composition. Building upon the principle of Bernoulli, thermal acceleration is effective when the enthalpy h is high enough (specifically, h ≫ c²). For a steady, relativistic flow, the principle is expressed by the law hγ = const, establishing that the flow accelerates when h decreases (as the flow expands). The process of thermal acceleration has been regularly seen in numerical simulations (e.g., Perucho & Martí 2007; Anglés-Castillo et al. 2021). By assuming that jets are hot at injection, we also found evidence for thermal acceleration at these scales: as highlighted in Sect. 3.2 and shown in Fig. 10, we obtain the expected anti-correlation between thermal and kinetic energy fluxes during the evolution of our models. For instance, in the case of FIH3b model, in which the magnetic flux is an order of magnitude lower than the internal one and irrelevant to the acceleration of the jet, the jet is clearly accelerated by the Bernoulli process alone. Interestingly, such model, and so thermal acceleration alone, is able to reach Lorentz factors higher than 2 on scales of ∼2 − 3 pc, matching thus the observed speed profile of NGC 315 (see Fig. 8).

Thermal acceleration can be relevant in relativistic jets at larger scales than those assigned to magnetic acceleration: if jets keep or gain internal energy on parsec scales, as it can be the case at recollimation shocks, it can act at the following conical expansions. In other words, even when jets are magnetically dominated at injection, once they are collimated and the magnetic acceleration stops, internal energy could still nourish the bulk acceleration. The detection of jet flow acceleration on large scales, for example, ∼100 pc in blazar jets (Homan et al. 2015), could be thus interpreted as an observational signature of thermal acceleration.

4.2. Force-free versus nonforce-free models

One of the main triggers of the different acceleration patterns in our models is the force-free versus nonforce-free nature of the flows. The reason lies in the cancellation of the field force in force-free models, eliminating the contribution of magnetic tension to jet collimation. Indeed, in the force-free configurations, the magnetic tension does not act as a collimating factor, and expansion is therefore enhanced, with half-opening angles reaching 7°. On the contrary, in the nonforce-free case, magnetic tension plays an active role in controlling expansion, which can only be temporarily overcome by considerable overpressure (as in, e.g., model NIH4). It is noticeable that the recollimation pattern in the nonforce-free model NIH4 is, as a result, very different from that in FIH4, the force-free case: the angle formed by the shock wave and the jet axis is clearly smaller than in FIH4, and therefore NIH4 does not develop a Mach disk.

Because acceleration in the outer layers only takes place in force-free jets, which are, in general, magnetically dominated, we could derive the conclusion that this type of velocity pattern (slow spine and fast sheath) is a signature of magnetic acceleration in the case of thin shear layers, although FI models show that this does not strictly exclude internal energy as an accelerating mechanism. Fast spines would correspond, in contrast, either to internal energy-dominated flows or strongly sheared (possibly wind-shielded, two-flow) jets.

4.3. The role of shear layers

Our simulations show that the different thicknesses of the shear layer can lead to a fundamental change in the acceleration pattern of force-free jets. Figure 7 shows that this may be caused by the difference in the opening angles of the region separating the inner jet and the shear layers (indicated by a white contour that stands for tracer f = 0.9). Initially, thicker layers produce smaller opening angles and limit recollimation to conical shocks, allowing acceleration along the spine, even if modulated by the succession of expansions and recollimations. We propose that the direct consequence of the thicker shear layer is the delaying of the wave propagation. As a result, the triple point, associated with the incident/reflected/Mach shocks in the jet-ambient interaction surface, which is visible in Fig. 7 (highlighted by the black squares) is found in model FMH1 at ∼0.7 pc while in model FMH1_m4 around 1.3 pc. The role of the width of the shear layer has been confirmed by further simulations (see Appendix D).

Comparing the maximum Lorentz factors reached by the different models (Fig. 5), we observe that FMH1_m4 reaches a slightly larger Lorentz factor than its thin-layered counterpart, FMH1. The most plausible explanation of this result is that the growth of the instabilities observed in model FMH1 dissipates some fraction of the kinetic energy of the flow, whereas the jet of model FMH1_m4 is shielded from its environment by a thick shear layer, which, together with further expansion can contribute to preserving its collimation (see, e.g., Perucho 2019).

Altogether, we can conclude that the role of shear layers at the collimation/acceleration region can be of great importance for jet evolution, with thick-layered jets producing collimated, stable, fast flows with velocity profiles consisting of a fast spine and a slower, shielding layer. On the contrary, thin-layered jets initially develop a slow, hot spine, surrounded by a sheath along which the flow is accelerated. The boundary between these two regions, spine and sheath, is prone to the development of instabilities and can eventually dissipate part of the gained kinetic energy.

We highlight, however, that this is a small-scale effect (see Sect. 3.2, e.g., model FMH1) and is expected to not be largely significant (see the small velocity differences between FMH1 and FMH1_m4, as aforementioned) within the simulated grid because the instabilities have not reached the nonlinear regime and have not developed mixing between the hot and slow spines and the faster and cold sheath.

4.4. Comparison with NGC 315

A one-to-one comparison between observations and simulations is, in general, not possible due to limitations of both. Radio maps can be affected by uncertainties in the calibration and imaging processes, as well as by time-variable and local effects that do not necessarily allow us to infer the physical conditions of the underlying plasma flow. In simulations, as discussed, some simplifications are unavoidable. Nonetheless, some general trends can be identified and discussed qualitatively, and we do this in the following by considering NGC 315.

Concerning the jet opening angles, Fig. 9 shows that most of our simulated models succeed in reproducing the terminal values at the end of the collimation region, as inferred from the 5 GHz observations at different epochs (θ ∼ 1 − 3°). These opening angles also agree with the average values observed in radio galaxies (e.g., Pushkarev et al. 2017). In contrast, we observe that the opening angle in the innermost jet, described at the highest radio frequencies (15, 22, and 43 GHz), does not match our synthetic profiles. The limitations suffered by numerical models forced us to inject super magnetosonic flows, as previously explained (Sect. 2.4). This also makes any comparison challenging at the compact scales: while the simulated jets open abruptly from the injection point at ∼0.3 pc due to initial overpressure, the high-frequency data from the jets in NGC 315 indicate that such an opening happens on smaller scales and that the jet at 0.3 pc is going through the final part of its collimation before transiting to the conical shape (Boccardi et al. 2021). Nonetheless, as mentioned, the observational and simulated data reconcile on parsec scales.

Except for the innermost jet region, an agreement between observations and simulations can be found in the fact that they both show a co-spatiality of the acceleration and collimation processes. The simulations show acceleration episodes between (0.5 − 1.0) pc and (1.5 − 2.0) pc (Fig. 8) along the region where the opening angle is also decreasing (Fig. 9). This is probably the last part of the accelerating region, which in NGC 315 extends within the inner ∼1 − 2 pc. The co-spatiality agrees with the expectation, for FM models, of a magnetic acceleration being relevant when the jet is evolving with a parabolic or quasi-parabolic shape (Komissarov et al. 2007). Nonetheless, as mentioned in Sect. 4.1, thermal acceleration alone could in principle reproduce the observed velocities on parsec scales.

A comparison between the simulated and observed Lorentz factors may indicate a better agreement for models with jet power ∼10⁴³ erg s⁻¹, consistent with the jet power estimated from observational constraints (see, e.g., Morganti et al. 2009; Ricci et al. 2022). Indeed, the average Lorentz factor of ∼2 on scales of ∼2 − 3 pc can be reproduced by models with such jet powers. On the contrary, models with higher jet power (∼10⁴⁴ erg s⁻¹) seems to fail to align with the observational data: while they manage to reach the Lorentz factor peak of γ ∼ 4 at the distance of ∼1.5 pc (observed in only one epoch), they exhibit terminal Lorentz factors exceeding the maximum observed values by more than a factor of two. Remarkable examples are the highly-magnetized models, such as FML1, FMH1c, and FMH0.3. As in the case of the opening angle, a meaningful comparison on subparsec scales is not possible: data at 22 GHz and 43 GHz show acceleration from very small velocities, out of the reach of our setups (see Sect. 2.4).

Another interesting aspect might be noticed from Fig. 8: the Lorentz factor evolution along the jet changes between epochs, even at the same frequency, for example, the 1996 5 GHz epochs, which are separated by five months. While taking into account that such variations can be a consequence of various local and observational effects, it is also possible that they are the result of intrinsic changes in the injection conditions. Indeed, the jet velocities and variability in Lorentz factors recovered from the family of models with jet powers ∼10⁴³ erg s⁻¹ could potentially match the observed time-dependent changes in the Lorentz factor. The simulations we have presented here differ from each other by only a few parameters, within a range of plausible values derived from observational results (see, e.g., Fig. 2). Thus, we could suggest that the different jet profiles inferred from the same frequencies but different epochs may arise from variations in the jet properties at the formation site, within timescales of a few months. This might show how jets are far from being subject to regular injection and that conditions and, consequently, acceleration patterns, can behave in an extremely dynamic way. Further evidence for time-variability of the physical properties at the jet base is provided by the observed time-dependence of the core shift effect (Plavin et al. 2019), as well as by the well-known flux density variability that characterizes the nuclear regions on timescales from hours to months.

4.5. Comparison with other AGN jets

As discussed, most of the force-free jets we have simulated are characterized by sheath acceleration and by a slow spine. This result appears at odds with observational constraints, which rather indicate the existence of a fast spine and a slow sheath. Such a velocity structure was suggested, for instance, to reconcile the observed spectral properties of FRI radio galaxies with their beamed parent population, the BL Lacs (Chiaberge et al. 2000), and appears required to explain the general anti-correlation between the observed Lorentz factor and the jet viewing angle (see, e.g., Homan et al. 2021, Fig. 12). Misaligned jets, which are seen from the side, are usually characterized by mildly relativistic speeds, as opposed to blazar jets, which are observed closer to the jet axis with Lorentz factors of the orders of tens (e.g., Lister et al. 2019). This is in agreement with the frequent observation of limb-brightening in misaligned jets (see, among others, the recent cases of Centaurus A and NGC 315 by Janssen et al. 2021; Park et al. 2021, respectively), which can be interpreted as due to the Doppler de-boosting of the fast spine. Note that this applies not only to the low-power FRI jets, but also to the powerful FRIIs (see the case of Cygnus A, presented by Boccardi et al. 2016).

Consequently, in this scenario, a match between our simulations and observations implies that jets would need to propagate with shear layers protecting the spine or with nonforce-free magnetic field configurations, since jets with fast spine are produced in both families of models. On the opposite side, while faster shear layers with hot spines are likely not suggested from VLBI observations, such a radial velocity structure is not ruled out. In this situation, the initial force-free, thin shear layer models are favored.

As recently shown by Boccardi et al. (2021), the properties of the outer jet sheath may vary depending on the accretion mode in the AGN. Specifically, an extended and possibly disk-launched jet sheath was observed in High Excitation Radio Galaxies (HERG), which are mostly characterized by FRII morphologies. On the other hand, Low Excitation Radio Galaxies (LERG), associated to both FRI and FRII morphologies, showed in comparison a narrow jet anchored in the very inner regions of the accretion flow. In future works, we aim to explore possible connections between HERG/LERG and the different acceleration patterns we observed in this paper.

Furthermore, in a future work, we plan to extrapolate synthetic spectral index maps from our models using the post-processing code published in Fromm et al. (2018). In this way, we will be able to compare the simulated spectral properties with the observations and explore whether the different acceleration patterns we recover, and especially the outer layer acceleration, may result in characteristic signatures that can be related to the underlying magnetic field properties (e.g., Kramer & MacDonald 2021).

4.6. A note of caution regarding our setups: Jet evolution

When jets are initially launched in AGN, they trigger a bow shock that propagates through the ambient medium. Because the jet flow advances faster within the formed jet channel than this forward bow shock, the jets get surrounded by overpressured, shocked jet plasma (known as backflow or cocoon). Therefore, collimation must be, in this initial phase, controlled by both the pressure of this shocked environment and the jet toroidal field. As the jet expands and the bow and reverse shocks propagate to large distances, the pressure in the cavity formed by the jet falls continuously with time (see, e.g., Perucho et al. 2014; Perucho 2019). Then, the gravitational potential of the galaxy, which acts in a dynamical time $\sim 1/\sqrt{G{\rho }_{\mathrm{DM}}}$ (where G is the gravitational constant and ρ_DM is the dark matter density in the galaxy), allows the galactic pressure and density profiles to be recovered (see, e.g., Perucho et al. 2014) and therefore, pressure should fall to the original ISM values of ∼10⁻¹⁰ dyn/cm².

The expected jet energy fluxes and magnetic field intensity in this region (F_j ∼ 10⁴³ − 10⁴⁴ erg s⁻¹ and B ∼ 0.1 − 1 G) result in extremely large overpressure of the jet with respect to the ISM. Once the jet is accelerated to relativistic velocities, the opening angle is proportional to the inverse of the Lorentz factor (see Komissarov 2012). The large initial overpressure poses a strong challenge to jet collimation and to numerical simulations aimed at studying this region. In our simulations, we decided to achieve overpressures as large as possible, which could allow free expansion, by introducing significant pressure gradients in distance. Once free expansion is allowed, the ambient pressure becomes relatively irrelevant, and this has allowed us to run our simulations. Establishing ambient media with realistic initial pressure causes the code to crash and makes such simulations impossible.

An obvious question arises: how is it then possible that the jets manage to settle in such ambient media with a pressure that can be orders of magnitude smaller than the jet pressure? We would hypothesize that it is because the jet has adapted to its environment in a continuous, relatively slow way during the time it takes for its terminal shock to reach kiloparsec scales, as opposed to simulations. In addition, this probably makes a nonforce-free magnetic field configuration necessary, for magnetic tension to play a significant role in jet collimation in the region where it is precisely still strong, namely the sub and trans-Alfvènic regions, until the jet is collimated by its own velocity. Furthermore, the collimating effect of the toroidal field also avoids the loss of internal energy via Bernoulli acceleration, allowing this energy budget to act at larger scales (see also Anglés-Castillo et al. 2021).

5. Conclusions

In this paper, we performed a numerical study of the jet acceleration in typical FR I radio galaxies using the well-known source NGC 315 as a starting model. We explored the evolution on subparsec and parsec scales of different axisymmetric jet models using a two-dimensional RMHD code, and compared our final results with the observational properties inferred on the same scales by means of VLBI observations. Our results are summarized here.

• All the simulated models show acceleration from initial Lorentz factor values of γ ∼ (1 − 2) up to a maximum of γ ∼ (2 − 5) within the simulated scales of (0.3 − 3.3) pc (Fig. 8). The flux budget evolution for the different simulated models reported in Fig. 10 clearly shows how the dominant energy flux at injection mainly drives the acceleration. Both the internal and magnetic energy channels are dissipated in turn to accelerate the jet, and therefore increase the kinetic energy of the beam. Moreover, we note that, in very specific cases, such as during the expansion phase, the internal energy can be dissipated to increase the jet speed, even in magnetically dominated jets. Our results on the acceleration mechanisms are mainly relevant in two different ways. (i) They confirm the current magnetic acceleration paradigm (Komissarov et al. 2007; Komissarov 2012). Indeed, in the case of cold, magnetically dominated jets, the magnetic energy is converted into kinetic energy, thus accelerating the bulk flow. (ii) They also expand the current view regarding thermal acceleration in internally dominated jets. In contrast to what is proposed by Vlahakis & Königl (2003a,b, 2004), where the authors claim that thermal acceleration only plays a role on compact scales before being overtaken by the magnetic acceleration, our findings indicate that when jets are thermodynamically relativistic (Perucho et al. 2017), thermal acceleration remains relevant even at parsec scales. Furthermore, thermal acceleration enables us to reach higher Lorentz factors than previously thought, aligning with those observed in AGN jets (in this specific case, in NGC 315). Remarkably, given that thermal acceleration can operate even when the jet undergoes conical expansion, our results suggest that jets can continue to accelerate on large scales, beyond the collimation region.

• We infer an explicit relation between the jet speed profiles and the intrinsic jet properties. The different velocity structures show a dependence on magnetic field configuration, that is, they differ depending on whether the configuration is force-free or nonforce-free, and/or depending on which energy flux dominates. When a Mach disk –that is, a strong internal planar shock– is formed due to the fast expansion, the acceleration deviates towards the outer and expanding layers. A mild acceleration on the spine is visible only toward the end of the grid of our simulation domain. On the contrary, when jets develop milder, conical shocks, the acceleration is concentrated in the central spine and smoothly decreases toward the external medium. The Mach disk is prominently visible in the force-free and magnetically dominated scenario. However, a thicker shear layer may delay the propagation and bouncing of the waves for long enough to avoid the formation of the Mach disk and to alter the downstream evolution of the jet, leading to a spine-accelerated structure. This result implies that the initial thickness of the shear layer, which is possibly associated with accretion disk-launched winds, is a relevant element in our understanding of the evolution and propagation of relativistic jets. While the acceleration focused on the outer layers can be a signature of magnetically dominated jets with thin sheaths, this does not necessarily exclude hot jets. In the nonforce-free configuration, the Mach disks are completely absent, and a classical expansion or recollimation pattern with spine acceleration is observed.

• The opening angle profiles shown in Fig. 9 highlight how the force-free models expand faster than their nonforce-free counterparts. This is attributed to the disengagement of the magnetic field force in the former scenario, which eliminates the contribution of the magnetic tension to the jet collimation. From the comparison of the opening angle (Fig. 9) and Lorentz factor (Fig. 8) profiles, we can distinguish two clear evolution cases: (i) the models that develop strong shocks and have the acceleration focused in the outer layers are those that show the fastest growth and highest maxima in the opening angle and (ii) those developing conical shocks with a downstream fast-spine structure are the ones with a shallower slope and lower maxima.

• As seen in Figs. 8 and 9, a number of models can reproduce both the Lorentz factor and the opening-angle profiles inferred for NGC 315 employing centimeter and millimeter VLBI observations. From the comparison with the observed jet speed in NGC 315, we suggest certain models are favored over others in their ability to model such a source. On the one hand, jets with power larger than 10⁴⁴ erg s⁻¹ show Lorentz factors at 3 pc that are remarkably larger than the observed one, giving an indication of a possible upper limit on the jet power of NGC 315. On the other hand, we highlight how different models with jet powers of ∼10⁴³ erg s⁻¹ match the different acceleration profiles obtained across the multiepoch observations at 5 GHz. This may suggest that the small changes in the jet properties at the injection may lead to consistent variations in the observed terminal Lorentz factor. Such variations can occur over timescales of a few months, indicating that the variations in the injection properties may drive the different observed speed profiles.

In conclusion, by comparing the results from 2D RMHD simulations with observational data, we suggest that thermal acceleration can play a relevant role in accelerating relativistic jets and should not be neglected when describing such a phenomenon. In the future, the inclusion of resistivity and a 3D grid will help to further understand the role of thermal acceleration in the context of the propagation of relativistic jets.

Acknowledgments

We would like to thank the referee for the comments, which highly improved the quality of the paper. We thank Eduardo Ros for his comments which improved the readability of the manuscript. Moreover, we thank Giancarlo Mattia for the useful discussion on the possible effects of magnetic reconnection. All the simulations were performed on the Cobra cluster of the Max Planck Society. L.R. and B.B. acknowledge the financial support of a Otto Hahn research group from the Max Planck Society. L.R. is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – project number 443220636. M.P. and J.M.M. acknowledge support by the Spanish Ministry of Science through Grants PID2019-105510GB-C31/AEI/10.13039/501100011033, PID2019-107427GB-C33, and from the Generalitat Valenciana through grants PROMETEU/2019/071 and ASFAE/2022/005. J.L.M. acknowledges support from the Generalitat Valenciana through grant ASFAE/2022/005.

References

Appendix A: Radial initial fluxes distribution

The convolution of the jet with the shear layer alters the initial analytical distribution of the fluxes, leading to the differences between the modeled and simulated values reported in Table 2.

In Fig. A.1, we show examples of the two main possibilities presented in this paper using models FMH1 and FMH1c. On the one hand, in the lower power models (left panel), the shear layer carries a relevant fraction of the total jet power, in the form of kinetic energy (in Fig. A.1, the kinetic energy absorbs the contribution of the rest mass). This is the consequence of the density in the jet increasing in the shear layer to match the external one. On the other hand, in the high-power, highly magnetized, models (right panel), the shear layer is magnetically dominated for the majority of its cross-section. As shown in Sects. 3 and 4, these differences do not lead to differences in the steady solutions (models FMH1 and FMH1c both develop the Mach disk and evolve with a similar global structure).

Fig. A.1. Relative initial radial distribution of the different energy channels for models FMH1 (left panel) and FMH1c (right panel). We highlight how the rest-mass contribution is absorbed into the kinetic energy. All the values are normalized to the maximum flux in each respective plot.

Appendix B: From initial conditions to equilibrium: pressure wave

In Fig. B.1 we highlight the evolution of the pressure wave at the three different time steps of t = 20 r_j/c, t = 40 r_j/c, and t = 70 r_j/c for model FMH1_m4. As explained in Sect. 3.1, since the initial jet is overpressured with respect to the external medium, once the simulation starts the jet abruptly expands towards it, giving the material radial velocities and pushing a fraction of it outside of the right boundary. As seen in Fig. B.2, which shows the axial pressure evolution of the ambient for different time steps, the passage of the pressure wave alters the pressure profile which oscillates around the original one. However, after a sufficient number of time steps, the lower and upper boundary conditions allow us to recover the initial, intended pressure ambient profile. Moreover, to test whether this would lead to differences in our final, steady solutions we run model FMH1_m4 with a larger grid of [n_x, n_z] = [100, 30], in order to avoid the ambient material to leave the grid. Fig. B.3 shows the Lorentz factor for model FMH1_m4 (same as Fig. 5 but cut at 1.2 pc) and the counterpart with the larger grid on the x-axis. No major differences in the final profiles are seen, proving how the lost material dragged out by the pressure wave does not affect our results and thus our conclusions.

Fig. B.1. Pressure maps at different times for the model FMH1_m4. Starting from left to right, the time steps are t = 20 rj/c, t = 40 rj/c, and t = 70 rj/c. The pressure wave moving from the left to the right boundary is clearly visible.

Fig. B.2. Ambient pressure evolution along the axial direction for the different time steps of t = 0 rj/c (red line), t = 70 rj/c (green line), t = 280 rj/c (blue line), t = 520 rj/c (orange line), and t = 1090 rj/c (black line). In the intermediate time steps, the oscillations in the pressure values are due to the passage of the wave. After enough time, such oscillations are absorbed by the boundary conditions and the original profile is recovered.

Fig. B.3. Lorentz factor maps for FMH1_m4 with the original grid (left plot) and the alternative grid of [nx, nz] = [100, 30]. This comparison shows how the material lost in the right boundary does not affect the steady solutions.

Appendix C: Resolution study

For our simulations, we chose a resolution of 20 cells/r_j. A first argument in support of this choice is given in Sect. 3.1 and in particular from the radial profiles shown in Fig. 4, which show how the radial evolution of the different physical quantities do not differ between 20 and 40 points. Further evidence in support of our chosen resolution is given in Fig. C.1 which shows the evolution of both the Lorentz factor (left panel) and opening angles (right panel) for two different models (FMH1 and NIH4) at the two discussed resolutions of 20 and 40 points. In order to save computational time, the simulations with 40 cells/r_j have a smaller grid of [n_x, n_z] = [30, 30], implying a physical distance of 0.3-1.2 pc. In the simulated distances, the two different resolutions lead to the same profiles in both tested models, confirming the validity of our choice of using the resolution of 20 cells/r_j.

Fig. C.1. Comparison of the profiles obtained from simulations at the resolution of 20 and 40 cells/rj. Left panel: Lorentz factor evolution up to 1.2 pc for models FMH1 and NIH4. Right panel: Half-opening angle obtained in the same models. Blue lines are for model NIH4 (continuous - 40 cells/rj; dotted - 20 cells/rj) while green lines are for model FMH1 (same line-style scheme). The comparison shows how the resolution of 20 cells/rj is enough to reach the saturation point, allowing us to use this resolution value for our simulations.

Appendix D: Further thicker shear layer models from Table 1

In this Appendix, we report in Table D.1 the other simulated models that we did not include in the main text for the sake of the readability of the paper. In detail, we report here the models for which we have tested the thicker shear layers with m = 2 and 4. As shown in Fig. D.1, in all the cases here reported the increased width of the shear layer with m = 4 is not sufficient to avoid the formation of the Mach disks (see Sect. 3). However, we highlight models FMH1c_m4 and FMH1c_m2. In the former, the shear layer with m = 4 is not enough to lead to the morphology change we discussed in the main text. Instead, the shear layer with m = 2 is thick enough to trigger the change in the acceleration pattern, from the outer layer to the internal spine. This model supports what we proposed in Sect. 4.3 concerning the role of the shear layer.

Fig. D.1. Lorentz factor maps for the other simulated models. The white contours represent the tracer at levels 0.2, 0.4, 0.6, and 0.8.

All Tables

All Figures

Fig. 1. Initial profiles along the radial distance in the different magnetic field configurations. We use the representative values of Bϕ = 0.24, Rm = 0.37, vj = 0.8, pa = 0.01, and Bz = 0.05 (for the nonforce-free configuration) in code units. Upper panels: nonforce-free configuration. Lower panel: force-free configurations. Left panels: Magnetic field distribution of the azimuthal (blue line) and axial (green line) components. Right panels: profile distribution of the gas pressure (blue line), magnetic pressure (green line), and total pressure (black line). The dotted black line indicates the shear layer with m = 8 normalized at the higher value in the plot for an easier representation.

In the text

Fig. 2. Expected magnetic field distribution between 0.1 pc and 3 pc. The upper and lower magnetic field limits are extrapolated from the core shift measurement under the different assumption (for the details see Ricci et al. 2022). The black vertical line highlights the starting point of zi = 0.3 pc of our simulations, while the horizontal dashed lines are the expected range of magnetic field at zi.

In the text

Fig. 3. Absolute value of radial velocity for the model FMH1_m4. The white contours highlight the tracer values at levels of 0.2, 0.4, 0.6, and 0.8 (the tracer levels are in decreasing order, i.e., from 0.8 to 0.2, from left to right). Outside of the jet structure delineated by the shear layer (white contours), the velocity is smaller than 0.01 everywhere in the grid, meaning that this can be considered a steady solution.

In the text

Fig. 4. Radial profiles of the different physical properties for model FMH1 at different distances from the injection (cell 1 in red, half grid in blue, and final cell in green) and at different times (initial in dotted lines and final in continuous lines). In the top panels, we show the radial profiles with a resolution of 10 points, in the middle the ones at a resolution of 20 points, and in the bottom ones the profiles with 40 points (here, the half grid is not shown due to the smaller simulated grid). The convergence in the profiles between 20 and 40 points led us to use resolutions of 20 points for our simulations.

In the text

Fig. 5. Lorentz factor maps for the simulated models. The white contours are representative of the tracer at levels of 0.2, 0.4, 0.6, 0.8. The results show the clear correlation between the the jet speed profiles and the intrinsic properties, with the main differences shown between force-free and nonforce-free along with magnetically dominated or internally dominated models. For details, see the text.

In the text

	Fig. 6. Pressure maps for the simulated models. The white contours are representative of the tracer at levels of 0.2, 0.4, 0.6, 0.8. As in Fig. 5 the different jet nature leads to the different pressure profiles in the spine and in the shear layer regions.
In the text

Fig. 7. Modulus of the gradient of the rest-mass density (in logarithmic scale) for models FMH1 (left panels) and FMH1_m4 (right panels). The upper panels show the entire grid, while in the bottom ones we show a zoom-in in the 0.3–0.9 pc region. The modulus of the gradient of the rest-mass density highlights the location of steep sound waves and shock waves. The black squares highlight the triple point associated with the incident/reflected/Mach shocks in the jet-ambient transition surface (for details, see Sect. 4.3). In the bottom panels, the contours represent the tracer at the level f = 0.9, showing the transition between the inner jet and the shear layer.

In the text

Fig. 8. Radially averaged Lorentz factor as a function of distance for all the simulated models. Top panel: magnetically dominated jet models. Bottom panel: thermally and kinetically dominated jet models. The different data points are the Lorentz factor for different epochs and frequencies of NGC 315 inferred by Ricci et al. (2022). The continuous lines represent the force-free magnetically dominated models, the dashed lines the force-free thermally dominated ones, the dotted lines are for the kinetically dominated model, and the dashed-dotted lines for the nonforce-free configurations.

In the text

Fig. 9. Half-opening angle as a function of distance for all the simulated models. Top panel: magnetically dominated jet models. Bottom panel: thermally and kinetically dominated jet models. The force-free magnetically dominated setups are shown in continuous lines, the force-free thermally dominated in dashed lines, the kinetically dominated in dotted lines, and the nonforce-free in dashed-dotted lines. The data points represent the intrinsic half-opening angle profiles of NGC 315 obtained by Ricci et al. (2022) for multi-frequency and multi-epoch VLBI observations.

In the text

Fig. 10. Evolution of the integrated energy fluxes along the axial direction for all the models reported in Table 1. Black, red, green, blue, and orange lines represent, respectively, the total, kinetic, internal, magnetic, and rest-mass energy fluxes. Depending on the dominating energy channel, its dissipation in favor of the jet kinetic flux is visible in all the models. An exception is the kinetically dominated model FKS1, in which both the thermal and magnetical channels are dissipated to accelerate the jet, even if not dominant.

In the text

	Fig. A.1. Relative initial radial distribution of the different energy channels for models FMH1 (left panel) and FMH1c (right panel). We highlight how the rest-mass contribution is absorbed into the kinetic energy. All the values are normalized to the maximum flux in each respective plot.
In the text

	Fig. B.1. Pressure maps at different times for the model FMH1_m4. Starting from left to right, the time steps are t = 20 rj/c, t = 40 rj/c, and t = 70 rj/c. The pressure wave moving from the left to the right boundary is clearly visible.
In the text

Fig. B.2. Ambient pressure evolution along the axial direction for the different time steps of t = 0 rj/c (red line), t = 70 rj/c (green line), t = 280 rj/c (blue line), t = 520 rj/c (orange line), and t = 1090 rj/c (black line). In the intermediate time steps, the oscillations in the pressure values are due to the passage of the wave. After enough time, such oscillations are absorbed by the boundary conditions and the original profile is recovered.

In the text

	Fig. B.3. Lorentz factor maps for FMH1_m4 with the original grid (left plot) and the alternative grid of [nx, nz] = [100, 30]. This comparison shows how the material lost in the right boundary does not affect the steady solutions.
In the text

Fig. C.1. Comparison of the profiles obtained from simulations at the resolution of 20 and 40 cells/rj. Left panel: Lorentz factor evolution up to 1.2 pc for models FMH1 and NIH4. Right panel: Half-opening angle obtained in the same models. Blue lines are for model NIH4 (continuous - 40 cells/rj; dotted - 20 cells/rj) while green lines are for model FMH1 (same line-style scheme). The comparison shows how the resolution of 20 cells/rj is enough to reach the saturation point, allowing us to use this resolution value for our simulations.

In the text

	Fig. D.1. Lorentz factor maps for the other simulated models. The white contours represent the tracer at levels 0.2, 0.4, 0.6, and 0.8.
In the text