11institutetext:  DiSAT, Università dell’Insubria, Via Valleggio 11, I-22100
Como, Italy
11email<EMAIL_ADDRESS>22institutetext: INAF – Osservatorio Astronomico
di Brera, via E. Bianchi 46, 23807 Merate, Italy 33institutetext: INAF –
Osservatorio Astrofisico di Torino, Strada Osservatorio 20, 10025 Pino
Torinese, Italy 44institutetext: Dipartimento di Fisica, Università degli
Studi di Torino, Via Pietro Giuria 1, 10125 Torino, Italy 55institutetext:
Dipartimento di Fisica, Universita‘ degli Studi di Genova, Via Dodecaneso 33,
I-16146 Genova, Italy 66institutetext: INAF - Istituto di Radioastronomia, Via
P. Gobetti 101, I-40129 Bologna, Italy 77institutetext: Dipartimento di Fisica
e Astronomia, Università di Bologna, Via P. Gobetti 93, I-40129 Bologna, Italy

# FR0 jets and recollimation-induced instabilities

A. Costa 1122 G. Bodo 33 F. Tavecchio 22 P. Rossi 33 A. Capetti 33 S.
Massaglia 44 A. Sciaccaluga 5522 R. D. Baldi 66 G. Giovannini 6677

###### Abstract

Context. The recently discovered population of faint FR0 radiogalaxies has
been interpreted as the extension to low power of the classical FRI sources.
Their radio emission appears to be concentrated in very compact (pc-scale)
cores, any extended emission is very weak or absent and VLBI observations show
that jets are already mildly or sub-relativistic at pc scales. Based on these
observational properties we propose here that the jets of FR0s are strongly
decelerated and disturbed at pc scale by hydrodynamical instabilities.

Aims. With the above scenario in mind, we study the dynamics of a low-power
relativistic jet propagating into a confining external medium, focusing on the
effects of entrainment and mixing promoted by the instabilities developing at
the jet-environment interface downstream of a recollimation shock.

Methods. We perform a 3D relativistic hydrodynamical simulation of a
recollimated jet by means of the state-of-the-art code PLUTO. The jet is
initially conical, relativistic (with initial Lorentz Factor $\Gamma$=5), cold
and light with respect to the confining medium, whose pressure is assumed to
slowly decline with distance. The magnetic field is assumed to be dynamically
unimportant.

Results. The 3D simulation shows that, after the first
recollimation/reflection shock system, a rapidly growing instability develops,
as a result of the interplay between recollimation-induced instabilities and
Richtmyer-Meshkov modes. In turn, the instabilities promote strong mixing and
entrainment that rapidly lead to the deceleration of the jet and spread its
momentum to slowly moving, highly turbulent external gas. We argue that this
mechanism could account for the peculiarities of the low-power FR0 jets. For
outflows with higher power, Lorentz factor or magnetic field, we expect that
the destabilizing effects are less effective, allowing the survival of the jet
up to the kpc scale, as observed in FRIs.

###### Key Words.:

galaxies: jets — radiation mechanisms: non-thermal — shock waves —-
instabilities

##  1 Introduction

Despite decades of efforts, the comprehension of relativistic jets ejected by
active galactic nuclei (AGN) is still rather sketchy (e.g. Blandford et al.,
2019). Challenges for current research concern the mechanisms able to
accelerate and collimate the jet in the central supermassive black hole (SMBH)
vicinity, the composition of the outflowing plasma (pair vs proton-electron,
matter vs magnetically dominated), the velocity structure of the flow, the
mechanism(s) accelerating particles to ultra-relativistic energies. Since
their discovery, one of the most active research topics has been the role of
the various instabilities in shaping the jet’s dynamical and dissipative
properties. A very broad range of instabilities can be at work: Kelvin-
Helmholtz, current driven, pressure driven, centrifugal, Rayleigh-Taylor (e.g.
Birkinshaw, 1996; Bodo et al., 2013, 2019; Kim et al., 2017, 2018; Begelman,
1998; Das & Begelman, 2019; Gourgouliatos & Komissarov, 2018), and these
different instabilities may lead to different outcomes in the jet dynamics and
energy dissipation. In particular, they are thought to play a prominent role
in low-power jets (classified as Fanaroff-Riley I (FRI) sources, Fanaroff &
Riley 1974), whose structure has been generally interpreted as the result of
instabilities and subsequent entrainment of external gas (e.g. Rossi et al.,
2008; Laing & Bridle, 2014; Perucho et al., 2014; Massaglia et al., 2016;
Rossi et al., 2020), eventually disrupting the jet at kpc scales. On the other
hand, powerful FRII jets seem to be much less prone to instabilities and can
reach distances up to Mpc scale (Willis et al., 1974), where they feed the
giant radio lobes.

