EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 682 (February 2024) A&A, 682 (2024) L19 Full HTML
Open Access
Issue
A&A
Volume 682, February 2024
Article Number L19
Number of page(s) 8
Section Letters to the Editor
DOI https://doi.org/10.1051/0004-6361/202348954
Published online 26 February 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. 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. 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 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 under-pressured 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.

2. 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 z0 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

ρ j , 0 ρ ext , 0 = 7.6 × 10 − 6 , p j , 0 p ext , 0 = 10 − 3 , $$ \begin{aligned} \frac{\rho _{j,0}}{\rho _{\rm ext,0}} = 7.6 \times 10^{-6} , \qquad \frac{p_{j,0}}{p_{\rm ext,0}} = 10^{-3}, \end{aligned} $$(1)

respectively, and the external pressure is

p ext , 0 ρ ext , 0 c 2 = 3 × 10 − 6 , $$ \begin{aligned} \frac{p_{\rm ext,0}}{\rho _{\rm ext,0} \, c^2}=3 \times 10^{-6}, \end{aligned} $$(2)

corresponding to a temperature of 3 × 107 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 z0 and we adopted a resolution of 35 points per initial jet radius r0 = 0.2z0. We set outflow conditions at all boundaries except at z = z0, where we fixed the injection of the jet for r < r0 and the environment static profiles. We ran the simulation up to tf = 368 z0/c (corresponding to 1840 r0/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, z0 = 1, and ρ0,ext = 1. The unit time is t0 = z0/c. Readers can refer to 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 Fig. 1, where we display a 3D volume rendering of the z component of the four-velocity Γvz at tf. At the base, the jet is relativistic (Γvz > 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 (Γvz ≤ 1, light blue), quickly becoming sub-relativistic. 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.

thumbnail 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 Γvz are opaque (gray), and which are transparent (black).

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 vz > 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 in 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 Γvz.

thumbnail 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 vz ≥ 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.

thumbnail Fig. 3.

Four-velocity Γvz in the x − y plane at different z* = 2.3,  4.7,  6.1,   and 13 for A, B, C, and D, respectively, defined in Fig. 2.

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 ⟨vz(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 tf. The top panel shows Φρ(z) = ∫xyρΓvz 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 tf.

thumbnail Fig. 4.

Entertainment and deceleration as functions of time and altitude. Upper panel: mass flux, Φρ(z, t) = ∫xyρΓvzdxdy, starting from Φρ(z0) = 5.3 × 10−6, as a 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.

The bottom panel of Fig. 4 displays the propagation velocity averaged on x − y planes. The strong decrease in ⟨vz⟩ 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 Γvz, 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 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).

thumbnail Fig. 5.

Histogram displaying the fraction of cells with a given jet velocity Γvz 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).

Figure 6 shows a pressure map (in units of ρext, 0c2, in logarithmic scale) in the x − z plane at y = 0, at the end of the simulation. The cold regions for which 𝒯 = p/ρc2 < 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, 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 (𝒯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 ∝ pB2 (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 ∝ p2. 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.

thumbnail Fig. 6.

Slice in the y = 0 plane of the logarithm of the pressure in code units, log(p). The colder regions, with 𝒯 < 0.1, are not displayed.

4. 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 (z0) 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 z0 = 1 pc and ρ0,ext = 1 mp 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 ∼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. 2015, 2019; Cheng & An 2018; Capetti et al. 2020; Cheng et al. 2021; Giovannini et al. 2023). The injected jet luminosity is:

L j = π ( z 0 θ j ) 2 v z , j ρ j , 0 c 2 h j , 0 Γ 2 $$ \begin{aligned}&L_j = \pi (z_0 \theta _j)^2 v_{z,j} {\rho }_{j,0} c^2h_{j,0} \Gamma ^2 \end{aligned} $$(3)

≃ 10 40 ( ρ ext , 0 1 m p cm − 3 ) ( z 0 1 pc ) 2 e r g s − 1 , $$ \begin{aligned}&\quad \simeq 10^{40}\left(\frac{{\rho }_{\rm ext,0}}{1\,m_{\rm p}\,\mathrm{cm} ^{-3}}\right)\left(\frac{z_0}{1\,\mathrm{pc} }\right)^2\,\mathrm erg\,s ^{-1}, \end{aligned} $$(4)

where h is the specific enthalpy and in the initially cold jet hj, 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 Lj ≤ 1040 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 ≤ 1040 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).

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

  1. Abolmasov, P., & Bromberg, O. 2023, MNRAS, 520, 3009 [NASA ADS] [CrossRef] [Google Scholar]
  2. Baldi, R. D. 2023, A&ARv, 31, 3 [NASA ADS] [CrossRef] [Google Scholar]
  3. Baldi, R. D., Capetti, A., & Giovannini, G. 2015, A&A, 576, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Baldi, R. D., Capetti, A., & Massaro, F. 2018, A&A, 609, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Baldi, R. D., Capetti, A., & Giovannini, G. 2019, MNRAS, 482, 2294 [NASA ADS] [CrossRef] [Google Scholar]
  6. Baldi, R. D., Giovannini, G., & Capetti, A. 2021, Galaxies, 9, 106 [NASA ADS] [CrossRef] [Google Scholar]
  7. Begelman, M. C. 1998, ApJ, 493, 291 [NASA ADS] [CrossRef] [Google Scholar]
  8. Birkinshaw, M. 1996, Ap&SS, 242, 17 [NASA ADS] [CrossRef] [Google Scholar]
  9. Blandford, R., Meier, D., & Readhead, A. 2019, ARA&A, 57, 467 [NASA ADS] [CrossRef] [Google Scholar]
  10. Boccardi, B., Perucho, M., Casadio, C., et al. 2021, A&A, 647, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Bodo, G., Mamatsashvili, G., Rossi, P., & Mignone, A. 2013, MNRAS, 434, 3030 [Google Scholar]
  12. Bodo, G., Mamatsashvili, G., Rossi, P., & Mignone, A. 2019, MNRAS, 485, 2909 [Google Scholar]
  13. Bodo, G., & Tavecchio, F. 2018, A&A, 609, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Capetti, A., Massaro, F., & Baldi, R. D. 2020, A&A, 633, A161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Casadio, C., MacDonald, N. R., Boccardi, B., et al. 2021, A&A, 649, A153 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Cheng, X. P., & An, T. 2018, ApJ, 863, 155 [NASA ADS] [CrossRef] [Google Scholar]
  17. Cheng, X., An, T., Sohn, B. W., Hong, X., & Wang, A. 2021, MNRAS, 506, 1609 [NASA ADS] [CrossRef] [Google Scholar]
  18. Das, U., & Begelman, M. C. 2019, MNRAS, 482, 2107 [NASA ADS] [CrossRef] [Google Scholar]
  19. Fanaroff, B. L., & Riley, J. M. 1974, MNRAS, 167, 31P [Google Scholar]
  20. Garofalo, D., & Singh, C. B. 2019, ApJ, 871, 259 [NASA ADS] [CrossRef] [Google Scholar]
  21. Ghisellini, G. 2011, in 25th Texas Symposium on Relativistic Astrophysics (Texas 2010), eds. F. A. Aharonian, W. Hofmann, & F. M. Rieger, AIP Conf. Ser., 1381, 180 [Google Scholar]
  22. Giovannini, G., Baldi, R. D., Capetti, A., Giroletti, M., & Lico, R. 2023, A&A, 672, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Gottlieb, O., Bromberg, O., Singh, C. B., & Nakar, E. 2020, MNRAS, 498, 3320 [NASA ADS] [CrossRef] [Google Scholar]
  24. Gottlieb, O., Nakar, E., & Bromberg, O. 2021, MNRAS, 500, 3511 [Google Scholar]
  25. Gourgouliatos, K. N., & Komissarov, S. S. 2018, Nat. Astron., 2, 167 [Google Scholar]
  26. Grandi, P., Torresi, E., Macconi, D., Boccardi, B., & Capetti, A. 2021, ApJ, 911, 17 [NASA ADS] [CrossRef] [Google Scholar]
  27. Hardcastle, M. J., Evans, D. A., & Croston, J. H. 2007, MNRAS, 376, 1849 [NASA ADS] [CrossRef] [Google Scholar]
  28. Heckman, T. M., & Best, P. N. 2014, ARA&A, 52, 589 [Google Scholar]
  29. Kim, J., Balsara, D. S., Lyutikov, M., & Komissarov, S. S. 2017, MNRAS, 467, 4647 [Google Scholar]
  30. Kim, J., Balsara, D. S., Lyutikov, M., & Komissarov, S. S. 2018, MNRAS, 474, 3954 [Google Scholar]
  31. Komissarov, S. S., & Falle, S. A. E. G. 1997, MNRAS, 288, 833 [Google Scholar]
  32. Laing, R. A., & Bridle, A. H. 2014, MNRAS, 437, 3405 [NASA ADS] [CrossRef] [Google Scholar]
  33. Lalakos, A., Tchekhovskoy, A., Bromberg, O., et al. 2023, arXiv e-prints [arXiv:2310.11487] [Google Scholar]
  34. Massaglia, S., Bodo, G., Rossi, P., Capetti, S., & Mignone, A. 2016, A&A, 596, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Matsumoto, J., & Masada, Y. 2013, ApJ, 772, L1 [NASA ADS] [CrossRef] [Google Scholar]
  36. Matsumoto, J., & Masada, Y. 2019, MNRAS, 490, 4271 [NASA ADS] [CrossRef] [Google Scholar]
  37. Matsumoto, J., Komissarov, S. S., & Gourgouliatos, K. N. 2021, MNRAS, 503, 4918 [NASA ADS] [CrossRef] [Google Scholar]
  38. Mignone, A., & Bodo, G. 2005, MNRAS, 364, 126 [Google Scholar]
  39. Mignone, A., Plewa, T., & Bodo, G. 2005, ApJS, 160, 199 [CrossRef] [Google Scholar]
  40. Mignone, A., Bodo, G., Massaglia, S., et al. 2007, in JENAM-2007, “Our Non-Stable Universe”, 96 [Google Scholar]
  41. Mizuno, Y., Gómez, J. L., Nishikawa, K.-I., et al. 2015, ApJ, 809, 38 [Google Scholar]
  42. Mukherjee, D., Bodo, G., Mignone, A., Rossi, P., & Vaidya, B. 2020, MNRAS, 499, 681 [Google Scholar]
  43. Nandi, P., Stalin, C. S., Saikia, D. J., et al. 2023, ApJ, 959, 116 [NASA ADS] [CrossRef] [Google Scholar]
  44. Park, J., Kino, M., Nagai, H., et al. 2024, A&A, in press, https://doi.org/10.1051/0004-6361/202347562 [Google Scholar]
  45. Pereira-Santaella, M., Álvarez-Márquez, J., García-Bernete, I., et al. 2022, A&A, 665, L11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Perucho, M., & Martí, J. M. 2007, MNRAS, 382, 526 [NASA ADS] [CrossRef] [Google Scholar]
  47. Perucho, M., Martí, J. M., Laing, R. A., & Hardee, P. E. 2014, MNRAS, 441, 1488 [NASA ADS] [CrossRef] [Google Scholar]
  48. Rohoza, V., Lalakos, A., Paik, M., et al. 2023, ApJ, submitted, [arXiv:2311.00018] [Google Scholar]
  49. Rossi, P., Mignone, A., Bodo, G., Massaglia, S., & Ferrari, A. 2008, A&A, 488, 795 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Rossi, P., Bodo, G., Massaglia, S., & Capetti, A. 2020, A&A, 642, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Russell, H. R., Fabian, A. C., McNamara, B. R., & Broderick, A. E. 2015, MNRAS, 451, 588 [NASA ADS] [CrossRef] [Google Scholar]
  52. Sadler, E. M., Ekers, R. D., Mahony, E. K., Mauch, T., & Murphy, T. 2014, MNRAS, 438, 796 [Google Scholar]
  53. Torresi, E., Grandi, P., Capetti, A., Baldi, R. D., & Giovannini, G. 2018, MNRAS, 476, 5535 [NASA ADS] [CrossRef] [Google Scholar]
  54. Venturi, G., Cresci, G., Marconi, A., et al. 2021, A&A, 648, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Wilkes, J. I., Danehy, P., Nowak, R., & Alderfer, D. 2008, in 38th Fluid Dynamics Conference and Exhibit, AIAA, 4389 [Google Scholar]
  56. Willis, A. G., Strom, R. G., & Wilson, A. S. 1974, Nature, 250, 625 [Google Scholar]

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

∂ ∂ t ( ρ Γ ρ h Γ 2 v ρ h Γ 2 − p ρ Γ f ) + ∇ · ( ρ Γ v ρ h Γ 2 v v + p I ρ h Γ 2 v ρ Γ f v ) = ( 0 f g f g · v 0 ) , $$ \begin{aligned} \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), \end{aligned} $$(A.1)

