© ESO, 2016

1. Introduction

Extragalactic radio sources have been classified into two categories (Fanaroff & Riley 1974) based upon their radio morphology: a first class of objects, named Fanaroff-Riley I (FR I), which is preferentially found in rich clusters and hosted by weak-lined galaxies, shows jet-dominated emission and two-sided jets at the kiloparsec scale that smoothly extend into the intracluster medium, where they form large-scale plumes or tails of diffuse radio emission. The second class, named Fanaroff-Riley II (FR II, or classical doubles), found in poorer environments and hosted by strong emission-line galaxies, presents lobe-dominated emission and one-sided jets at the kpc scale that abruptly terminate into hot-spots of emission.

In addition to morphology, FR I and FR II radio sources have been distinguished based on power: objects below ~${\mathrm{10}}^{\mathrm{25}}{\mathit{h}}_{\mathrm{70}}^{\mathrm{2}}$ W Hz^-1 str^-1 at 178 MHz were typically found to be FR I sources. A perhaps more illuminating criterion was found by Ledlow & Owen (1996), who plotted the radio luminosity at 1.4 GHz against the optical absolute magnitude of the host galaxy: they found the bordering line of FR I to FR II regions to correlate as ${\mathit{L}}_{\mathrm{R}}\mathrm{\propto }{\mathit{L}}_{\mathrm{opt}}^{\mathrm{1.7}}$, that is, in a luminous galaxy more radio power is required to form a FR II radio source. This correlation is important since it can be interpreted as an indication that the environment may play a crucial role in determining the source structure. Furthermore, a class of hybrid sources has been discovered that shows FR I structure on one side of the radio source and FR II morphology on the other (Gopal-Krishna & Wiita 2000). These arguments describe the basic question of the origin of the FR I and FR II dichotomy, whether the shapes are intrinsic or moulded by the environment (Gopal-Krishna & Wiita 2000; Wold et al. 2007; Kawakatu et al. 2009; Massaglia 2003).

Recent studies based on the cross correlation of wide-area optical and radio surveys (Best & Heckman 2012) unveiled yet another class of compact radio-galaxies representing the majority of the local radio-loud AGN population. Because they lack the prominent extended radio structures characteristic of the other FR classes, they were dubbed FR 0 (Baldi & Capetti 2009, 2010; Sadler et al. 2014; Baldi et al. 2015). Nonetheless, FR 0s often show radio jets, but they extend only a few kpc at most (Baldi et al. 2015). This suggests that in most radio-loud AGN, the jets fail to propagate to (or become too faint to be detected at) radii exceeding the size of their host.

The distorted, diffuse, and plume-like morphologies of FR I sources led researchers to model them as turbulent flows (Bicknell 1984, 1986; Komissarov 1990a,b; De Young 1993), while the characteristics of FR II, such as their linear structure and the hot-spots at the jet termination, are associated with hypersonic flows. The difference in morphology between the two classes therefore reflects a difference in how the jet energy is dissipated during jet propagation: in the first case, the jet gradually dissipates its energy and is characterized by entrainment of the ambient material, while in the second case, the jet retains its velocity and dumps all its energy at its termination, forming the observed hot-spots. While the dynamics of jets with high Mach number has been widely studied by means of numerical simulations in 2D and 3D by, for example, Massaglia et al. (1996), Zanni et al. (2003), Hardcastle & Krause (2013), Hardcastle & Krause (2014), 3D simulations of transonic jets have been carried out following the evolution of instabilities to turbulence (Basset & Woodward 1995; Hardee et al. 1995; Loken 1997; Bodo et al. 1998), and simulations of turbulent jets that include the jet head propagation arelimited to Nawaz et al. (2014, 2016), who investigated the properties of the jet in Hydra A. With this paper, we therefore intend to start a systematic study of these flows. We perform three-dimensional simulations of the propagation of jets with low Mach number in a stratified medium, which is meant to model the interstellar-intracluster transition. In this first paper we perform hydrodynamic simulations and neglect the effects of the magnetic field.

An important point to consider is that both FR I and FR II radio sources show evidence of relativistic flows at the parsec scale, and therefore a deceleration to sub-relativistic velocities must occur in FR Is between the inner region and the kiloparsec scale (Giovannini et al. 1994, 2001; Laing & Bridle 2014). Several models have been proposed for the deceleration mechanism (Bicknell 1994, 1995; Komissarov 1994; Bowman et al. 1996; De Young 2005) and numerical simulations of these processes have been performed (Perucho & Martí 2007; Rossi et al. 2008; Tchekhovskoy & Bromberg 2016). In this paper, however, we assume that the deceleration has already taken place, and we consider a scale where the jet is non-relativistic. Moreover, we stress that to study the transition to turbulence and the turbulent structure of these flows, it is essential to perform the simulations in three dimensions. The behavior of 2D jets with the same physical parameters is completely different from their 3D counterparts.