It has been recently realized that, among local extragalactic jetted sources,
the largest fraction is composed by low-power objects, with radio morphology
characterized by a compact core with virtually no extended emission at kpc
scale (see Baldi 2023, and references therein). Since the radio properties
make this population the natural low-power extension of the FRI classical
radio galaxies (Fanaroff & Riley, 1974), they have been dubbed “FR0s”
(Ghisellini, 2011; Sadler et al., 2014). VLBI studies of FR0s (Cheng & An,
2018; Cheng et al., 2021; Baldi et al., 2021; Giovannini et al., 2023)
reported that the complex jet structure and the large number of two-sided
structures are strong evidence that, contrary to what is observed in FRI radio
galaxies, FR0 jets in VLBI images are mildly or sub-relativistic, with bulk
velocity on the order of $0.5c$ or less at parsec scales. Moreover many FR0
jets are complex and display substructures on pc scales, hinting for a strong
interaction with the surrounding interstellar medium. In turn, this suggests
that low-power jets of FR0s are possibly able to efficiently remove cold gas
from the nucleus of the host galaxy, thus influencing the accretion onto the
SMBH (Baldi, 2023). In fact, observational evidence continues to mount that
low-power jetted AGNs in general can deposit a large amount of jet energy in
the interstellar medium through shocks and turbulence (e.g., Venturi et al.
2021; Pereira-Santaella et al. 2022; Nandi et al. 2023). The lack of a
relevant extended emission and the morphology observed at VLBI scales indicate
that the deceleration and the dissipation of the FR0 jet power occur close to
the SMBH, at pc scale, therefore challenging any scenario involving mechanisms
operating in a smooth and gradual way. Furthermore, the compact core implies
the existence of a localized region of intense dissipation close to the AGN
core.

Building on this observational base, we adopt here a scenario for the dynamics
of FR0 jets whose main actor is a recollimation shock. Specifically, we assume
that a (weakly magnetized, low-power) jet expands and rapidly becomes
underpressured with respect to the external gas. In such situation a
recollimation shock structure develops into the jet (e.g. Komissarov & Falle,
1997; Bodo & Tavecchio, 2018), accompanied by growing instabilities that
decelerate and perturb the jet and quickly destroy the flow. In this framework
the shock plays the double role of ensuring both the dissipation (and the
subsequent emission) of part of the jet kinetic power, and the excitation of
instabilities which rapidly decelerate and disrupt the jet.

Recent simulations show that recollimation shocks promote the formation of
instabilities at the jet/external medium interface, that can have a strong
impact on the flow. Several studies, in particular, concentrate on the
Rayleigh-Taylor (RTI, Matsumoto & Masada 2013, 2019; Gottlieb et al. 2021),
centrifugal (CFI, induced by the effective gravity resulting from the motion
of the plasma along curved streamlines in the recollimation region
Gourgouliatos & Komissarov 2018) and Richtmyer-Meshkov (RMI; triggered by the
passage of the reflected shock at the jet/external medium interface, e.g.
Matsumoto & Masada 2013, 2019; Gottlieb et al. 2021) instabilities. Here we
will focus on the case of a “light” jet, i.e. stable against the RTI (e.g.
Abolmasov & Bromberg, 2023). Instabilities excited by recollimation can be
damped by a sufficiently intense magnetic field with a suitable geometry
(Matsumoto et al., 2021; Gottlieb et al., 2020), but in this preliminary
exploration we assume a pure hydrodynamical (HD) jet, a choice suitable to
model a flow in which the magnetic field is not dynamically important.

##  2 Simulation

As mentioned in the Introduction, we consider a scenario in which the jet,
after an initial phase of free expansion, becomes underpressured with respect
to the ambient medium. In such situation, if we assume axisymmetry, the jet is
characterized by a series of recollimation and reflection shocks (Komissarov &
Falle, 1997), as confirmed by two-dimensional cylindrical simulations (e.g.
Mizuno et al., 2015). Our aim is to investigate the stability of such
configuration, when the axisymmetry constraint is relaxed, in agreement with
the real jet images (see Boccardi et al., 2021, and references therein). To
this aim we first perform two-dimensional simulations starting with a conical
jet (as in Bodo & Tavecchio 2018) and let the system evolve until a steady
state is reached. We then use this steady state as initial condition for the
3D simulations.

The jet, whose initial opening angle is $\theta_{j}=0.2$, is relativistic,
with a Lorentz factor $\Gamma_{j}=5$ at injection, and propagates through a
surrounding isothermal medium at rest, with density and pressure that decay
along $z$ with power law profiles of index $\eta=0.5$. The jet is injected at
a distance $z_{0}$ from the cone vertex. We simulate a ”light” jet, which is
under-dense and under-pressured with respect to the confining gas. The values
of the density and pressure ratios between jet and ambient, at the jet base,
are respectively:

$\frac{\rho_{j,0}}{\rho_{ext,0}}=7.6\times
10^{-6}\,,\qquad\frac{p_{j,0}}{p_{ext,0}}=10^{-3}\,,$ (1)

and the external pressure is:

$\frac{p_{ext,0}}{\rho_{ext,0}\,c^{2}}=3\times 10^{-6},$ (2)

corresponding to a temperature of $3\times 10^{7}$K.