where ρ,  Γ,  h,  v, and  p represent the rest-frame number density, the Lorentz factor, the proper specific enthalpy, the three-velocity in units of c, and the pressure of the fluid, respectively. The acceleration term fg is the specific external force three vector, set to fg = ∇pext(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)

h = 5 2 T + 9 4 T 2 + 1 , $$ \begin{aligned} h = \frac{5}{2}\mathcal{T} +\sqrt{\frac{9}{4}\mathcal{T} ^2+1}, \end{aligned} $$(A.2)

with 𝒯 = p/(ρc2), which approximates the Synge EoS of a single-species relativistic perfect fluid (Mignone et al. 2005). We adopted 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 was run in 2D cylindrical coordinates (r, z) in a domain [0, 20]×[0.5, 30], where the lengths are in units of z0, which represents the distance from the jet-launching site. The grid is uniform with 1000 × 3000 points in [0, 1.5]×[0.5, 20] and geometrically stretched with 400 × 700 grid points in the outer regions.

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

v z = 1 − 1 Γ j 2 z R , $$ \begin{aligned} v_z = \sqrt{1- \frac{1}{\Gamma _j^2}} \frac{z}{R}, \end{aligned} $$(A.3)

v r = 1 − 1 Γ j 2 r R , $$ \begin{aligned} v_r = \sqrt{1- \frac{1}{\Gamma _j^2}} \frac{r}{R}, \end{aligned} $$(A.4)

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