The plan of the paper is the following: in Sect. 2 we describe the numerical setup and show the equations we solve, Sect. 3 outlines the observational framework and constrains the physical parameters, and in Sect. 4 we present and discuss the obtained results. The conclusions are drawn in Sect. 5.

2. Numerical setup

2.1. Equations and method of solution

For our purposes, we solve the Euler equations of gas dynamics, which can be written in quasi-linear form as

The quantities ρ, P, and v are the density, pressure, and velocity, respectively, and Γ = 5/3 is the ratio of the specific heats. The jet and the external material are distinguished using a passive tracer, f, set equal to unity for the injected jet material and equal to zero for the ambient medium. We note that the choice of Γ = 5/3 is consistent with a supersonic nonrelativistic jets with Mach numbers >4 (see discussion below).

Equations (1)–(4) were solved using the piecewise parabolic method (PPM) and the HLLC Riemann solver of the PLUTO code (Mignone et al. 2007, 2012).

2.2. Initial and boundary conditions

The 3D simulations were carried out on a Cartesian domain with coordinates in the range x ∈ [−L/ 2,L/ 2], y ∈ [0,L_y] and z ∈ [−L/ 2,L/ 2] (lengths are expressed in units of the jet radius r_j and y is the direction of jet propagation). At t = 0, the domain is filled with a perfect gas at uniform pressure but spherically stratified according to a King-like profile in density: $\mathit{\rho }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{\eta }}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}{\mathrm{(}\mathit{r}\mathit{/}{\mathit{r}}_{\mathrm{c}}\mathrm{)}}^{\mathit{\alpha }}}\mathit{,}$(5)where the density is expressed in units of the jet density ρ_j, $\mathit{r}\mathrm{=}\sqrt{{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}}$ is the spherical radius, r_c is the core radius, and η is the ratio ρ_j/ρ_c between the jet density and the core density. We set r_c = 40r_j throughout and considered different values of the parameters α and η. The ambient temperature increases with radius as $\mathit{T}\mathrm{\propto }\mathrm{1}\mathrm{+}\mathrm{(}\mathit{r}\mathit{/}{\mathit{r}}_{\mathrm{c}}\mathrm{)}\alpha$(6)to maintain the pressure uniform.

We imposed zero-gradient boundary conditions on all computational boundaries with the exception of the injection region located at y = 0. Here we prescribed inside the unit circle a constant cylindrical inflow directed along the y direction, and in all but one case (see Table 1), in pressure equilibrium with the ambient. The inflow jet values are given by $\begin{array}{c}\\ {\mathit{\rho }}_{\mathrm{j}}\mathrm{=}\mathrm{1}\mathit{,}\\ {\mathit{v}}_{\mathit{y,}\mathrm{j}}\mathrm{=}\mathit{M,}\\ {\mathit{p}}_{\mathrm{j}}\mathrm{=}\mathrm{1}\mathit{/}\mathrm{\Gamma }\mathit{,}\\ {\mathit{f}}_{\mathrm{j}}\mathrm{=}\mathrm{1}\mathit{,}\end{array}$(7)where velocity is measured in units of the jet sound speed on the axis, c_sj, therefore M represents the jet Mach number, and pressure is measured in units of ${\mathit{\rho }}_{\mathrm{j}}{\mathit{c}}_{\mathrm{sj}}^{\mathrm{2}}$. Outside the jet nozzle, reflective boundary conditions hold. To avoid sharp transitions, we smoothly joined the injection and reflective boundary values in the following way: $\mathit{Q}\mathrm{\left(}\mathit{x,z,t}\mathrm{\right)}\mathrm{=}{\mathit{Q}}_{\mathit{r}}\mathrm{\left(}\mathit{x,z,t}\mathrm{\right)}\mathrm{+}\frac{{\mathit{Q}}_{\mathrm{j}}\mathrm{-}{\mathit{Q}}_{\mathit{r}}\mathrm{\left(}\mathit{x,z,t}\mathrm{\right)}}{\cosh \left[\mathrm{(}\mathit{R}\mathit{/}{\mathit{R}}_{\mathrm{s}}{\mathrm{)}}^{\mathit{n}}\right]}\mathit{,}$(8)where Q = { ρ,v_x,ρv_y,v_z,p,f } are primitive flow quantities with the exception of the jet velocity, which was replaced by the y-momentum, while $\mathit{R}\mathrm{=}\sqrt{{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}}$ is the cylindrical radius. We note that Q_r(x,z,t) are the corresponding time-dependent reflected values, while Q_j are the constant injection values given by Eqs. (7). In Eq. (8) we set R_s = 1 and n = 6 for all variables except for the density, for which we used R_s = 1.4 and n = 8. This choice ensured monotonicity in the first (ρv_y) and second $\mathrm{\left(}\mathit{\rho }{\mathit{v}}_{\mathit{y}}^{\mathrm{2}}\mathrm{\right)}$ fluid outflow momenta (Massaglia et al. 1996). We did not explicitly perturb the jet at its inlet, differently from Mignone et al. (2010). The growing non-axially symmetric modes originate from numerical noise.