We perform the simulations with the relativistic hydrodynamical (RHD) module
of the state-of-the-art code PLUTO (Mignone et al., 2007). The computational
box covers the domain $[-5,5]\times[-5,5]\times[1,30]$ in units of $z_{0}$ and
we adopt a resolution of 35 points per initial jet radius $r_{0}=0.2z_{0}$. We
set outflow conditions at all boundaries except at $z=z_{0}$, where we fix the
injection of the jet for $r<r_{0}$ and the environment static profiles. We run
the simulation up to $t_{f}=368\,z_{0}/c$ (corresponding to $1840\,r_{0}/c$),
when a quasi steady-state is reached. We complement the RHD equations with the
equation for the evolution of a passive tracer, that is set unity for the
injected jet material and to zero for the ambient medium. In this way we can
study in detail the mixing process between jet and ambient. In our simulations
we adopt units so that $c=1$, $z_{0}=1$ and $\rho_{0,ext}=1$. The unit time
will be $t_{0}=z_{0}/c$. See the Appendix for more details on the numerical
setup.

##  3 Results

In the 3D simulations, relaxing the axisymmetry constraint, the dynamics is
strongly modified by instabilities that show a very fast growth. The global
structure of the perturbed jet is presented in figure 1, where we display a 3D
volume rendering of the $z$ component of the 4-velocity $\Gamma v_{z}$ at
$t_{f}$. At the base, the jet is relativistic ($\Gamma v_{z}>3$, yellow-red),
and it is possible to distinguish the first compression stage, followed by an
expansion and a second compression phase. The jet later decelerates while it
entrains external material ($\Gamma v_{z}\leq 1$, light blue), becoming
quickly sub-relativistic.

Figure 1: 3D volume rendering of the $z$ component of the jet 4-velocity . The
black-grey bar next to the colorbar shows which values of $\Gamma v_{z}$ are
opaque (grey), and which are transparent (black).

The jet instability seems to be ascribable to a combination of the Kelvin-
Helmholtz instability (KHI), the CFI, and the RMI instabilities. In the slow,
entrainment regions of Fig. 1, it is possible to find signatures of KHI-
induced helical deformation, but it is the RMI that develops important non-
linear perturbations.

More details on the dynamics can be obtained from Figures 2 and 3. In Fig. 2
we show two-dimensional slices, at $y=0$, of the distribution of the Lorentz
factor (right) and of the tracer of the external material, moving with
$v_{z}>0.1$ (left). The tracer value represents the fraction of external
material present in a given computational cell. The stationary axisymmetric
jet profile is over-plotted on the left in white with a contour of the jet
tracer. The four white horizontal lines on the right panel indicate the
locations of the $z=z^{*}$ cuts shown in Fig. 3, where we represent the $z$
component of the four-velocity $\Gamma v_{z}$.

Figure 2: The two pictures display the maps of the Lorentz factor (right
panel) and of the tracer of the medium material (left panel) that is
accelerated above a threshold of $v_{z}\geq 0.1$ in the $x-z$ plane, at $y=0$.
The white curve on the left is a contour of the 2D stationary jet. White
horizontal lines on the right indicate 4 different distances $z*$ at which we
evaluate the $x-y$ cuts in Fig. 3. Figure 3: The 4-velocity $\Gamma v_{z}$ in
the $x-y$ plane at different $z*=2.3,\,4.7,\,6.1,\,13$ respectively for A, B,
C, and D defined inn Fig. 2.

The jet is decelerated by the first strong recollimation shock (located at the
boundary of the yellow region in Fig. 2), that reaches the axis at $z\simeq
2.6$ and is reflected, reaching an anti-node at $z\simeq 3.8$. During the
expansion stage, after $z\simeq 3$, some external material starts to be
entrained at the jet-environment contact discontinuity (CD) (see Fig. 2),
because of small scale perturbations induced by the recollimation
instabilities. A secondary effect of these perturbations is that the
reflection shock becomes able to cross the CD where it is corrugated, and
excites the RMI (Matsumoto & Masada, 2013, 2019). The jet finally becomes
unstable to RMI approximately at the anti-node, as it can be seen in Fig. 3,
panel B, a region particularly favourable to the starting of instabilities
(Wilkes et al., 2008). After the anti-node, the jet is recollimated through a
second shock, but the downstream region is unstable, and after the second
recollimation point, at $z\simeq 5$, the jet is not able to expand again. Only
a portion of the material in the jet core continues to propagate at a
relativistic velocity, with $\Gamma\simeq 3$ up to $z\simeq 10$, while
entraining external medium through the turbulent interaction between the two
fluids. Finally, at $z>10$, the jet becomes sub-relativistic. These different
stages are also displayed in Fig. 3: panel A shows the transversal jet
structure before the first recollimation point, where we clearly distinguish
the unshocked jet portion in yellow and the shocked portion in orange; panels
B and C show the development of perturbations and the progressive deceleration
of the jet core surrounded by a slow mixing layer; in panel D we see that all
the jet has become sub-relativistic.