ρ j ( r , z , t = 0 ) = ρ j ( 0 , z 0 , 0 ) ( R R 0 ) − 2 , $$ \begin{aligned} \rho _j(r,z,t=0) = \rho _{j}(0,z_0,0) \left(\frac{R}{R_0}\right)^{-2}, \end{aligned} $$(A.5)

and pressure

p j ( r , z , t = 0 ) = p j ( 0 , z 0 , 0 ) ( R R 0 ) − 2 γ , $$ \begin{aligned} p_j(r,z,t=0) = p_{j}(0,z_0,0) \left(\frac{R}{R_0}\right)^{-2\gamma }, \end{aligned} $$(A.6)

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

ρ ext ( z , t = 0 ) = ρ ext ( z 0 , 0 ) ( z z 0 ) − η , $$ \begin{aligned} \rho _{ext}(z,t=0)=\rho _{ext}(z_0,0)\left(\frac{z}{z_0}\right)^{-\eta }, \end{aligned} $$(A.7)

p ext ( z , t = 0 ) = p ext ( z 0 , 0 ) ( z z 0 ) − η , $$ \begin{aligned} p_{ext}(z,t=0) = p_{ext}(z_0,0) \left(\frac{z}{z_0}\right)^{-\eta }, \end{aligned} $$(A.8)