The physical domain was covered by N_x × N_y × N_z computational zones that are not necessarily uniformly spaced. For domains with a large physical size, we employed a uniform grid resolution in the central zones around the beam (typically for | x | , | z | ≤ 6) and a geometrically stretched grid elsewhere.In the central zones we have a resolution of 10 grid points per jet radius for the cases of M = 4,10, and 13 for the case of M = 40. As shown in Anjiri et al. (2014), a resolution of 12 points per jet radius produces results that agree well, within a few percent, with computations performed by doubling this resolution. The comparison between different resolutions was carried out considering both the morphological and the energy budget aspects. The complete set of parameters for our simulations is given in Table 1.

3. Observational background and parameter constraints

To obtain physical values from the results of our simulations, we set constraints derived from observational data to the values of the fundamental units of length r_j, density ρ_j, and velocity c_sj.

As discussed in the introduction, FR I jets propagate at relativistic velocities close to the central AGNs at the parsec scales, while they are non-relativistic or weakly relativistic at the kiloparsec scales (Giovannini et al. 1994, 2001). This means that a jet-braking process must have taken place between the inner and outer scales (Rossi et al. 2008; Laing & Bridle 2014). Our goal here is not to study this deceleration process, however, but to consider the jet propagation at scales starting from a few hundred pc, where the jet is already non-relativistic, to the typical linear size of FR I radio sources, which range from about 10 up to about hundred of kiloparsec (Parma et al. 1987, for the B2 catalog sources).

At these spatial scales, the jet propagates in the galactic interstellar medium for the first few kiloparsec before exiting into the intracluster medium. The jet radius in most FR I sources is at the limits of the angular resolution of radio interferometers such as the VLA, ~0.1−0.3′′ (Parma et al. 1987). We assumed a fiducial value for the jet radius of r_j = 100 pc. A consistent limit on the jet radius is also determined by requirements on numerical resolution: since we require the grid physical extension to be about ten kiloparsec, the jet radius cannot be much smaller than 100 pc with the given grid size if we aim to cover it with a number of points not much smaller than 10. In the following, the values of the spatial coordinates are given in kiloparsec within this assumption for the initial jet radius.

Fig. 1 3D isocontours of the tracer distribution for case A (M = 4, η = 10-2, ℒkin = 1.1 × 1042 erg s-1) at t = 740 time units, 5.7 × 107 yr. The size of the computational box is 6.4 × 12 × 6.4 kpc.

The typical galactic core radius r_c attains values of about a few kiloparsec, and we assumed r_c = 4 kpc, which means 40 jet radii. X-ray observations have shown (Posacki et al. 2013) that within the galactic core, the average temperature of the ambient medium is T_c ~ 0.1−0.3 keV, with a particle density of n_c ~ 1 cm^-3 (Balmaverde et al. 2008). As the jet enters the intragroup/intracluster medium, the ambient gas is hotter by about one order of magnitude (Sun et al. 2009; Connor et al. 2014) and the density drops by approximately the same amount. This is consistent with assuming the parameter α in Eq. (5) to be between 1 and 2.

The jet parameters can be expressed as functions of the external density, ρ_c, the external temperature, T_c, the Mach number M, the density ratio η, and the ratio between jet and ambient pressure as follows:

$\begin{array}{ccc}& & {\mathit{\rho }}_{\mathrm{j}}\mathrm{=}\mathrm{0.01}\left(\frac{\mathit{\eta }}{\mathrm{0.01}}\right)\hspace{0.17em}\left(\frac{{\mathit{\rho }}_{\mathrm{c}}}{\mathrm{1}{\mathrm{cm}}^{-3}}\right)\hspace{0.17em}{\mathrm{cm}}^{-3}\mathit{,}\\ & & {\mathit{v}}_{\mathrm{j}}\mathrm{=}\mathrm{5.1}\mathrm{\times }{\mathrm{10}}^{\mathrm{8}}\mathit{M}\hspace{0.17em}{\left(\frac{{\mathit{T}}_{\mathrm{c}}}{\mathrm{0.2}\mathrm{keV}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}{\left(\frac{\mathit{\eta }}{\mathrm{0.01}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}{\left(\frac{{\mathit{P}}_{\mathrm{j}}}{{\mathit{P}}_{\mathrm{c}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{cm}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{.}\end{array}$

We can also derive the jet kinetic power as $\begin{array}{ccc}{\mathrm{ℒ}}_{\mathrm{kin}}& \mathrm{=}& \mathrm{1.1}\mathrm{\times }{\mathrm{10}}^{\mathrm{42}}{\left(\frac{{\mathit{r}}_{\mathrm{j}}}{\mathrm{100}\mathrm{pc}}\right)}^{\mathrm{2}}\hspace{0.17em}\left(\frac{{\mathit{\rho }}_{\mathrm{c}}}{\mathrm{1}{\mathrm{cm}}^{-3}}\right)\hspace{0.17em}{\left(\frac{{\mathit{T}}_{\mathrm{c}}}{\mathrm{0.2}\mathrm{keV}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}\\ & & \mathrm{\times }{\left(\frac{\mathit{M}}{\mathrm{4}}\right)}^{\mathrm{3}}\hspace{0.17em}{\left(\frac{\mathit{\eta }}{\mathrm{0.01}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}{\left(\frac{{\mathit{P}}_{\mathrm{j}}}{{\mathit{P}}_{\mathrm{c}}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{erg}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{.}\end{array}$(11)The computational time unit τ is the sound travel time over the initial jet radius, which in physical units corresponds to $\begin{array}{ccc}\mathit{\tau }& \mathrm{=}& \frac{{\mathit{r}}_{\mathrm{j}}}{{\mathit{c}}_{\mathrm{sj}}}\mathrm{=}\mathrm{7.7}\mathrm{\times }{\mathrm{10}}^{\mathrm{4}}\left(\frac{{\mathit{r}}_{\mathrm{j}}}{\mathrm{100}\mathrm{pc}}\right)\hspace{0.17em}{\left(\frac{{\mathit{T}}_{\mathrm{c}}}{\mathrm{0.2}\mathrm{keV}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\\ & & \mathrm{\times }{\left(\frac{\mathit{\eta }}{\mathrm{0.01}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}{\left(\frac{{\mathit{P}}_{\mathrm{j}}}{{\mathit{P}}_{\mathrm{c}}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{yr}\mathit{.}\end{array}$(12)Our typical simulation covers about 100 to 1000 time units, corresponding to about 10⁷−10⁸ yr of the source lifetime, while it has traveled more than 10 kpc in the ambient medium.

Various methods for converting the kinetic jet power ℒ_kin into observed radio luminosity L_r have been proposed; we are here mainly interested in an estimate of ℒ_kin corresponding to the power at which the FR I/FR II transition occurs, ℒ_trans. Willott et al. (1999) obtained a relation according to which ℒ_kin scales as ${\mathit{L}}_{\mathrm{r}}^{\mathrm{6}\mathit{/}\mathrm{7}}$, leading to ℒ_trans ~ 1.4 × 10⁴³ erg s^-1. Birzan et al. (2004) estimated ℒ_kin from the observations of the X-ray cavities inflated by radio AGN. The value of ℒ_trans based on their calibration support the above estimate, returning values in the range ℒ_trans ~ 3 × 10⁴²−5 × 10⁴³ erg s^-1 at the FR I/FR II power transition.

We performed various simulations that we describe in the following section, starting at ℒ_kin = 1.1 × 10⁴² erg s^-1, within the expected range of FR I power, and then increased the jet power into the FR II regime.

Fig. 2 Results for case A (M = 4, η = 10-2, ℒkin = 1.1 × 1042 erg s-1) at t = 740 time units, 5.7 × 107 yr. Left: cut in the (y,z) plane of the logarithmic density distribution in units of jet density. The spatial units are in jet radii (100 pc) and the image extends over 12 kpc in the jet direction and over 6 kpc in the transverse direction. Right: maximum pressure on a (x,z) plane at a given position y along the jet as a function of y. The pressure is plotted, from now on, in units of the ambient value (instead of the computational units ${\mathit{\rho }}_{\mathrm{j}}{\mathit{c}}_{\mathrm{sj}}^{\mathrm{2}}$).

Fig. 3 Results for a 2D simulation with the same jet parameters as case A and at t = 500 time units, 4 × 107 yr.

4. Results and discussion

We have carried out a series of simulation runs to explore the effects of the different parameters on the jet propagation, evolution, and resulting morphologies.

Although we simulated jets propagating in a stratified medium, which reproduces the real ambient medium of radio galaxies, we started our analysis with an initial simulation (case A in Table 1) in which the ambient density was constant. This was motivated by the need of comparing our results with those obtained in previous works. The jet has a Mach number M = 4 and a density contrast of η = 10^-2 corresponding to ℒ_kin = 1.1 × 10⁴² erg s^-1, well below the FR I/FR II power transition. The outcome of the simulation is presented in the 3D image shown in Fig. 1. After an initial collimated phase, the jet disrupts and expands into the ambient medium. We present an image of the jet tracer and not of emission, but it can be envisaged that this diffuse gas distribution will give rise to an FR I radio morphology.

This is confirmed by Fig. 2, in which we show the maximum pressure on a transverse (x,z) plane at a given position y along the jet as a function of y. The pressure rises along the jet, reaching its maximum at r ~ 40, and then it steadily decreases. At the jet head there is only a small increase, ~20%, indicating a weak shock.

We simulated a jet described by the same parameters, but in 2Dcylindrical geometry (axisymmetric), see Fig. 3, and the results are radically different. The jet remains collimated over its whole length. The pressure displays several local maxima associated with the jet internal shocks and, more importantly, with its head.

This can be explained by recalling that the properties of turbulence in two and three dimensions are very different because in two dimensions we do not have the energy cascade to small scales. Consequently, the entrainment properties and the jet evolution become strongly diversified. In particular, the FR I morphology resulting from the propagation of a low-power jet arises and can be explored only with high-resolution 3D simulations. To extract physical information from our simulations in a more realistic framework, all the cases that follow were carried out in the stratified medium described in Sect. 2.

In the top panel of Fig. 4 we show a cut in the (y,z) plane of the logarithmic density distribution for case B of Table 1 (M = 4, η = 10^-2, the same parameters as for case A, but including the external density stratification) at 5.2 × 10⁷ yr of the jet lifetime and at a traveled distance of 12 kpc. The general features are very similar to those shown by case A, without ambient gas stratification. No sign of the head Mach disk is observed, and the jet smoothly mixes with the ambient medium. This is confirmed by the cuts of the tracer distribution and of the longitudinal velocity of the matter belonging to the external medium,that is, the quantity (1−f) × v_y in the central and bottom panels of Fig. 4.We also note that the bow shock has disappeared from the ambient medium since the head propagation has become subsonic at these distances. The jet remains well collimated out to 30–40 jet radii, then it rapidly widens at large distances. This is the same value of the core radius of the external gas distribution; however, this is observed also in case A, where no stratification is present, and it should be considered as a coincidence. This distance is most likely related to the growth length of instabilities that cause the transition to turbulence. At the same location it decelerates and entrains external material that is also globally accelerated. This behavior is preserved qualitatively out to the jet head, with just a further widening and deceleration.

Fig. 4 Results for our reference case B (M = 4, η = 10-2, ℒkin ~ 1.1 × 1042 erg s-1), at t = 640 time units, 4.9 × 107 yr). The spatial units are in jet radii (100 pc) and the image extends over 12 kpc in the jet direction and over 6 kpc in the transverse direction. Cuts in the (y,z) plane of (top) the logarithmic density distribution in units of jet density, (center) tracer distribution, and (bottom) distribution of the longitudinal velocity of the external medium in units of csj.

As first pointed out by Bodo et al. (1994), who studied the nonlinear temporal evolution of jet velocity shear instabilities (Kelvin-Helmholtz instabilities) in 2Dcylindrical symmetry, and then by Bodo et al. (1998) in the 3D extension of this investigation, the evolution of unstable modes follows a trend that is essentially unchanged in the general features and that develops in three phases: (1) a first linear phase where the unstable modes grow according to the linear theory until internal shocks start to form; (2) this is followed by an acoustic phase where the growth of the internal shocks is accompanied by a global deformation of the jet, which drives shocks into the external medium; these shocks carry both momentum and energy away from the jet and transfer them to the external medium; (3) eventually a final mixing phase where, as a consequence of the shock evolution, turbulent mixing between jet and the ambient medium begins. In Fig. 5(top panel) the cut of the pressure distribution of the external medium, i.e. (1−f) × P, is suitable to show how these three phases develop in space.In Fig. 5 (bottom panel) we show a detail of the total pressure in the central and initial section of the domain that is relevant for the linear and, partly, the acoustic phase. Oblique shocks internal to the jet are observed.

Fig. 5 Same as in Fig. 4, but for the pressure distribution of the external medium. The three phases of the evolution are labeled (top panel).The bottom panel shows the total pressure distribution in the initial part of the jet with internal shocks.

The main properties obtained for case B are also reproduced in case C, which differs only for the parameter α, which describes the density profile of the ambient medium (see Fig. 6). We only note a slightly slower (by ~10%) advance of the jet head, which is due to the slower density decrease of the external gas.

In the cases considered so far, we assumed that the jet is initially in pressure equilibrium with the ambient medium. This choice is motivated by the fact that an overpressured (underexpanded) jet develops a recollimation shock that leads to equal pressure between the internal and external gas (Belan et al. 2010). This probably occurs at radii smaller than those probed by our simulations, within the initial 100 pc of the jet length, and the jet is therefore already in pressure equilibrium at injection. It is nonetheless useful to assess the main features of the evolution of an overpressured jet. We simulated this with P_jet/P_amb = 10, as case D. The general morphology is similar to that of case A (see Fig. 7), but, as expected, the jet shows a recollimation shock at r < 10 where it has already expanded to twice its initial radius. This rapid spread causes the jet disruption to occur at smaller distances, r ~ 20, than in case A. Nonetheless, the jet advances at a very similar speed.

Fig. 6 Cut in the (y,z) plane of the logarithmic density distribution for case C, differing from the reference case B for the value of α that describes the external medium, at t = 740 time units, 5.7 × 107 yr.

Fig. 7 Cut in the (y,z) plane of the logarithmic density distribution for case D, overpressured jet, at t = 720 time units, 1.7 × 107 yr.

Fig. 8 Cut in the (y,z) plane of the logarithmic density distribution for case E, the highest jet power (M = 40, η = 10-2, ℒkin = 1.1 × 1045), at t = 36 time units, 3 × 106 yr.

Fig. 9 Maximum pressure on transverse planes vs y for case B at t = 640 time units (left panel), case E at t = 36 (right panel).

This investigation is characterized by two main parameters, which are the density ratio η and the Mach number M, as defined in Sect. 2. The cases presented in detail are characterized by the lowest value of the kinetic power. Exploring the parameter plane, we reach higher power and note a transition between FR I-like to FR II-like morphologies. This becomes evident by comparing the density distributions of Figs. 4 and 8 relative to the cases B and E, respectively. Case E differs from case B only for its Mach number, 40 instead of 4, but this translates into a 1000 higher kinetic power. Case E is characterized by the presence of a shocked region at the jet head, with the consequent cocoon that shrouds and separates the shocked region by the external unperturbed medium; these are characteristics typical of FR II sources.

A more quantitative representation of the differences between these two simulations is obtained by examining the plots of the maximum pressure as functions of the longitudinal distance given in Fig. 9. In the low-power case B (left panel) the pressure reaches its maximum value at ~40 jet radii and then it steadily decreases with only a slight increase at the jet head. In the high-power case E (center panel) the maximum pressure shows a strong increase toward the jet termination point, in addition to small-scale peaks (which probably are the locations of the internal shocks).

At which Mach number does this transition occur, however? Maintaining η constant, we performed a simulation with a more moderate Mach number in case F, M = 10 (corresponding to ℒ_kin = 1.7 × 10⁴³). Figure 10 shows the cocoon and the shocked region at the jet head, as confirmed by the maximum pressure behavior in Fig. 9 (right panel), where shocks are still clearly visible. Thus, at η = 0.01 the transition between FR I and FR II morphologies occurs between M = 4 and M = 10.

Finally, we considered case G with the same Mach number (M = 4) as in the reference case, but a lower density (η = 10^-3). The jet kinetic power is ℒ_kin = 3.5 × 10⁴², ~3 times higher than case B. It is counterintuitive that a lighter jet has a higher jet power. However, from the condition of pressure equilibrium, a lighter jet is also hotter and therefore has a higher velocity for the same Mach number. The resulting density distribution and pressure profile (see Fig. 11) indicates that this case corresponds to a FR I morphology.

From our simulations we recover a separation between FR I and FR II morphologies that occurs at a jet power that agrees remarkably well with the observations.

Fig. 10 Left: cut in the (y,z) plane of the logarithmic density distribution for case F (M = 10, η = 10-2, ℒkin = 1.7 × 1043) at t = 98 time units, 7.5 × 106 yr; right: maximum pressure on transverse planes along the jet.

Fig. 11 Left: cut in the (y,z) plane of the logarithmic density distribution for case G (M = 4, η = 10-3, ℒkin = 3.5 × 1042) at t = 1100 time units, 8.5 × 107 yr; right: maximum pressure on transverse planes along the jet.

Fig. 12 Jet head position as a function of time (solid line) compared with the theoretical estimate for a uniform medium (dashed line) for cases B, E, and F in the left, center, and right panel, respectively.

Another important difference between the cases producing FR I and FR II is the advance speed of the jet. The progress of the jet into the ambient medium is shown in Fig. 12 for case B (left panel), where the solid line represents the position of the jet head as a function of time and the dashed line, given as a reference, is obtained by the longitudinal momentum conservation at the jet head in a uniform medium (Martí et al. 1994): ${\mathit{y}}_{\mathrm{head}}\mathrm{=}\frac{\mathit{M}}{\mathrm{1}\mathrm{+}\mathrm{1}\mathit{/}\sqrt{\mathit{\eta }}}\hspace{0.17em}\mathit{t}\mathit{.}$(13)We note that, notwithstanding the longitudinal decrease of the ambient medium, which favors an increase in jet propagation speed, the jet slows down because the Mach disk is disrupted. The jet advance shows a knee at t ~ 100 time units where the velocity drops by a factor ~2. After the knee the velocity approaches a constant value of v_head ~ 500 km s^-1.

Conversely, in cases E and F (Fig. 12, center and right panel, respectively) the head progress rate increases when it reaches r ~ 20 as a result of the drop in the external density. This is because, unlike what is seen in case B, the jet remains well collimated. A subtle difference is nonetheless present as in case F the advance speed appears to decrease at r ≳ 80, where indeed the jet starts to loose its coherence, while it remains constant for the more powerful jet simulated with case E.

5. Summary and conclusions

We have performed three-dimensional numerical simulations of turbulent jets, which generate morphologies typical of FR I radio sources. The jet propagates in a stratified medium that is meant to model the interstellar intracluster transition. FR I radiosources are known to be relativistic at the parsec scale, therefore a deceleration to sub-relativistic velocities must occur between this scale and the kiloparsec scale. In this paper we did not model the deceleration process, but we assumed that deceleration already occurred and considered scales where the jet is sub-relativistic.We did not consider buoyancy effects, by not taking into account the galaxy gravitational potential, but they are expected to be negligible up to the distances reached by the jet head in our simulations (of about ten kpc), while they may become more relevant at lager distances, as the head further slows down. In this first analysis we also neglected the effect of the magnetic field. The parameters governing the jet evolution are therefore the Mach number M and the initial jet-to-ambient density ratio η, which, by constraining the values of the external density and temperature through observational data, can be combined to give the jet kinetic power.

The estimated jet kinetic power of the transition between FR I and FR II is 10⁴³ erg s^-1, and we investigated a series of cases below this threshold. These low-power cases have M = 4 and different values of density ratio, and they all give rise to turbulent structures typical of FR I sources. The jet power, instead of being completely deposited at the termination through a series of terminal shocks, as in FR II sources, is gradually dissipated by the turbulence. We showed that three-dimensionality is an essential ingredient for the occurrence of the transition to turbulence. Two-dimensional simulations with the same parameters lead to FR II like behavior with energy dissipation concentrated at the jet termination. Increasing the Mach number to 40 and consequently the kinetic power well above the FR I – FR II threshold, we obviously recover well-collimated jets that dump all their energy at the termination shocks. At intermediate Mach numbers (M = 10), with a kinetic power around the transition value, we find characteristics typical of FR II sources, even though the energy deposition at the jet termination starts to become more gradual and the morphology acquires some of the FR I properties.

The simulations presented show that in FR Is the jet energy is transferred to the ISM in part inducing, through entrainment, a global low-velocity outflow; the remaining power is instead dissipated through acoustic waves. The energy released by active nuclei is thought to play a fundamental role in the evolution of their host galaxies (Fabian 2012); in particular, in radio-loud AGN, the kinetic energy carried by their jets is transferred to the surrounding medium, which leads to the so-called radio-mode feedback. The FR I jets, although of lower power with respect to those of FR IIs, are extremely important from the point of view of feedback. This is because the FR I jets remain confined within the central regions of the host over longer timescales (and possibly for their whole lifetime) which exceeds 10⁷ yr for our reference case B. Furthermore, the entire host is affected by the radio-mode feedback, while in more powerful radiosources a smaller volume (immediately surrounding the jets) is involved. Finally, as a result of the steep radio luminosity function (e.g. Mauch & Sadler 2007), the less powerful FR I radio sources are much more common than the FR II, and they are then (potentially) able to affect the general evolution of massive elliptical galaxies. Clearly, further simulations are required to quantitatively assess the effects of feedback in FR Is.

Acknowledgments

The numerical calculations were performed at CINECA in Bologna, Italy, with the help of an ISCRA grant. P.R. and G.B. acknowledge support by PRIN-INAF Grant, year 2014.

References

All Tables

All Figures

	Fig. 1 3D isocontours of the tracer distribution for case A (M = 4, η = 10-2, ℒkin = 1.1 × 1042 erg s-1) at t = 740 time units, 5.7 × 107 yr. The size of the computational box is 6.4 × 12 × 6.4 kpc.
In the text

Fig. 2 Results for case A (M = 4, η = 10-2, ℒkin = 1.1 × 1042 erg s-1) at t = 740 time units, 5.7 × 107 yr. Left: cut in the (y,z) plane of the logarithmic density distribution in units of jet density. The spatial units are in jet radii (100 pc) and the image extends over 12 kpc in the jet direction and over 6 kpc in the transverse direction. Right: maximum pressure on a (x,z) plane at a given position y along the jet as a function of y. The pressure is plotted, from now on, in units of the ambient value (instead of the computational units ${\mathit{\rho }}_{\mathrm{j}}{\mathit{c}}_{\mathrm{sj}}^{\mathrm{2}}$).

In the text

	Fig. 3 Results for a 2D simulation with the same jet parameters as case A and at t = 500 time units, 4 × 107 yr.
In the text

Fig. 4 Results for our reference case B (M = 4, η = 10-2, ℒkin ~ 1.1 × 1042 erg s-1), at t = 640 time units, 4.9 × 107 yr). The spatial units are in jet radii (100 pc) and the image extends over 12 kpc in the jet direction and over 6 kpc in the transverse direction. Cuts in the (y,z) plane of (top) the logarithmic density distribution in units of jet density, (center) tracer distribution, and (bottom) distribution of the longitudinal velocity of the external medium in units of csj.

In the text

	Fig. 5 Same as in Fig. 4, but for the pressure distribution of the external medium. The three phases of the evolution are labeled (top panel).The bottom panel shows the total pressure distribution in the initial part of the jet with internal shocks.
In the text

	Fig. 6 Cut in the (y,z) plane of the logarithmic density distribution for case C, differing from the reference case B for the value of α that describes the external medium, at t = 740 time units, 5.7 × 107 yr.
In the text

	Fig. 7 Cut in the (y,z) plane of the logarithmic density distribution for case D, overpressured jet, at t = 720 time units, 1.7 × 107 yr.
In the text

	Fig. 8 Cut in the (y,z) plane of the logarithmic density distribution for case E, the highest jet power (M = 40, η = 10-2, ℒkin = 1.1 × 1045), at t = 36 time units, 3 × 106 yr.
In the text

	Fig. 9 Maximum pressure on transverse planes vs y for case B at t = 640 time units (left panel), case E at t = 36 (right panel).
In the text

	Fig. 10 Left: cut in the (y,z) plane of the logarithmic density distribution for case F (M = 10, η = 10-2, ℒkin = 1.7 × 1043) at t = 98 time units, 7.5 × 106 yr; right: maximum pressure on transverse planes along the jet.
In the text

	Fig. 11 Left: cut in the (y,z) plane of the logarithmic density distribution for case G (M = 4, η = 10-3, ℒkin = 3.5 × 1042) at t = 1100 time units, 8.5 × 107 yr; right: maximum pressure on transverse planes along the jet.
In the text

	Fig. 12 Jet head position as a function of time (solid line) compared with the theoretical estimate for a uniform medium (dashed line) for cases B, E, and F in the left, center, and right panel, respectively.
In the text