FR0 jets and recollimation-induced instabilities

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.}{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. 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. The 3D simulation shows that, after the first recollimation/reflection shock system, a rapidlygrowing 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.


Introduction
Despite decades of efforts, the comprehension of relativistic jets ejected by active galactic nuclei (AGNs) is still rather unclear (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 versus proton-electron, matter versus magnetically dominated), the velocity structure of the flow, and 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, and Rayleigh-Taylor (e.g., Birkinshaw 1996;Bodo et al. 2013Bodo et al. , 2019;;Kim et al. 2017Kim et al. , 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 kiloparsec scales.On the other hand, powerful FRII jets seem to be much less prone to instabilities and can reach distances up to the megaparsec 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 are composed by low-power objects, with a radio morphology characterized by a compact core with virtually no extended emission at the kiloparsec 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).with very-long-baseline interferometry (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 a bulk velocity on the order of 0.5c or less at parsec scales.Moreover many FR0 jets are complex and display substructures on parsec 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 the parsec 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, for this work we adopt 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 a situation, a recollimation shock structure develops in 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 and external medium interface, which can have a strong impact on the flow.Several studies, in particular, concentrate on the Rayleigh-Taylor instability (RTI, Matsumoto & Masada 2013, 2019;Gottlieb et al. 2021), the centrifugal instability (CFI, induced by the effective gravity resulting from the motion of the plasma along curved streamlines in the recollimation region; Gourgouliatos & Komissarov 2018), and the Richtmyer-Meshkov instability (RMI, triggered by the passage of the reflected shock at the jet and external medium interface, e.g., Matsumoto & Masada 2013, 2019;Gottlieb et al. 2021).Here we focus on the case of a "light" jet: 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); however, 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.

Simulation
As mentioned in the Introduction, we consider a scenario in which the jet, after an initial phase of free expansion, becomes under-pressured with respect to the ambient medium.In such a situation, if we assume axisymmetry, the jet is characterized by a series of recollimation and reflection shocks (Komissarov & Falle 1997), as confirmed by 2D cylindrical simulations (e.g., Mizuno et al. 2015).Our aim is to investigate the stability of such a 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 performed 2D simulations starting with a conical jet (as in Bodo & Tavecchio 2018) and let the system evolve until a steady state was reached.We then used this steady state as an initial condition for the 3D simulations.
The jet, whose initial opening angle is θ j = 0.2, is relativistic, with a Lorentz factor Γ j = 5 at injection, and it propagates through a surrounding isothermal medium at rest, with a density and pressure that decay along z with power law profiles with an index η = 0.5.The jet was injected at a distance z 0 from the cone vertex.We simulated 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, and the external pressure is corresponding to a temperature of 3 × 10 7 K.
We performed 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] × [−5, 5] × [1, 30] in units of z 0 and we adopted 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 fixed the injection of the jet for r < r 0 and the environment static profiles.We ran the simulation up to t f = 368 z 0 /c (corresponding to 1840 r 0 /c), when a quasi steady state was reached.We complemented the RHD equations with the equation for the evolution of a passive tracer, which was set to unity for the injected jet material and to zero for the ambient medium.In this way we could study in detail the mixing process between jet and ambient.In our simulations, we adopted units so that c = 1, z 0 = 1, and ρ 0,ext = 1.The unit time is t 0 = z 0 /c.Readers can refer to the appendix for more details on the numerical setup.

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 Fig. 1, where we display a 3D volume rendering of the z component of the four-velocity Γv z at t f .At the base, the jet is relativistic (Γ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 (Γv z ≤ 1, light blue), quickly becoming subrelativistic.The jet instability seems to be ascribable to a combination of the Kelvin-Helmholtz instability (KHI), the CFI, and the RMI.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 nonlinear perturbations.
More details on the dynamics can be obtained from Figs. 2 and 3.In Fig. 2 we show 2D 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 L19, page 2 of 8  The jet is decelerated by the first strong recollimation shock (located at the boundary of the yellow region in Fig. 2), which reaches the axis at z 2.6 and is reflected, reaching an antinode at z 3.8.During the expansion stage, after z 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 antinode, as it can be seen in Fig. 3 (panel B), a region particularly favorable to the starting of instabilities (Wilkes et al. 2008).After the antinode, the jet is recollimated through a second shock, but the downstream region is unstable, and after the second recollimation point, at z 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 Γ 3 up to z 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 of the jet has become sub-relativistic.
Figure 4 compares the entrainment of external material with the overall jet deceleration, by plotting the mass flux Φ ρ (z, t) and the average velocity v z (z, t) 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 Φ ρ (z) = xy ρΓv z dxdy.For small values of z, the flux is small since the jet density is low (ρ = 7.6 × 10 −6 ); however, it starts increasing at z 3, where there are the first signatures of entrainment.After the second recollimation point, at z 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 .L19, page 3 of 8  [140,180,220,260,300,340,368].Lower panel: averaged propagation velocity.More details can be found in the appendix.
The bottom panel of Fig. 4 displays the propagation velocity averaged on x−y planes.The strong decrease in v z up to z 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 antinode at z 3.8, the average velocity decreases as a result of both the second recollimation shock and of the start of the entrainment process.After z 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 a histogram of Γv z , the z component of the fourvelocity distribution, in the x − y plane, at three different positions in z.Black refers to the jet base, blue corresponds to the region just after the second recollimation point (see also Fig. 3, panel B), where there still is a fast spine, and red represents the region where the flow is completely sub-relativistic (see also Fig. 3, panel D).
Figure 6 shows a pressure map (in units of ρ 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 T = p/ρc 2 < 0.1, such as the injection and the environment, are left black in the image.Initially the jet is cold and its energy is mainly kinetic; however, Fig. 5. Histogram displaying the fraction of cells with a given jet velocity Γv z at the end of the simulation, calculated on three z = z * planes, representing the initial jet, the jet after the antinode (as in Fig. 3B) that has entrained material, and the jet further in z, where it has become sub-relativistic (as in Fig. 3D).
when the fluid crosses the shocks, its thermal energy increases, resulting in large values for the temperature in the shock downstream.After the first recollimation shock, the jet reaches pressure balance with the cold environment (T ext,0 = 3 × 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 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 ∝ 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 scales as j ∝ 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 instead may cause much weaker diffuse emission.

Discussion
We have reported the results of a simulation of a conical jet, under-pressured 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 that 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, which represent the jet distance from the engine (z 0 ) and the external density (ρ 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 ρ ext,0 = 1 m p 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 L19, page 4 of 8 Fig. 6.Slice in the y = 0 plane of the logarithm of the pressure in code units, log(p).The colder regions, with T < 0.1, are not displayed.of ∼3 pc, where it crosses powerful standing shocks that are sites of particle acceleration and nonthermal emission, creating bright spots that result in the observed core at parsec scales (Baldi et al. 2021).The jet is mildly relativistic up to a distance on the order of 10 pc, transferring its momentum to entrained external material and it finally becomes sub-relativistic on larger scales.This is in agreement with the observed properties of FR0s (e.g., Baldi et al. 2015Baldi et al. , 2019;;Cheng & An 2018;Capetti et al. 2020;Cheng et al. 2021;Giovannini et al. 2023).The injected jet luminosity is: where h is the specific enthalpy and in the initially cold jet h j,0 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 ≤ 10 40 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. 2024;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 evidence 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 ≤ 10 40 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 a higher density and/or magnetization, ensuring the stability required to survive the recollimation instabilities and reach the kiloparsec 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 are not directly observable at the present time, such as the SMBH and accreting disk spin (Garofalo & Singh 2019;Grandi et al. 2021;Giovannini et al. 2023;Lalakos et al. 2023).

Appendix A: Numerical setup
The simulations presented in this Letter were 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, which solves the system of relativistic hydrodynamics equations where ρ, Γ, h, v, and p represent the rest-frame number density, the Lorentz factor, the proper specific enthalpy, the threevelocity in units of c, and the pressure of the fluid, respectively.The acceleration term f g is the specific external force three vector, set to f g = ∇p ext (t = 0), in order to maintain the static ambient medium in dynamical equilibrium.We also evolved a passive tracer f , initially set to zero for the surrounding medium, and to one for the injected relativistic jet, to track the evolution of the jet material and to study the mixing between the two fluids.We closed the set of equations with the Taub-Matthews equation of state (EoS)

Fig. 1 .
Fig. 1. 3D volume rendering of the z component of the jet four-velocity.The black-gray bar next to the colorbar shows which values of Γv z are opaque (gray), and which are transparent (black).

Fig. 2 .
Fig. 2. Two pictures displaying 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 ≥ 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 four different distances z * at which we evaluate the x−y cuts in Fig. 3.