Figure 4: Upper panel: mass flux, $\Phi_{\rho}(z,t)=\int_{xy}\rho\Gamma
v_{z}dxdy$, starting from $\Phi_{\rho}(z_{0})=5.3\times 10^{-6}$, as function
of $z$ at different late times $t=[140,180,220,260,300,340,368]$. Lower panel:
averaged propagation velocity. More details can be found in the appendix.
Figure 5: The histogram displays the fraction of cells with a given jet
velocity $\Gamma v_{z}$ at the end of the simulation, calculated on three
$z=z*$ planes, representing the initial jet, the jet after the anti-node (as
in Fig. 3B) that has entrained material, and the jet further in $z$, where it
has become sub-relativistic (as in Fig. 3D). Figure 6: The picture shows the
$y=0$ slice of the logarithm of the pressure in code units, $\mbox{log}(p)$.
The colder regions, with $\mathcal{T}<0.1$, are not displayed.

Figure 4 compares the entrainment of external material with the overall jet
deceleration, by plotting the mass flux $\Phi_{\rho}(z,t)$ and the average
velocity $\langle v_{z}(z,t)\rangle$ in the $z$ direction, as functions of the
altitude $z$, where the different lines refer to different times, from black
to red, at $t_{f}$. The top panel shows $\Phi_{\rho}(z)=\int_{xy}\rho\Gamma
v_{z}\,dxdy$. For small values of $z$ the flux is small, since the jet density
is low ($\rho=7.6\times 10^{-6}$), but it starts increasing at $z\simeq 3$,
where there are the first signatures of entrainment. After the second
recollimation point, at $z\simeq 5$, the flux grows significantly, as a result
of the entrainment of the heavy external medium. The increase in the flux is
almost linear (Fig. 2). After increasing with time at the beginning of the
simulation, as shown from curves from $t=140$ to $t=260$, the mass flux seems
to converge to a stationary profile, showing little dispersion from $t=300$ to
$t_{f}$.

The bottom panel of Fig. 4 displays the propagation velocity averaged on $x-y$
planes. The strong decrease of $\langle v_{z}\rangle$ up to $z\simeq 3$ is due
to the recollimation shock (see also the red region in Fig. 3, panel A); then
the jet accelerates again as it expands. After the anti-node at $z\simeq 3.8$,
the average velocity decreases as a result of both the second recollimation
shock and of the starting of the entrainment process. After $z\simeq 5$, the
jet keeps slowing down, because part of its momentum is transferred to the
heavier entrained material, reaching sub-relativistic velocities in a smooth
way (see also Fig. 3, panel D). The curves for different times almost overlap,
with variations $<10\%$, indicating that the jet has reached an almost quasi
steady-state. The transition from relativistic to sub-relativistic velocities
is also displayed in Fig. 5, showing an histogram of $\Gamma v_{z}$, the $z$
component of the four-velocity distribution, in the $x-y$ plane, at three
different positions in $z$. Black refers to the jet base, blue refers to the
region just after the second recollimation point (see also Fig. 3, panel B),
where there still is a fast spine, and red to the region where the flow is
completely sub-relativistic (see also Fig. 3, panel D).

Fig. 6 shows a pressure map (in units of $\rho_{ext,0}c^{2}$, in logarithmic
scale) in the $x-z$ plane at $y=0$, at the end of the simulation. The cold
regions for which $\mathcal{T}=p/\rho c^{2}<0.1$, like the injection and the
environment, are left black in the image. Initially the jet is cold and its
energy is mainly kinetic but, when the fluid crosses the shocks, its thermal
energy increases, resulting in large values of temperature in the shock
downstream. After the first recollimation shock, the jet reaches pressure
balance with the cold environment ($\mathcal{T}_{ext,0}=3\times 10^{-6}$, from
Eq. 2); pressure is then further increased downstream of the reflection shock,
for $z>2.6$. After the expansion, we observe another highly pressurized region
after the second recollimation point at $z\simeq 5$. In these regions we
expect the presence of non-thermal particles accelerated at the shocks and
emitting synchrotron radiation. The synchrotron emissivity can be
qualitatively estimated as $j\propto pB^{2}$ (e.g. Bodo & Tavecchio, 2018),
where $p$ and $B$ are the pressure and the magnetic field strength in the jet.
Our simulation is purely hydrodynamical, hence we are not able to infer the
magnetic field strength. Assuming that the magnetic field energy density is a
fraction of the thermal pressure, the emissivity will scale as $j\propto
p^{2}$. We therefore expect that the regions downstream of the reflection
shocks are the brightest components, with an emissivity larger by a factor of
$10-100$ with respect to the other regions, which may originate much weaker
diffuse emission.

##  4 Discussion

We have reported the results of a simulation of a conical jet, underpressured
with respect to the environment, in which the development of an oblique
recollimation shock promotes the growth of instabilities. The instability
evolution leads to vigorous turbulent mixing with the external gas which is
entrained by the jet, and the consequent spreading of momentum results in a
rapid jet deceleration. Hydrodynamical simulations are scale invariant in
principle; here we set the length and density scales, that represent the jet
distance from the engine ($z_{0}$) and the external density ($\rho_{ext,0}$),
to values inferred from the typical properties of low-power radio galaxies
(e.g. Heckman & Best, 2014; Russell et al., 2015; Boccardi et al., 2021;
Casadio et al., 2021). Assuming $z_{0}=1$ pc and
$\rho_{ext,0}=1\,m_{p}\mbox{cm}^{-3}$ the simulation is suitable to reproduce
the observed properties of FR0 sources (e.g. Baldi, 2023). The simulated jet
is confined at a distance of $\sim 3$ pc, where it crosses powerful standing
shocks that are sites of particle acceleration and non-thermal emission,
creating bright spots that will result in the observed core at pc scales
(Baldi et al., 2021). The jet is mildly relativistic up to a distance of the
order of $10$ pc, transferring its momentum to entrained external material, to
finally become sub-relativistic on larger scales. This is in agreement with
the observed properties of FR0s. (e.g. Baldi et al., 2015, 2019; Cheng & An,
2018; Capetti et al., 2020; Cheng et al., 2021; Giovannini et al., 2023). The
injected jet luminosity is:

$\displaystyle L_{j}=$
$\displaystyle\pi(z_{0}\theta_{j})^{2}v_{z,j}\rho_{j,0}c^{2}h_{j,0}\Gamma^{2}$
(3) $\displaystyle\simeq$ $\displaystyle
10^{40}\left(\frac{\rho_{ext,0}}{1\,m_{p}\,\mbox{cm}^{-3}}\right)\left(\frac{z_{0}}{1\,\mbox{pc}}\right)^{2}\,\mbox{erg
s${}^{-1}$},$ (4)

where $h$ is the specific enthalpy and in the initially cold jet
$h_{j,0}\simeq 1$. This jet power is consistent with typical values derived
for FR0s (Baldi et al., 2018). Nevertheless, this simulation should be
interpreted as a case study of a light and cold jet. We expect lighter setups,
with $L_{j}\leq 10^{40}\,\mbox{erg s}^{-1}$ (thus describing low-power FR0
jets), to be even more unstable.

It is well known that the accretion flow and the environment play a key role
in determining the confinement properties of jets, and their transition to a
conical and cylindrical shape (e.g. Perucho & Martí, 2007; Park et al., 2023;
Rohoza et al., 2023). In particular, FR0s are hosted by giant elliptical
galaxies, with hints of a hot corona (Hardcastle et al., 2007; Torresi et al.,
2018) from X-ray spectroscopic data. Observations of low-power radio galaxies
suggest that jets propagate with a conical geometry at scales of $1-20$ pc in
a stratified medium (e.g. Russell et al., 2015; Boccardi et al., 2021; Casadio
et al., 2021), with evidences of recollimation shocks (for BL Lac Casadio et
al., 2021).

In terms of luminosity, FR0s can be interpreted as the low-power tail of the
FRI class of radio-galaxies, with $L\leq 10^{40}\,\mbox{erg s}^{-1}$ (Baldi et
al., 2018). In the scenario proposed here, the initial opening angle of the
jet, the jet-environment density and pressure ratios play the most important
roles in causing the confinement and triggering the instabilities downstream
of the shocks. We expect that FRIs, and especially FRIIs, are characterized by
more powerful jets, with higher density and/or magnetization, ensuring the
stability required to survive the recollimation instabilities and reach the
kpc scale. This is consistent with the view proposed in Gourgouliatos &
Komissarov (2018). The discussion of the role of the different parameters will
be addressed in a forthcoming paper (Costa et al. in prep.). More
fundamentally, these parameters could be connected to some of the key
properties of the central engine, some not directly observable at present,
such as the SMBH/accreting disk spin (Garofalo & Singh, 2019; Grandi et al.,
2021; Giovannini et al., 2023; Lalakos et al., 2023).

## Acknowledgments

We thank P. Coppi for useful discussions. We aknowledge financial support by a
INAF Theory Grant 2022 (PI F. Tavecchio) and the PRIN 2022 (2022C9TNNX)
project. We acknowledge support by CINECA, through ISCRA and Accordo Quadro
INAF-CINECA, and by PLEIADI, INAF - USC VIII, for the availability of HPC
resources. This work has been funded by the EU - Next Generation EU.