where the power law index is η = 0.5. The exact value of η is not decisive, as long as it is smaller than two, so that the higher pressure of the external medium confines the relativistic jet.

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

ρ j ( 0 , z 0 , 0 ) ρ ext ( z 0 , 0 ) = 7.6 × 10 − 6 , p j ( 0 , z 0 , 0 ) p ext ( z 0 , 0 ) = 10 − 3 , $$ \begin{aligned} \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}, \end{aligned} $$(A.9)

while the external pressure is defined through the temperature: 𝒯ext(z0, 0) = 3 × 10−6. In our simulations, we set ρext, 0 as the density unit.

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

q = q ext + ( q j − q ext ) sech [ ( r z θ q ) α q ] $$ \begin{aligned} q = q_{ext} + \left(q_j-q_{ext}\right){{\,\mathrm{sech}\,}}\left[\left(\frac{r}{z\theta _q}\right)^{\alpha _q}\right] \end{aligned} $$(A.10)

in the inlet regions (Mukherjee et al. 2020). The index αq and the angle θq, which define the width and the radial scale of the smoothing, depend on the specific quantity q, and they are different for different values to avoid artificial local extrema in the energy and/or momentum (Abolmasov & Bromberg 2023). We set θΓ = 0.16 and αΓ = 8 for the Lorentz factor, and θρ, p = 0.29 and αρ, p = 10 for the density and the pressure.