## References

  * Abolmasov & Bromberg (2023) Abolmasov, P. & Bromberg, O. 2023, MNRAS, 520, 3009
  * Baldi (2023) Baldi, R. D. 2023, A&A Rev., 31, 3
  * Baldi et al. (2015) Baldi, R. D., Capetti, A., & Giovannini, G. 2015, A&A, 576, A38
  * Baldi et al. (2019) Baldi, R. D., Capetti, A., & Giovannini, G. 2019, MNRAS, 482, 2294
  * Baldi et al. (2018) Baldi, R. D., Capetti, A., & Massaro, F. 2018, A&A, 609, A1
  * Baldi et al. (2021) Baldi, R. D., Giovannini, G., & Capetti, A. 2021, Galaxies, 9, 106
  * Begelman (1998) Begelman, M. C. 1998, ApJ, 493, 291
  * Birkinshaw (1996) Birkinshaw, M. 1996, Ap&SS, 242, 17
  * Blandford et al. (2019) Blandford, R., Meier, D., & Readhead, A. 2019, ARA&A, 57, 467
  * Boccardi et al. (2021) Boccardi, B., Perucho, M., Casadio, C., et al. 2021, A&A, 647, A67
  * Bodo et al. (2013) Bodo, G., Mamatsashvili, G., Rossi, P., & Mignone, A. 2013, MNRAS, 434, 3030
  * Bodo et al. (2019) Bodo, G., Mamatsashvili, G., Rossi, P., & Mignone, A. 2019, MNRAS, 485, 2909
  * Bodo & Tavecchio (2018) Bodo, G. & Tavecchio, F. 2018, A&A, 609, A122
  * Capetti et al. (2020) Capetti, A., Massaro, F., & Baldi, R. D. 2020, A&A, 633, A161
  * Casadio et al. (2021) Casadio, C., MacDonald, N. R., Boccardi, B., et al. 2021, A&A, 649, A153
  * Cheng et al. (2021) Cheng, X., An, T., Sohn, B. W., Hong, X., & Wang, A. 2021, MNRAS, 506, 1609
  * Cheng & An (2018) Cheng, X. P. & An, T. 2018, ApJ, 863, 155
  * Das & Begelman (2019) Das, U. & Begelman, M. C. 2019, MNRAS, 482, 2107
  * Fanaroff & Riley (1974) Fanaroff, B. L. & Riley, J. M. 1974, MNRAS, 167, 31P
  * Garofalo & Singh (2019) Garofalo, D. & Singh, C. B. 2019, ApJ, 871, 259
  * Ghisellini (2011) Ghisellini, G. 2011, in American Institute of Physics Conference Series, Vol. 1381, American Institute of Physics Conference Series, ed. F. A. Aharonian, W. Hofmann, & F. M. Rieger, 180–198
  * Giovannini et al. (2023) Giovannini, G., Baldi, R. D., Capetti, A., Giroletti, M., & Lico, R. 2023, A&A, 672, A104
  * Gottlieb et al. (2020) Gottlieb, O., Bromberg, O., Singh, C. B., & Nakar, E. 2020, MNRAS, 498, 3320
  * Gottlieb et al. (2021) Gottlieb, O., Nakar, E., & Bromberg, O. 2021, MNRAS, 500, 3511
  * Gourgouliatos & Komissarov (2018) Gourgouliatos, K. N. & Komissarov, S. S. 2018, Nature Astronomy, 2, 167
  * Grandi et al. (2021) Grandi, P., Torresi, E., Macconi, D., Boccardi, B., & Capetti, A. 2021, ApJ, 911, 17
  * Hardcastle et al. (2007) Hardcastle, M. J., Evans, D. A., & Croston, J. H. 2007, MNRAS, 376, 1849
  * Heckman & Best (2014) Heckman, T. M. & Best, P. N. 2014, ARA&A, 52, 589
  * Kim et al. (2017) Kim, J., Balsara, D. S., Lyutikov, M., & Komissarov, S. S. 2017, MNRAS, 467, 4647
  * Kim et al. (2018) Kim, J., Balsara, D. S., Lyutikov, M., & Komissarov, S. S. 2018, MNRAS, 474, 3954
  * Komissarov & Falle (1997) Komissarov, S. S. & Falle, S. A. E. G. 1997, MNRAS, 288, 833
  * Laing & Bridle (2014) Laing, R. A. & Bridle, A. H. 2014, MNRAS, 437, 3405
  * Lalakos et al. (2023) Lalakos, A., Tchekhovskoy, A., Bromberg, O., et al. 2023, arXiv e-prints, arXiv:2310.11487
  * Massaglia et al. (2016) Massaglia, S., Bodo, G., Rossi, P., Capetti, S., & Mignone, A. 2016, A&A, 596, A12, Mas16
  * Matsumoto et al. (2021) Matsumoto, J., Komissarov, S. S., & Gourgouliatos, K. N. 2021, MNRAS, 503, 4918
  * Matsumoto & Masada (2013) Matsumoto, J. & Masada, Y. 2013, ApJ, 772, L1
  * Matsumoto & Masada (2019) Matsumoto, J. & Masada, Y. 2019, MNRAS, 490, 4271
  * Mignone & Bodo (2005) Mignone, A. & Bodo, G. 2005, MNRAS, 364, 126
  * Mignone et al. (2007) Mignone, A., Bodo, G., Massaglia, S., et al. 2007, in JENAM-2007, “Our Non-Stable Universe”, 96–96
  * Mignone et al. (2005) Mignone, A., Plewa, T., & Bodo, G. 2005, ApJS, 160, 199
  * Mizuno et al. (2015) Mizuno, Y., Gómez, J. L., Nishikawa, K.-I., et al. 2015, ApJ, 809, 38
  * Mukherjee et al. (2020) Mukherjee, D., Bodo, G., Mignone, A., Rossi, P., & Vaidya, B. 2020, MNRAS, 499, 681
  * Nandi et al. (2023) Nandi, P., Stalin, C. S., Saikia, D. J., et al. 2023, arXiv e-prints, arXiv:2307.13263
  * Park et al. (2023) Park, J., Kino, M., Nagai, H., et al. 2023, arXiv e-prints, arXiv:2311.08647
  * Pereira-Santaella et al. (2022) Pereira-Santaella, M., Álvarez-Márquez, J., García-Bernete, I., et al. 2022, A&A, 665, L11
  * Perucho & Martí (2007) Perucho, M. & Martí, J. M. 2007, MNRAS, 382, 526
  * Perucho et al. (2014) Perucho, M., Martí, J. M., Laing, R. A., & Hardee, P. E. 2014, MNRAS, 441, 1488
  * Rohoza et al. (2023) Rohoza, V., Lalakos, A., Paik, M., et al. 2023, arXiv e-prints, arXiv:2311.00018
  * Rossi et al. (2020) Rossi, P., Bodo, G., Massaglia, S., & Capetti, A. 2020, A&A, 642, A69
  * Rossi et al. (2008) Rossi, P., Mignone, A., Bodo, G., Massaglia, S., & Ferrari, A. 2008, A&A, 488, 795
  * Russell et al. (2015) Russell, H. R., Fabian, A. C., McNamara, B. R., & Broderick, A. E. 2015, MNRAS, 451, 588
  * Sadler et al. (2014) Sadler, E. M., Ekers, R. D., Mahony, E. K., Mauch, T., & Murphy, T. 2014, MNRAS, 438, 796
  * Torresi et al. (2018) Torresi, E., Grandi, P., Capetti, A., Baldi, R. D., & Giovannini, G. 2018, MNRAS, 476, 5535
  * Venturi et al. (2021) Venturi, G., Cresci, G., Marconi, A., et al. 2021, A&A, 648, A17
  * Wilkes et al. (2008) Wilkes, J. I., Danehy, P., Nowak, R., & Alderfer, D. 2008, in 38th Fluid Dynamics Conference and Exhibit, AIAA, 4389
  * Willis et al. (1974) Willis, A. G., Strom, R. G., & Wilson, A. S. 1974, Nature, 250, 625

##  Appendix A Numerical setup

The simulations presented in this letter are performed with the PLUTO (Mignone
et al. 2007) modular code that solves the set of conservation equations of
fluid dynamics. In particular we employed the RHD module, that solves the
system of relativistic hydrodynamics equations

$\frac{\partial}{\partial t}\left(\begin{matrix}\rho\Gamma\\\ \rho
h\Gamma^{2}\mathbf{v}\\\ \rho h\Gamma^{2}-p\\\ \rho\Gamma
f\end{matrix}\right)+\mathbf{\nabla}\cdot\left(\begin{matrix}\rho\Gamma\mathbf{v}\\\
\rho h\Gamma^{2}\mathbf{v}\mathbf{v}+p\mathbf{I}\\\ \rho
h\Gamma^{2}\mathbf{v}\\\ \rho\Gamma
f\mathbf{v}\end{matrix}\right)=\left(\begin{matrix}0\\\ \mathbf{f}_{g}\\\
\mathbf{f}_{g}\cdot\mathbf{v}\\\ 0\end{matrix}\right),$ (5)

where $\rho,\,\Gamma,\,h,\,\mathbf{v},\,p,$ respectively represent the rest-
frame number density, the Lorentz factor, the proper specific enthalpy, the
3-velocity in units of $c$, and the pressure of the fluid. The acceleration
term $\mathbf{f}_{g}$ is the specific external force three-vector, set to
$\mathbf{f}_{g}=\nabla p_{ext}(t=0)$, in order to maintain in dynamical
equilibrium the static ambient medium. We also evolve a passive tracer $f$,
initially set to $0$ for the surrounding medium, and to $1$ for the injected
relativistic jet, to track the evolution of the jet material and to study the
mixing between the two fluids. We close the set of equations with the Taub-
Matthews equation of state:

$h=\frac{5}{2}\mathcal{T}+\sqrt{\frac{9}{4}\mathcal{T}^{2}+1},$ (6)

with $\mathcal{T}=p/(\rho c^{2})$, that approximates the Synge EoS of a
single-species relativistic perfect fluid (Mignone et al. 2005). We adopt a
linear reconstruction scheme with a second order Runge-Kutta method for time
integration, and the HLLC Riemann solver (Mignone & Bodo 2005).

###  A.1 2D setup

The preliminary axisymmetric simulation is run in 2D cylindrical coordinates
$(r,z)$ in a domain $[0,20]\times[0.5,30]$, where the lengths are in units of
$z_{0}$, that represents the distance from the jet launching site. The grid is
uniform with $1000\times 3000$ points in $[0,1.5]\times[0.5,20]$ and
geometrically stretched with $400\times 700$ grid points in the outer regions.

As initial condition, at $t=0$, in the region $r/z<0.2$, we have an
axisymmetric conical outflow, with opening angle $\theta_{j}=0.2$, with
Lorentz factor $\Gamma_{j}=5$, more precisely the velocity components are
defined as

$v_{z}=\sqrt{1-\frac{1}{\Gamma_{j}^{2}}}\frac{z}{R},$ (7)
$v_{r}=\sqrt{1-\frac{1}{\Gamma_{j}^{2}}}\frac{r}{R},$ (8)

where $R$ is the spherical radius defined as $R=\sqrt{r^{2}+z^{2}}$. Density
and pressure decay adiabatically as

$\rho_{j}(r,z,t=0)=\rho_{j}(0,z_{0},0)\left(\frac{R}{R_{0}}\right)^{-2},$ (9)

and pressure