We used outflow conditions at the right (r = 20) and upper (z = 30) boundaries, while we used reflective conditions at the axis r = 0. At the lower boundary z = zL = 0.5, for r > θjzL, we extended the initial analytical pressure and density profiles in the ghost cells to ensure dynamical equilibrium. We kept these values constant. The jet was injected in the region defined by 0.5 < z < 1 and 0 < r < θjz, where we set the velocities, density, and pressure defined in Eqs. A.3, A.4, A.5 and A.6. In addition, the fluxes of the Riemann solver were set to zero and hence the fluid variables have remained unchanged in this region. This was done to avoid spurious effects at the lower boundary. The 2D case was evolved until a steady state was reached, up to tf = 3000 in units of z0/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 × 550 × 1850 grid points, covering the physical domain [ − 5, 5]×[−5, 5]×[1, 30]. As in the 2D case, the grid is uniform only in the inner region, [ − 1, 1]×[−1, 1], × [1, 20], with 300 × 300 × 1500 cells, and stretched outside. We used outflow conditions at all boundaries, except at the lower one. In the 3D simulation, we set the lower boundary at z = z0 = 1. For x 2 + y 2 < θ j z 0 $ \sqrt{x^2+y^2} < \theta_j \, z_0 $, we have injection conditions with the jet parameters, and for x 2 + y 2 > θ j z 0 $ \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 ) = { 0 if v z ( x , y , z , t ) < 0.1 1 if v z ( x , y , z , t ) ≥ 0.1 $$ \begin{aligned} g(x,y,z,t) = {\left\{ \begin{array}{ll} 0\quad \text{ if} v_z(x,y,z,t) < 0.1\\ 1\quad \text{ if} v_z(x,y,z,t)\ge 0.1 \end{array}\right.} \end{aligned} $$(A.11)

which selects the jet section with a threshold on the propagation velocity vz. The average velocity ⟨vz⟩ is then the average of the distribution of the velocity on the x − y jet section, calculated as

⟨ v z ( z , t ) ⟩ = A − 1 ∫ xy v z ( x , y , z , t ) g ( x , y , z , t ) d x d y , $$ \begin{aligned} \langle v_z(z,t)\rangle = A^{-1}\int _{xy}v_z(x,y,z,t)g(x,y,z,t)dx dy\,, \end{aligned} $$(A.12)

where A = ∫xyg(x, y, z, t)dxdy.

All Figures

thumbnail 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 Γvz are opaque (gray), and which are transparent (black).
In the text
thumbnail 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 vz ≥ 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.
In the text
thumbnail Fig. 3.

Four-velocity Γvz in the x − y plane at different z* = 2.3,  4.7,  6.1,   and 13 for A, B, C, and D, respectively, defined in Fig. 2.
In the text
thumbnail Fig. 4.

Entertainment and deceleration as functions of time and altitude. Upper panel: mass flux, Φρ(z, t) = ∫xyρΓvzdxdy, starting from Φρ(z0) = 5.3 × 10−6, as a 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.
In the text
thumbnail Fig. 5.

Histogram displaying the fraction of cells with a given jet velocity Γvz 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).
In the text
thumbnail Fig. 6.

Slice in the y = 0 plane of the logarithm of the pressure in code units, log(p). The colder regions, with 𝒯 < 0.1, are not displayed.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.