$p_{j}(r,z,t=0)=p_{j}(0,z_{0},0)\left(\frac{R}{R_{0}}\right)^{-2\gamma},$ (10)

where $\gamma$ is the adiabatic index derived from the EoS, which in the case
of a cold gas yields the classical $\gamma=5/3$. For $r/z>0.2$ there is a
static medium, whose density $\rho_{ext}$ and pressure are power law functions
of the altitude $z$:

$\displaystyle\rho_{ext}(z,t=0)=\rho_{ext}(z_{0},0)\left(\frac{z}{z_{0}}\right)^{-\eta},$
(11) $\displaystyle
p_{ext}(z,t=0)=p_{ext}(z_{0},0)\left(\frac{z}{z_{0}}\right)^{-\eta},$ (12)

where the power law index is $\eta=0.5$. The exact value of $\eta$ is not
decisive, as long as it’s smaller than $2$, so that the higher pressure of the
external medium confines the relativistic jet.

We simulate a cold jet, that is light compared to the environment, so we
define its density and pressure with respect to the external values, as

$\frac{\rho_{j}(0,z_{0},0)}{\rho_{ext}(z_{0},0)}=7.6\times 10^{-6},\\\
\frac{p_{j}(0,z_{0},0)}{p_{ext}(z_{0},0)}=10^{-3},$ (13)

while the external pressure is defined through the temperature:
$\mathcal{T}_{ext}(z_{0},0)=3\times 10^{-6}$. In our simulations we set
$\rho_{ext,0}$ as the density unit.

In order to avoid numerical noise at the contact discontinuity, the initial
condition is smoothed at the jet-environment boundary. We smoothed the Lorentz
factor, the density and the pressure, with functions of the type

$q=q_{ext}+\left(q_{j}-q_{ext}\right)\operatorname{sech}\left[\left(\frac{r}{z\theta_{q}}\right)^{\alpha_{q}}\right]$
(14)

in the inlet regions (Mukherjee et al. 2020). The index $\alpha_{q}$ and the
angle $\theta_{q}$, that define the width and the radial scale of the
smoothing, depend on the specific quantity $q$, and are different for
different values to avoid artificial local extrema in the energy and/or
momentum (Abolmasov & Bromberg 2023). We set $\theta_{\Gamma}=0.16$ and
$\alpha_{\Gamma}=8$ for the Lorentz factor, and $\theta_{\rho,p}=0.29$,
$\alpha_{\rho,p}=10$ for the density and the pressure.

We use outflow conditions at the right ($r=20$) and upper ($z=30$) boundaries,
while we use reflective conditions at the axis $r=0$. At the lower boundary
$z=z_{L}=0.5$, for $r>\theta_{j}z_{L}$, we extend in the ghost cells the
initial analytical pressure and density profiles, for ensuring dynamical
equilibrium. We keep these values constant. The jet is injected in the region
defined by $0.5<z<1$ and $0<r<\theta_{j}z$, where we set the velocities,
density and pressure defined in Eqs. 7, 8, 9 and 10. In addition the fluxes of
the Riemann solver are set to zero and hence the fluid variables remain
unchanged in this region. This is done to avoid spurious effects at the lower
boundary.

The 2D case is evolved until a steady state is reached, up to $t_{f}=3000$ in
units of $z_{0}/c$.

###  A.2 3D setup

The initial condition for the 3D simulation is the axisymmetric steady-state
reached in the 2D simulations described above. The 2D results, obtained in the
cylindrical coordinates $(r,z)$, are projected on the Cartesian coordinates
$(x,y,z)$. The computational domain is made of $550\times 550\times 1850$ grid
points, covering the physical domain $[-5,5]\times[-5,5]\times[1,30]$. Like in
the 2D case, the grid is uniform only in the inner region,
$[-1,1]\times[-1,1],\times[1,20]$, with $300\times 300\times 1500$ cells, and
stretched outside. We use outflow conditions at all boundaries, except at the
lower one. In the 3D simulation we set the lower boundary at $z=z_{0}=1$. For
$\sqrt{x^{2}+y^{2}}<\theta_{j}\,z_{0}$ we have injection conditions with the
jet parameters, and for $\sqrt{x^{2}+y^{2}}>\theta_{j}\,z_{0}$ we set density
and pressure constant in time and following the initial profiles, as we have
done in 2D.

###  A.3 Calculation of the average of the propagation velocity

In order to average the velocity of the moving material, plotted in 4, we
defined a step function $g(x,y,z)$ as

$g(x,y,z,t)=\begin{cases}0\quad\mbox{if }v_{z}(x,y,z,t)<0.1\\\ 1\quad\mbox{if
}v_{z}(x,y,z,t)\geq 0.1\end{cases}$ (15)

that selects the jet section with a threshold on the propagation velocity
$v_{z}$. The average velocity $\langle v_{z}\rangle$ is then the average of
the distribution of the velocity on the $x-y$ jet section, calculated as

$\langle v_{z}(z,t)\rangle=A^{-1}\int_{xy}v_{z}(x,y,z,t)g(x,y,z,t)dxdy\,,$
(16)

where $A=\int_{xy}g(x,y,z,t)dxdy$.