How to Turn Jets into Cylinders near Supermassive Black Holes in 3D General Relativistic Magnetohydrodynamic Simulations

Accreting supermassive black holes ( SMBHs ) produce highly magnetized relativistic jets that tend to collimate gradually as they propagate outward. However, recent radio interferometric observations of the 3C 84 galaxy reveal a stunning, cylindrical jet already at several hundred SMBH gravitational radii, r  350 r g . We explore how such extreme collimation emerges via a suite of 3D general relativistic magnetohydrodynamic simulations. We consider an SMBH surrounded by a magnetized torus immersed in a constant-density ambient medium that starts at the edge of the SMBH sphere of in ﬂ uence, chosen to be much larger than the SMBH gravitational radius, r B = 10 3 r g . We ﬁ nd that radiatively inef ﬁ cient accretion ﬂ ows ( e.g., M87 ) produce winds that collimate the jets into parabolas near the black hole. After the disk winds stop collimating the jets at r  r B , they turn conical. Once outside r B , the jets run into the ambient medium and form back ﬂ ows that collimate the jets into cylinders some distance beyond r B . Interestingly, for radiatively ef ﬁ cient accretion, as in 3C 84, the radiative cooling saps the energy out of the disk winds; at early times, they cannot ef ﬁ ciently collimate the jets, which skip the initial parabolic collimation stage, start out conical near the SMBH, and turn into cylinders already at r ; 300 r g , as observed in 3C 84. Over time, the jet power remains approximately constant, whereas the mass accretion rate increases; the winds grow in strength and start to collimate the jets, which become quasi-parabolic near the base, and the transition point to a nearly cylindrical jet pro ﬁ le moves outward while remaining inside r B .


Introduction
Highly energetic gas outflows from accreting supermassive black holes (SMBHs) span orders of magnitude in distance; for instance, relativistic collimated outflows, or jets, in M87 reach distances of 10 kpc (Biretta et al. 1999), whereas the X-ray jet in OJ 287 extends out to 1 Mpc (Marscher & Jorstad 2011).These outflows may impact the formation of galaxies by injecting energy and momentum into the ambient gas through the process known as active galactic nucleus (AGN) feedback (see Fabian 2012;Parker et al. 2017).Understanding the AGN feedback is crucial for modeling galaxy evolution (Weinberger et al. 2018); jet activity can impact star formation (Ehlert et al. 2023), affect the cooling flows in galaxy clusters (Weinberger et al. 2023), and transport metal-enriched gas in galaxy cores (Kirkpatrick et al. 2009).Jets can produce high-energy (TeV) flares and accelerate hadronic cascades leading to multimessenger-electromagnetic and neutrino-emission, opening a novel window into the poorly understood dissipative processes around black holes (BHs; Mannheim 1995;Murase et al. 2012;Ansoldi et al. 2018;Keivani et al. 2018;Fang et al. 2023).
To meaningfully interpret these observations, it is important to understand how jets form and interact with the ambient medium.
Very long baseline interferometry observations have shed light on the properties of jet propagation, such as its acceleration and shape.Asada & Nakamura (2012) have found that the jet in M87 follows a parabolic collimation profile that transitions to conical at r ; 10 5 r g , where r g = GM/c 2 is the BH gravitational radius.This shape transition happens at the edge of the sphere of the gravitational influence of the SMBH, or the Bondi radius, defined as = r GM c s B 2 , where c s is the sound speed in the ambient medium.Asada & Nakamura (2012) and Levinson & Globus (2017) suggested that this transition from parabolic to conical in M87 might occur when the confining pressure profile changes.A study of a sample of 56 radio-loud AGNs shows the quasi-parabolic jet morphology inside the BH sphere of influence, which, due to the resolution constraints, may fit a parabolic-to-conical profile (Algaba et al. 2017).Boccardi et al. (2021) find that the jet in NGC 315 exhibits a parabolic shape and becomes conical well inside the Bondi radius; they discuss that the winds launched by a thick accretion disk might be responsible for this shape transition.High-resolution images of 3C 273 show that the jet exhibits a parabolic shape up to 10 7 r g before expanding conically well outside of the BH sphere of gravitational influence (Okino et al. 2022), which they expect to be approximately of the same order of magnitude as r B .The jet in Cygnus A features a parabolic shape up to r ; 2.4 × 10 4 r g ; 0.05r B , beyond which it turns into a cylinder; this suggests different environmental conditions than M87 (Boccardi et al. 2016).Although jets follow different collimation profiles at large radii, they typically have a parabolic shape near the BH.
Surprisingly, recent observations of 3C 84 have revealed an almost cylindrical jet collimation profile at 350r g  r  8000r g , where the jet maintains a near-constant cylindrical radius, R jet ≈ 250r g (Giovannini et al. 2018).They conclude that a parabolic expansion starting at the base of the jet at the BH event horizon cannot produce such a wide jet so close to the BH.Savolainen et al. (2023) note that this is a recently restarted jet, which is about 10 yr old, and deduce that the cylindrical jet shape needs a nearly flat pressure profile of the ambient medium, which can be due to a "leftover" cocoon, the shocked ambient medium by the previous jet.They report the emission around the restarted jet on (sub)parsec scales from the putative cocoon, which can be responsible for the collimation.The cylindrical jet collimation profile so close to the BH challenges the commonly used models that typically feature the parabolic expansion of the jets near the BH.
Significant progress has been made in describing the jet formation and propagation analytically and numerically (e.g., De Villiers et al. 2005;McKinney 2006;McKinney & Narayan 2007;Lyubarsky 2009;Tchekhovskoy et al. 2011;Chatterjee et al. 2019Chatterjee et al. , 2023; for reviews, see Blandford et al. 2019;Davis & Tchekhovskoy 2020).Bromberg et al. (2011) developed an analytical model of the relativistic unmagnetized jet propagating in the ambient medium.This model has been numerically verified in Harrison et al. (2018) and recently generalized analytically and calibrated numerically to jet propagation in an expanding medium (Gottlieb & Nakar 2022).This model creates a framework for understanding the jet collimation profile; they characterize the jet shape based on the jet injection opening angle θ 0 and the ratio of jet to rest-mass energy densities, L ˜, at the location of the jet head.
Yet, to understand the shape of the jet, it is crucial to account for the magnetic fields, since they can modify jet collimation (Komissarov et al. 2009) and give rise to 3D magnetic kink instabilities and dissipation (Lyubarskii 1999;Narayan et al. 2009;Bromberg & Tchekhovskoy 2016;Lalakos et al. 2023).Tchekhovskoy & Bromberg (2016) simulated relativistic magnetized jets launched by a spinning, perfectly conducting magnetized sphere.They immersed the sphere into a power-law density profile, mimicking the physical system outside the Bondi radius.However, their model did not include the central spinning BH that can launch jets magnetically via BH frame dragging (Blandford & Znajek 1977).They also did not account for the accretion disk winds, which collimate the jets.Chatterjee et al. (2019) carried out high-resolution 2D general relativistic magnetohydrodynamic (GRMHD) simulations with a large spatial separation (5 orders of magnitude) and measured the jet shape out to the Bondi radius scales.They found that the jet shape resembles the parabolic profile in M87 remarkably well throughout its entire length.Due to numerical constraints of 2D, Chatterjee et al. (2019) adopt a "standard and normal evolution" (SANE; Narayan 2012) state configuration, where large-scale magnetic fields are subdominant, and jets are less variable (Narayan 2012).In contrast, a system where magnetic flux accumulates on the BH to the point it becomes dynamically important and obstructs accretion is known as a magnetically arrested disk (MAD) state (Narayan et al. 2003;Tchekhovskoy et al. 2011;Chatterjee & Narayan 2022).Horizon-scale observations of M87 (Event Horizon Telescope Collaboration et al. 2019a) are consistent with both SANE and MAD magnetic fields, but the polarization measurements favor MAD due to the presence of organized poloidal magnetic fields (Event Horizon Telescope Collaboration et al. 2021).However, describing magnetic fields in AGN requires high-resolution polarimetric observations and inference from the numerical models, so the magnetic field remains a free parameter in the simulations.
A complete simulation of the full system in 3D, including jet launching by a rapidly spinning BH and propagation beyond the Bondi scales, is a numerically challenging task due to the large spatial and temporal separations (Lalakos et al. 2022(Lalakos et al. , 2023)).Ressler et al. (2021) simulated the formation and propagation of jets in the systems with a nonrotating constant ambient medium, r B = 100r g , and uniform magnetic fields tilted with respect to the BH spin axis.They found that the jets expand parabolically and then break apart due to the kink instability.Lalakos et al. (2022) studied jet propagation out to r B = 10 3 r g at the highest scale separation to date and found that the jets, which started out parabolically near the BH, collimated into cylinders near the Bondi radius, i.e., at large distances from the BH.What causes jets to turn cylindrical near the BH remains a mystery.The nature of the parabolic-toconical shape transition at the Bondi radius also remains unclear; it has never been seen in GRMHD simulations that included the ambient medium with a well-defined r B .
In this Letter, we investigate the dependence of the jet collimation profile on the ambient medium density and the thickness of the accretion disk via 3D GRMHD simulations.In Section 2, we describe our computational approach and simulation setup.In Section 3, we discuss how the jet shape depends on the properties of the accretion disk and ambient medium.We compare the results with the observations of 3C 84 and other AGN and discuss future work in Section 4. We conclude in Section 5. We adopt units such that G = M = c = 1, where M is the mass of the BH; units of time then become r g /c = 1.

Problem Setup and Computational Approach
We run the simulations using the GPU-accelerated GRMHD code H-AMR (Liska et al. 2022b).We use spherical coordinates r, θ, and f and a grid that is uniform in r log , θ, and f.We set the inner and outer boundaries of the computational grid at r in = 0.8r H and r out = 10 5 r g , respectively, where = + r a r 1 1 H 2 g ( ) is the event horizon radius and a is the dimensionless BH spin.Since we want to properly evolve the disk, which can be described as a nonrelativistic gas, we adopt an ideal gas law equation of state with the polytropic index of γ = 5/3.Note that relativistic outflows are typically magnetically dominated and have low gas pressure, so the choice of the exact value for the polytropic index does not affect the outflow evolution significantly.
All simulations have the resolution of the base grid, (N r , N θ , N f ) = (288,256,192), and use three levels of adaptive mesh refinement (AMR), which results in a maximum effective resolution of (2304,2048,1536) in the jets.We use the AMR refinement criterion that activates depending on the number of cells in the θ-direction in a jet or a cocoon (described in Gottlieb et al. 2022;Lalakos et al. 2023); we identify a jet based on the maximum Lorentz factor (as defined in Section 2.2) and cocoon based on a proxy for entropy (see Appendix B in Lalakos et al. 2023).We resolve the jet radius with at least five cells out to r = 10,000r g .

Initial Conditions
We start with a rapidly spinning BH of spin a = 0.9375 surrounded with an equilibrium hydrodynamic torus (Fishbone & Moncrief 1976).We place the inner edge of the torus at r disk,in = 20r g and the pressure maximum at = r r 41 disk,max g .These parameters result in a torus with the outer edge just inside r = 1000r g .We normalize the density such that r = max 1.To initialize the magnetic field in the disk, we set the toroidal component of the magnetic vector potential to r = - ) , which results in an initial poloidal magnetic field configuration.In this configuration, the magnetic flux surfaces follow the lines of constancy of density.We normalize the disk magnetic fields such that the ratio of the maximum thermal pressure to maximum magnetic pressure is = p p max max 100 gas mag .Beyond the outer edge of the torus, we place a constantdensity unmagnetized ambient medium; it starts at r = 1000r g and extends to the outer boundary of the grid.We set the internal energy of the gas such that the Bondi radius equals r B = 10 3 r g .We carry out a suite of three nonradiative simulations, M7, M6, and M5, which we name according to the ambient medium density; e.g., model M7 has the density ρ amb = 10 −7 in the code units.Since the luminosity, L, of 3C 84 is about 0.4% of its Eddington luminosity, L Edd (Plambeck et al. 2014), it is close enough to the standard 1% Eddington ratio above which cooling becomes significant.In fact, already at such low luminosities, the disk becomes radiatively efficient (Ryan et al. 2017;Liska et al. 2023aLiska et al. , 2023b)).To account for this, we also carry out another simulation, M5C, where we introduce a prescription to approximate radiative cooling; here, "C" in the name stands for cooling.We decrease the internal energy of the disk over the local Keplerian timescale to a target temperature profile chosen to correspond to a target disk scale height (Noble et al. 2009) that we choose to be H/R = 0.1.

Jet Definitions
In our analysis of the large-scale jet structure, we demarcate the jets using the parameter m r = -T u t r r ( ), where -T t r ( ) is the total energy flux in the radial direction, n m T is the stressenergy tensor, ρu r is the radial rest-mass flux density, and u r is the contravariant radial component of the proper velocity fourvector.The μ parameter can be thought of as the maximum possible Lorentz factor the jet could accelerate to if all of its energy flux were to convert into kinetic energy flux.To avoid division by 0 where u r vanishes (e.g., at the jet stagnation surface), we use an alternative form of μ given by the total specific energy of a fluid (Chatterjee et al. 2019;Lalakos et al. 2022): μ = −u t (h + σ + 1), where u t is the temporal component of the covariant four-velocity and a conserved quantity for a point mass; )is the specific (nonrelativistic) gas enthalpy, which excludes the rest-mass contribution; and magnetization is approximated as s r = p c 2 m 2 ( ) (see Lalakos et al. 2023).To isolate the jet regions, we examine the angular distribution of μ and set the cutoff value, μ jw , which demarcates the boundary between the BH-powered jets and disk-powered winds, such that at μ < μ jw , μ drops sharply to ∼1, which corresponds to nonrelativistic gas.We consider the regions at μ > μ jw , u r > 0 to be part of the jets and μ μ jw , u r > 0 to be part of the disk winds; we will specify μ jw values for each figure.In the above, we only consider regions with positive radial velocity, i.e., those with the flow moving away from the BH; this allows us to exclude the backflows and focus on the jets and winds.
Close to the BH, jets are highly magnetized, and the expression for the total specific energy, μ, is dominated by the contribution of the magnetization term, σ.This allows us to use a simplified jet definition to compute the jet power near the jet base, σ > 1, which focuses on the highly magnetized material.Conversely, we compute the disk winds power as the energy flux in the regions with σ < 1.

Jet Collimation Dependence on the Ambient Medium and Disk Cooling
Figure 1 shows 2D density snapshots in the xz-plane in our simulations at t = 20k; each column corresponds to a different simulation.The top row shows the large-scale structure of the jets and their cocoons, while the bottom row shows the respective zoom-in on the jet base; white dashed contours show jet boundaries.We define the jet as μ 5 for model M6 and μ 4 for models M7, M5, and M5C. Figure 1(a) shows that in our lowest ambient density model M7, by this time, the jet, which is launched by the BH, collides with and drills a hole in the constant-density ambient medium.This interaction redirects the highly magnetized jet material at the jet head into backflows.These low-density (dark blue) regions are separated from the jet by the denser (light blue) disk winds.As the backflows return to the jet base, they collimate the jets into cylinders at r  4000r g .As Figure 1(b) shows, the jets transition to cylindrical geometry only well outside of the Bondi radius; they expand laterally closer to the BH, as we discuss below.Figures 1(c) and (e) show that models with higher ambient medium density, M6 and M5, exhibit more prominent backflows.The jets propagate slower as they struggle to drill through a denser ambient medium and become narrower due to a stronger collimation by a more energetic cocoon.In contrast, at distances within the Bondi sphere, which coincides with the characteristic size of the disk, the jets exhibit a shape similar to model M7, as we see in Figures 1(d) and (f).This suggests that the disk winds determine the collimation profile of the jets near the BH, whereas the ambient medium impacts the jet farther away.Figure 1(f) shows that even for the strongest backflows present in model M5, the disk and its winds prevent the backflows from reaching the jet base; the jet appears to turn cylindrical only after the disk ends.Thus, the fact that jet shape differences in models M7, M6, and M5 emerge only at large distances outside the torus indicates that the disk is a key factor in setting the collimation near the BH.
To understand how disk cooling impacts jet evolution, we consider model M5C with disk cooling.Figure 1(g) shows that the jet in model M5C struggles to drill through the ambient medium, develops nonaxisymmetric features, and creates even more prominent backflows than in model M5, suggesting that the jet becomes weaker and more prone to instabilities. Figure 1(h) shows that the cooling causes the disk to collapse toward the midplane; this allows the backflows to reach smaller distances.The jet becomes wider, expands laterally at smaller distances, and develops a cylindrical shape closer to the BH, well inside the Bondi radius.Movies of M5 and M5C time evolution are available here. 9ince the jet behavior qualitatively depends on the strength of the jet with respect to the ambient medium, we quantify how strong the jet is. Figure 2 shows the BH mass accretion rate measured at r = 5r g ,  M; the accumulated absolute magnetic flux on the BH calculated at the BH event horizon, Φ BH ; and the power of the jet, P jet , at r = 5r g .We calculate the absolute magnetic flux on the BH as Φ BH = 0.5Φ|B r |dA θf , where the integral is over the entire sphere with r = r H and the factor of 0.5 converts it to one hemisphere.To compute the power of the jet, we compute the energy outflow rate, , in the highly magnetized regions, where the magnetization exceeds unity, s r , where p m is the magnetic pressure and ρ is the fluid frame mass density.All jet quantities are computed only for the top jet, so we constrain θ, 0 θ π/2; the bottom jet shows qualitatively similar behavior.
Figure 2(a) shows that models M7, M6, and M5 reach a steady state by t ; 15k when the mass accretion rate reaches a quasi-steady state.With the mass accretion rate staying approximately constant and the magnetic flux decreasing, as seen in Figure 2(b), the power of the jets in Figure 2(c) decreases as expected (Tchekhovskoy et al. 2011).Since models M7, M6, and M5 exhibit very similar behavior in  M, Φ BH , and P jet , we conclude that the differences in the ambient medium density do not affect the jet dynamics, which instead is entirely determined by the accretion disk.Over time, as we see in Figure 2(c), the jet in our highest ambient density model M5 becomes progressively weaker until it disappears completely at t ; 60k.The Bondi radius is shown with a dotted circle; note that each simulation has a different spatial scale.Panels (a), (c), and (e) show that as ρ amb increases, the jet propagates more slowly and becomes visually less stable.In panel (a), where ρ amb is lower, the cocoon is less prominent because the jet is stronger with respect to the ambient medium, and the jet creates a weaker shock.The jet interaction with the ambient medium redirects the material from the jet head and creates backflows, the low-density regions separated from the jet by denser disk winds as seen in panels (a), (c), (e), and (g).The backflows are the most prominent in the highest-density simulations (e) and (g), since it is harder for the jet to drill through the ambient medium.The disk cooling causes the jet to become weaker and more prone to instabilities, and panel (g) shows that the jet appears even less stable.The bottom row is a zoom-in between z = 0r g and 4000r g ; the white dashed lines show jet boundaries defined as μ 5 for model M6 and μ 4 for models M7, M5, and M5C (see Section 2.2).Despite the difference in the large-scale jet morphology between models M7, M6, and M5 at r ; 3000r g , panels (b), (d), and (f) show a similar jet morphology at the base, which is determined by the torus and its winds.With the cooling present, which impacts the torus, the jet in panel (h) opens up at smaller distances and appears wider.We notice that the backflows reach closer to the BH, since the disk scale height decreases due to cooling.
Unlike the changes in the ambient medium density, disk cooling significantly impacts the near-BH behavior of the system.Figure 2(a) shows that the mass accretion rate in model M5C flattens out at t ; 15k-30k at a level similar to models M7, M6, and M5; afterward, it starts increasing until it saturates at t ; 50k.We will use the similarity in mass accretion rate across all models and compare the jet collimation profiles at t = 20k (Section 3.2).
As we will see in Section 3.3, at t  30k, more of the disk in model M5C loses thermal support, collapses toward the midplane, and moves closer to the BH.Although this results in a significant increase in the mass accretion rate, Figure 2(c) shows that the power of the jet approximately flattens out at t  15k.This leads to jet efficiency, ), decreasing in time.This allows us to study the jet shape evolution for a range of η jet values (Section 3.3).

Jet Radius Comparison
Figure 3 shows the jet radius, R jet , as a function of distance from the jet base, z.We calculate the jet radius as where θ jet is the jet opening angle, calculated from the area subtended by the jet: Figure 3 depicts the jet radius as a function of distance from the BH.For model M7, the jet expands parabolically up to r ; 200r g ; this profile is characteristic of jets launched via the Blandford-Znajek mechanism due to the collimation by the disk winds (Blandford & Znajek 1977; see also McKinney & Narayan 2007;Chatterjee et al. 2019).When the jet outruns the disk winds, they cannot keep collimating the jet further, and it expands laterally and transitions into a conical shape.However, such a free, ballistic expansion cannot last forever in the presence of the ambient medium.Eventually, at r ; 5000r g , due to the constant-pressure environment formed by the backflows, the jet gets collimated into a cylinder, as predicted in the asymptotic solution by Lyubarsky (2009).Here, we demonstrate for the first time that such a shape transition naturally emerges in jets selfconsistently launched by a rotating accreting BH.
Similar to model M7, the jets in models M6 and M5 exhibit a parabolic shape near the jet base, inside the Bondi radius.Figure 3 shows that their shape also transitions to the conical shape at r ; 200r g .However, Figure 3 also reveals that the jets in models M6 and M5 turn into cylinders closer to the BH, at r ; 3000r g and 1500r g , respectively; a denser ambient medium results in stronger backflows that can collimate the jets closer to the BH.In addition to collimating jets into cylinders closer to the BH, models M6 and M5 feature narrower cylindrical jets because the jets become less powerful relative to the denser ambient medium.
Interestingly, Figure 3 shows that the jet in model M5C expands right from the outset and becomes significantly wider than that in model M5 on scales of a few hundred r g , which is the characteristic size of the torus.We associate such rapid jet expansion at small radii with the disk cooling that weakens the disk outflows and suppresses the jet collimation.However, already at r ; 300r g , the jet in model M5C becomes cylindrical due to the collimation by the backflows.In fact, in model M5C, the backflows manage to get closer to the jet base and turn the jet into a cylinder on a smaller scale than in other models due to the cooling-reduced scale height of the disk and weakened disk outflows.The transition into a cylindrical shape happens well inside the Bondi radius and may relate to the observations of the cylindrical jet in 3C 84 at similar distances, r  350r g (Giovannini et al. 2018)

Time Evolution of the Radiatively Efficient Model, M5C
Because the disk winds are crucial for jet collimation inside the Bondi radius, and disk cooling impacts the disk structure, we study how the cooling impacts the mass flow rates and the power of the outflows.Figure 4(a) shows the magnitudes of mass inflow and outflow rates,  M in and  M out .We define inflow and outflow based on the sign of radial velocity, u r ; inflow has negative and outflow has positive radial velocities.We adopt the sign convention that the mass accretion rate is positive if directed into the BH and negative if directed away from the BH.Based on this definition, the net mass accretion rate equals Figure 4(a) shows that at t  30k, the mass inflow follows the mass outflow, resulting in a relatively low net mass accretion rate, as seen in Figure 2(a) at t  30k.At t ; 30k, the mass inflow rate starts increasing relative to the mass outflow rate.We can interpret this increase in  M in as follows.As the disk loses thermal pressure support due to radiative cooling and collapses, more mass comes closer to the BH, and the mass accretion timescale decreases;  M in and  M both increase, as seen in Figures 2(a) and 4(a).In Figure 4(b), we compute the jet power by considering the energy flux in the regions with the magnetization σ > 1 (see Section 2.2); the rest of the energy flux we attribute to the winds.At t ; 30k, Figure 4(b) shows that the disk winds become more powerful and start dominating the jets.As seen in Figure 4(c), the jet energy efficiency  h = P Mc 2 ( ) decreases as expected due to a nearly flat jet power and an increasing mass accretion rate.
To understand how the jet shape changes as a function of jet and wind energy efficiencies, in Figure 5, we analyze the time series of model M5C snapshots, taken from t = 20k to 50k, after which the jet disappears, suppressed by the winds and the ambient medium.Figure 5(a) shows that the outer edge of the torus reaches r ; 500r g , and only the inner regions of the torus are significantly impacted by the cooling.Because the inner disk regions collapse toward the midplane and the winds weaken due to radiative cooling, the jet expands nearly conically near the BH, as seen in Figure 3.At later times, cooling starts to impact the outer parts of the disk as well, and the sequence of panels (b), (c), and (d) in Figure 5 shows how the disk reduces in size over time.
Figure 5(e) shows jet collimation profiles time-averaged over the time interval of Δt = 2k centered around the times shown (see legend); the shaded regions show one standard deviation.In Figure 5(e), while computing the collimation profile of M5C, we use μ 3 for the jet definition, since it becomes progressively weaker and, therefore, its values of μ decrease.Figure 4 shows that the winds become stronger than the jet at t  30k. Figure 5(e) shows that as a result, the winds more strongly collimate the jets at 30k  t  50k at r  500r g than at t = 20k.Note that the rise of the wind power at t ; 30k pushes the backflows away from the jet base and causes the jet to become quasi-parabolic near the BH and M values are similar and the P jet values are within a factor of 2 between model M5C and the rest of the models (M7-M5).Jets in models M7, M6, and M5 have a similar parabolic collimation profile extending out to r ; 200r g , beyond which they expand conically.Outside of the sphere of the BH gravitational influence, r  r B , the jets in models M7, M6, and M5 become cylindrical; as we increase the ambient medium density, the transition point moves inward.The jet in model M5C initially expands conically up to around 300r g and then becomes nearly cylindrical.
Figure 4.The disk cooling in model M5C initially allows the jets to energetically dominate the disk wind, but over time, the mass accretion rate increases, and the winds become stronger, whereas the jet power slowly decreases because it is limited by the amount of available large-scale vertical magnetic flux in the system.All quantities are measured at r = 15r g to capture the disk winds and the jet but exclude backflows.Panel (a) shows the inflow, outflow, and net mass accretion rates.Until t ; 30k, the mass inflow and outflow rates are comparable.At t  30k, due to cooling, the disk material comes closer to the BH, and its accretion timescale decreases.More material starts falling in while approximately the same amount of material flies out, which results in an increased inward net mass flux.In panel (b), at t ; 30k, the power of the wind becomes larger than the power of the jet, and the winds become more prominent.The power of the jet reaches an approximately constant value at t ; 15k and slowly decreases over time.Panel (c) shows the jet and wind efficiencies η j and η w .The efficiency of the jet drops drastically due to the increasing mass accretion rate, which allows us to analyze the jet shape at different jet efficiencies.
turn into a cylinder farther out.The jet profile at later times shows increasingly more variability since it becomes less stable.Despite the stronger collimation near the BH, the jet still ends up turning cylindrical inside r B , albeit at larger distances, r ; 400r g -700r g .

Comparison with Observations
Because the density unit in our simulation is free, to map our scale-invariant simulations to real systems, we calculate the jet stability parameter in physical and code units as derived in Bromberg & Tchekhovskoy (2016) and Tchekhovskoy & Bromberg (2016), where P jet is the jet power, n is the ambient medium number density, m p is the proton mass, r is the distance to the jet head, γ jet is the Lorentz factor of the jet fluid, and K is a constant factor.Qualitatively, Equation (3) involves the power of the jet in the numerator and the density of the ambient medium times the length scale squared in the denominator.We consider the length scale equal to the Bondi radius since it is the order of magnitude of the jet length in the system.Therefore, we can compute a dimensionless ratio of the quantities that determines the jet stability: By computing the simulated and observed values of λ, we can make quantitative comparisons of jet power normalized by their ambient density between simulations and observations.In the code units, show that at t = 30k and 40k, the jet becomes quasi-conical near the BH because the jet power remains approximately constant, but the wind power increases in time.Panel (d) shows that t = 50k, the jet assumes a quasi-parabolic shape.Interestingly, throughout the simulation, the jet profile remains robustly quasi-cylindrical at r  700r g , i.e., well inside the Bondi radius.At even later times, t  60k, the jet becomes unstable (not shown).
For 3C 84, we use the mass of the BH, M = 8 × 10 8 M e (Scharwächter et al. 2013).From model FL from Fujita et al. (2016b), we adopt the ambient medium number density, n = 0.05 cm −3 ; the Bondi radius, r B = 8.6 pc; and the power of one jet, P jet = 2.8 × 10 43 erg s −1 and compute the stability parameter, λ 3C 84 ∼ 0.02, using Equation (4).Because for 3C 84, L/L Edd ; 0.4% (Plambeck et al. 2014), which is sufficiently close to the canonical 1% divide between radiatively inefficient and efficient flows, we expect the accretion flow in 3C 84 to be to some degree radiatively efficient.Hence, our model M5C provides the most relevant comparison.Indeed, λ 3C 84 is very close to that in model M5C, λ M5C = 0.02 (the stability parameters for jets in M6 and M7 are 10 and 100 times higher, respectively).Thus, our radiatively efficient model M5C offers insight into what can turn 3C 84 jets into cylinders so close to the BH.If the disk outflows in 3C 84 are weakened by the radiative cooling of the accretion disk, as in our model M5C, they will provide weak collimation to the jet and let it expand ballistically, i.e., essentially radially (Figure 5).The weakened and weakly collimated disk winds are unable to put up a fight against the powerful backflows either; produced by the jet interaction with the ambient medium, the backflows push the winds out of the way and reach deep inside the Bondi radius, where they collimate the jets into cylinders.The significant impact of the backflows on the jet shape deep inside r B in our model M5C suggests that they may be important to the jet collimation in real systems well inside the Bondi radius.In future work, we will compute synthetic images of jets and backflows to determine their observational signatures and compare them to the images of the 3C 84 jet and cocoon.It will also be interesting to quantify the mixing of backflows with winds, which can potentially impact the observational signatures of the system.
Interestingly, whereas our jets collimate into cylinders at about the same distance as in 3C  (Giovannini et al. 2018).Savolainen et al. (2023) note that the jet in 3C 84 has recently restarted and is only 10 yr old.This duration translates into t = 76k t g , which is comparable to the duration of our model M5C and is about 4 times as long as t = 20k, at which we measured the jet width (Section 3.2).If the mass accretion rate in our model M5C did not increase in time, we would expect that by t = 76k, our simulated jets would gradually expand laterally, and their width would more closely resemble that observed in 3C 84.We can further constrain our choice of parameters by comparing the ratio of the thermal energy of the cocoon with the energy delivered by the jet in 3C 84, denoted E c and E, respectively.Savolainen et al. (2023) . To compute the total thermal energy of the cocoon, E c sim , in our M5C model, we use the entropy proxy as a criterion that separates the cocoon from the ambient medium, r = g S u g ˜(see Appendix B in Lalakos et al. 2023); we choose the entropy proxy cutoff value of > S 1.85 ˜.We also exclude jets based on the μ cutoff but include backflows and disk winds.By integrating the jet power over time, we compute the cumulative energy delivered by the jet , comparable to the observed value of . We note that the jet injection hydrodynamic simulations presented in Savolainen et al. (2023) also observe the jets turning into cylinders due to the collimation by the cocoon.However, their simulations cannot explore the jet shape close to the BH because they do not include the BH and prescribe the jet injection at a large distance from the BH.For M87, we adopt the mass of the BH, M = 6.5 × 10 9 M e (Event Horizon Telescope Collaboration et al. 2019b); the number density at r = 1 kpc, n ; 0.2 cm −3 ; the Bondi radius, r B = (0.11 ± 0.02) kpc (Russell et al. 2015);10 and the jet power, P jet ∼ 10 44 erg s −1 (Bicknell & Begelman 1996;Owen et al. 2000).Plugging these values into Equation (4) gives λ M87 = 10 −4 .The accretion flow in M87 is highly sub-Eddington, L/L Edd ∼ 10 −6 (e.g., Broderick & Tchekhovskoy 2015), and thus has a low radiative efficiency.Thus, we perform the comparison against our nonradiative models, M7-M5.We find that the jet stability parameter in model M5 at t ; 20k is λ ; 0.05, i.e., a factor of ;500 larger than that in the M87 jet.Thus, our jets are much more powerful than in M87.Despite this, we can try to draw conclusions from our models for the case of M87.In our models M7, M6, and M5, the jets consistently develop a shape transition from parabolic to conical before becoming cylindrical around the Bondi radius.This suggests that in M87, the disk winds would collimate the jets out to a fraction of the Bondi radius, after which the jets would first turn conical and, shortly outside r B , cylindrical.However, the transition to cylindrical is not observed in M87.It is possible that for weaker jets, 3D magnetic kink instabilities can modify the picture found in this work; at smaller jet powers, the jets will be less stable, are more likely to break apart due to the magnetic kink instability, and are unable to propagate beyond a critical distance (Tchekhovskoy & Bromberg 2016).Once the jets become stuck, their exhaust outruns the jet head; as a result, the backflows do not form and cannot collimate the jets into cylinders.Whether this leads to a conical shape beyond the Bondi radius, as observed in M87, will be the topic of future work.Note that because P jet in model M5 decreases with increasing time such that by t = 60k the stability parameter drops down to λ M5,60k  10 −3 , late-time evolution of model M5 might provide a better quantitative comparison to M87.We leave this to future work.
For NGC 315, the inferred BH mass ranges from M = 3 × 10 8 M e (Gu et al. 2007) to M = 3.4 × 10 9 M e (Bettoni et al. 2003).Following Boccardi et al. (2021), we adopt an intermediate BH mass, M = 1.3 × 10 9 M e (Satyapal et al. 2005).We take the Bondi radius to be r B = (51 ± 25) pc (Inayoshi et al. 2020;Boccardi et al. 2021), the ambient medium number density to be n = 2.8 × 10 −1 cm −3 (Worrall et al. 2007;Ricci et al. 2022), and the jet power to be P jet ; 1.4 × 10 44 erg s −1 (Ricci et al. 2022).This yields a stability parameter, λ NGC 315 ∼ 2 × 10 −4 -2 × 10 −3 , 1-2 orders of magnitude below our models M5 and M5C, as evaluated at t = 20k.Because the accretion flow in NGC 315 is in the lowluminosity, radiatively inefficient regime, L/L Edd ; 5 × 10 −4 = 1 (Gu et al. 2007), our energy-conserving model M5 provides the most appropriate comparison.The jet in our model M5 transitions from a quasi-parabolic to a quasi-conical shape at a factor of a few smaller radius than r B .This is very different than what is observed; observations indicate that the parabolicto-conical jet transition in NGC 315 happens at subparsec scales, at   ŕ r 5 10 0.6 tr 3 g pc, which is much smaller than the Bondi scale, r B ; 4 × 10 5 r g ; 50 pc (Boccardi et al. 2021).In our models, the near-BH jet collimation is largely controlled by the disk winds, with the jets transitioning from parabolic to conical on the length scale comparable to the radial extent of our thick accretion torus.Thus, the observed parabolic-toconical transition distance can allow us to infer the radial extent of the inner thick radiatively inefficient accretion flow (RIAF) in NGC 315:    ŕ r r r 5 10 0.01 RIAF tr 3 g B .Note that similar to M87, at late times, the stability parameter of our M5 model drops down to λ M5,60k  10 −3 .Hence, the late-time evolution of our M5 model might provide a quantitative comparison to NGC 315.We leave such a detailed comparison to future work.
For Cygnus A, we use the BH mass of M = 2.5 × 109 M e (McNamara et al. 2011), the ambient medium number density of n = 3 cm −3 (Fujita et al. 2016a), the Bondi radius of r B = (0.015-0.12) kpc (Fujita et al. 2016a;Nakahara et al. 2019; see footnote 10), and the jet power of P jet = 3.9 × 10 45 erg s −1 (McNamara et al. 2011).This results in λ CygA = 2 × 10 −4 − 0.015.Since Cygnus A has L/L Edd ∼ 1% (Reynolds et al. 2015), disk cooling might be significant; we will compare it with both cooled and energy-conserving runs.At t = 20k, we have λ M5,20k = 0.05 for model M5 and λ M5C = 0.02 for model M5C.We conclude that λ in the jets in models M5 and M5C is at the upper range of the inferred λ range for Cygnus A. At late times, t ; 50k, the jets in model M5C exhibit a similar jet collimation profile to that observed in Cygnus A: the jets expand quasi-parabolically at small radii and switch to cylindrical inside the Bondi radius (Boccardi et al. 2016).Interestingly, recent observations indicate that the Cygnus A jet radius seemingly discontinuously increases by an order of magnitude across r = r B and thereafter continues expanding parabolically out to kiloparsec scales (Nakahara et al. 2019).The nature of the jet radius jump and kiloparsec-scale parabolic expansion is presently unclear.The last jet shape snapshot, shown in purple in Figure 5, reveals signs of jet recollimation at r ∼ 4r B , i.e., pinching followed by bounce; this feature might be related to the observed discontinuous increase in Cygnus A jet radius at r ∼ r B but does not explain the resumption of parabolic jet shape at r ?r B .Longer-term simulations will be required to investigate the larger-scale jet collimation properties.

Limitations and Future Work
Whereas our simulations have a high resolution that is enough to resolve the jet radius out to r = 10 4 r g with five cells across the jet radius, we have simplified the setup to make it less numerically expensive.We start with a torus in hydrostatic equilibrium, which later launches winds that are dynamically important to the jet collimation.However, as Penna et al. (2013) point out, the Fishbone & Moncrief (1976) prescription of the torus sets the Bernoulli parameter, which determines whether the material is bound and impacts the amount of the disk material that leaves the system via winds.Penna et al. (2013) remark that the final Bernoulli parameter appears to keep the memory of the initial conditions and depends on the initial Bernoulli parameter value prescribed by the torus solution.If a more realistic torus has more tightly bound material, then we expect to see weaker winds and even wider jets inside the Bondi radius.Note that the cooling function decreases the internal energy of the gas and makes it more tightly bound, thereby helping to erase the memory of the initial torus parameters.To investigate how the initial prescription of a torus impacts the power dynamics between disk winds and backflows and, overall, jet evolution, studies of the jet collimation profile for different initial and selfconsistently formed tori are needed at the largest scale separations (e.g., Ressler et al. 2021;Lalakos et al. 2022Lalakos et al. , 2023)).
In addition to the accretion disk structure, we assumed a magnetic field configuration, which results in the SANE disk and jets.MAD jets are also possible and may have different collimation profiles; they tend to be more powerful and show more variability, so for the same ambient density and jet power, we expect less collimated jets, less prominent backflows, and stronger time dependence (Igumenshchev et al. 2003;Narayan et al. 2003;Tchekhovskoy et al. 2011;Chatterjee & Narayan 2022).
For our torus in M5C, we use a predefined cooling function (Noble et al. 2009) to cool the disk toward the target scale height, H/R = 0.1.The cooling function reduces the internal energy of the disk but does not include radiation-driven outflows (Huang et al. 2023), which may strengthen the disk winds and more strongly collimate the jets.In our study, the inner parts of the disk launch the winds and, therefore, contribute significantly to the jet collimation, and the disk in M5C reaches the target scale height by t = 20k out to r ; 50r g .Because the timescale of our cooling function equals the Keplerian timescale at each radius, the cooling introduces time dependence into the accretion disk evolution.For instance, by t ; 2k, the cooling function has already cooled the disk past its pressure maximum, thereby robbing the disk of thermal pressure support and "pulling the rug" from predominantly pressure-supported outer regions of the initial torus.This causes the outer regions of the torus to fall toward the BH.Free-falling outer layers of the torus, located at r = 10 3 r g , would reach the BH at ~t r r r c 30k ff g 3 2 g ( ) , where for simplicity we neglected factors of order unity.This is comparable to the time at which the mass accretion rate starts to increase (Figure 4).
In nature, we expect that radiatively efficient accretion flows allow the jets to expand laterally without the associated increase in the mass accretion rate.This is because the coolinginduced collapse of our initial torus is the result of our torus being predominantly thermal-pressure-supported at large radii.In hindsight, if we chose a different, more natural initial condition-that of an accretion disk in which centrifugal and thermal pressure forces are comparable to each other at all radii -the cooling would shrink the disk by at most a factor of a few and would not lead to the dramatic increase in the accretion rate we see in our simulation M5C.Although the increase in mass accretion rate and disk wind power would probably not occur in nature, it allows the ratio of jet to disk power in our simulation to sample a range of values (Figure 4) and enables us to study the dependence of jet shape on this ratio.The global structure of the disk can impact the supply of the magnetic fields and, thereby, the jet power.As the disk becomes thinner, it could lose its ability to advect large-scale magnetic fields inward as theorized by Lubow et al. (1994), which would make the jets weaker.To fully understand the impact of disk cooling on the shape of the jets, 3D radiation-transport GRMHD simulations of sub-Eddington accretion disks are needed (e.g., Liska et al. 2022aLiska et al. , 2023aLiska et al. , 2023b)).
In this work, we explored a single value for the constant density of the ambient medium in the presence of cooling in model M5C.Whereas this model allows us to interpret the observations of 3C 84, the value of the number density at the Bondi radius, n, which we need to compare λ obs and l sim , is poorly constrained.We have taken the value of n from model FL in Fujita et al. (2016b), which adopts the observed value of n = 0.05 cm −3 at kiloparsec scales (Fabian et al. 2006) and assumes a flat density profile down to the Bondi radius.In particular, Fujita et al. (2016b) find that their model FL is preferable to their model EX, which assumes a power-law density profile and results in a higher density at the Bondi radius than model FL.More studies for a range of ambient density values are needed to better understand the jet shape in systems with radiatively efficient accretion flows.We expect that denser ambient media will result in more prominent backflows that may be able to turn jets into cylinders at even smaller distances from the BH.
To be able numerically to follow the jet shape well outside the Bondi radius, we reduce the scale separation and choose r B = 10 3 r g .However, the Bondi radius of an SMBH is of the order of 10 5 r g -10 6 r g (e.g., Russell et al. 2015).This assumption implies that the ambient medium in our simulations has a higher internal energy, making it harder for jets to propagate.The opening angle of the jets at r B is set by the collimation of the jets by disk winds.For more realistic values of the Bondi radius, the winds will have a larger range of distances over which to collimate the jets; the jets will become more collimated and, potentially, faster.If numerically possible, simulations with larger values of the Bondi radius will shed light on the scale separation dependence of the problem.
In our simulations, we adopt constant ambient density and pressure profiles to evaluate the impact of the ambient medium on the jet collimation at zeroth order.However, the radial ambient density profiles in AGN typically follow a power law; for instance, in M87, the density around the Bondi radius scales as ρ ∝ r −1 (Russell et al. 2015).Such a density profile is effectively equivalent to a less dense ambient medium, resulting in jets turning into cylinders at larger r; we expect the jets to propagate faster and be less collimated.Depending on the radial power-law index of the ambient medium distribution, the backflows might not form at all, since, at the jet head, the ambient medium can be negligible compared to the jet and unable to redirect the jet material and might not be able to turn the jets into cylinders.Moreover, we expect that in real systems, the density profile is highly inhomogeneous due to the previous jet activity and external factors (Fabian 2012).As the ambient medium plays a key role in the jet collimation, in future work, we will explore how different ambient medium density profiles impact the shape of the jet.
In Section 3.2, we measure the jet shapes in models M7, M6, and M5 at t ; 20k.The jet shapes show little variability at this time.This makes our measurements robust.However, we expect that, as time progresses, the radii of transitions between the parabolic and conical regions, as well as between the conical and cylindrical regions, will move outward as Figure 5 shows for model M5C.However, in model M5, the jet becomes energetically dominated by the disk winds, which causes the jets to be stifled and thereby set a characteristic timescale of jet evolution.In future work, we will study the time dependence of the jet shape break radius to understand how it maps onto the observations of jets whose age is known (e.g., the 3C 84 jet).

Summary
In this Letter, we present 3D GRMHD simulations of BHlaunched jets interacting with disk winds and the constantdensity ambient medium at the largest Bondi-to-BH scale separation to date, r B = 10 3 r g .We find that when we increase the ambient medium density by 2 orders of magnitude, the jets exhibit parabolic shapes and become conical at r ; 200r g , independent of the ambient medium density; this is because the radial extent of the accretion disk determines the near-BH collimation profile of the jets.Formed due to the jet interaction with the ambient medium, backflows create a constant-pressure jet environment at larger distances, r  r B , e.g., r ; 6r B for model M7, and turn the jets into cylinders.The conical-tocylindrical transition point moves inward for higher ambient medium densities or less powerful jets relative to the environment; they struggle to drill through the ambient medium and form more prominent backflows that collimate the jets into cylinders at smaller radii.For instance, when the jet-to-ambient energy ratio is lower by factors of 10 and 100, the jets turn into cylinders at r ; 3r B and r ; r B , respectively.
We considered AGN, such as M87 and NGC 315, that have RIAFs.We found that although our total energy-conserving models tended to feature more powerful jets than in these systems, they allowed us to interpret the observations qualitatively.We speculate that the jet in M87 might fall apart due to the kink instability before forming backflows; therefore, it does not have the cylindrical collimation we have seen in models M7, M6, and M5.Our simulation results suggest that the parabolic-to-conical jet shape transition occurs where the thick disk ends.Thus, the surprisingly small distance of the jet shape transition in NGC 315 (much smaller than r B ) might be related to the size of the thick accretion disk in this system.
However, none of these total energy-conserving (nonradiative) models result in cylindrical jets well inside r B , unlike what is seen in Cygnus A and 3C 84.Because Cygnus A and 3C 84 accretion flows are expected to be radiatively efficient, we explored if disk cooling can affect the jet collimation profile by prescribing a cooling function in our model M5C.Radiative cooling saps energy from the disk winds and weakens them.This enables the jet in our model M5C to energetically dominate the disk winds, conically expand in a ballistic fashion, and become exceptionally wide close to the BH.The weaker disk outflows allow jet backflows, formed by the jetambient-medium interaction, to reach deep inside the Bondi radius and collimate the jets into cylinders already at r ; 300r g , similar to the 3C 84 jet.Over time, the disk winds become more powerful and cause the jet to become quasi-parabolic near the base and cylindrical farther out but still inside the Bondi radius, similar to Cygnus A. Summing up, whereas our model M5C is simplified, it offers a controlled experiment that reveals a qualitatively new, long-sought behavior: conical jet expansion near the BH and backflows collimating jets into cylinders well inside r B .

Figure 1 .
Figure1.A constant ambient medium collimates the jets, which become cylindrical, with a higher density resulting in stronger collimation and narrower jets; surprisingly, the disk cooling allows the jets to become cylindrical closer to the BH.Movies of M5 and M5C are available here (see footnote 9).Panels show 2D density slices in the xz-plane at t = 20k with each column corresponding to a different simulation.The top row is a zoom-out to show the full structure of the cocoon.The Bondi radius is shown with a dotted circle; note that each simulation has a different spatial scale.Panels (a), (c), and (e) show that as ρ amb increases, the jet propagates more slowly and becomes visually less stable.In panel (a), where ρ amb is lower, the cocoon is less prominent because the jet is stronger with respect to the ambient medium, and the jet creates a weaker shock.The jet interaction with the ambient medium redirects the material from the jet head and creates backflows, the low-density regions separated from the jet by denser disk winds as seen in panels (a), (c), (e), and (g).The backflows are the most prominent in the highest-density simulations (e) and (g), since it is harder for the jet to drill through the ambient medium.The disk cooling causes the jet to become weaker and more prone to instabilities, and panel (g) shows that the jet appears even less stable.The bottom row is a zoom-in between z = 0r g and 4000r g ; the white dashed lines show jet boundaries defined as μ 5 for model M6 and μ 4 for models M7, M5, and M5C (see Section 2.2).Despite the difference in the large-scale jet morphology between models M7, M6, and M5 at r ; 3000r g , panels (b), (d), and (f) show a similar jet morphology at the base, which is determined by the torus and its winds.With the cooling present, which impacts the torus, the jet in panel (h) opens up at smaller distances and appears wider.We notice that the backflows reach closer to the BH, since the disk scale height decreases due to cooling.

Figure 2 .
Figure2.Whereas the ambient medium does not impact the power of the jet near its launching region, disk cooling makes a profound difference; it causes the mass accretion to increase continuously and the jets to become progressively weaker relative to the mass inflow.This allows us to study the dependence of the jet shape on the jet efficiency. M and P jet are measured at r = 5r g ; Φ BH is measured at the BH event horizon.Panel (a) shows that models M7, M6, and M5 reach a steady state by t ; 15k; the mass accretion rate  M remains constant and approximately the same across the three simulations, M7, M6, and M5, because the torus and not the ambient medium dominates the BH accretion.In our model M5C, M flattens out at t ; 15k-30k and then rapidly increases at 30k  t  50k.Panel (b) shows a similar smooth decline in Φ BH across models M7, M6, and M5.In contrast, Φ BH in model M5C remains constant at 10k  t  40k before rising slightly at t ; 40k.Panel (c) shows the power of the jets, which is defined as the energy flux in highly magnetized regions where the magnetization σ > 1.The jet power is approximately the same for models M7, M6, and M5 and follows the expected trend, µ F P jet BH 2 .As we see in Figures4 and 5, at t  35k, the jet in model M5 starts to become dominated by the disk winds and disappears at t ; 60k.The power of the jet in model M5C starts flattening out at t ; 15k.The jet in model M5C is less powerful yet lives longer than the jet in model M5.

Figure 3 .
Figure 3. First self-consistent demonstration of the jet shape transition from parabolic to conical to cylindrical in GRMHD simulations (models M7, M6, and M5) at the largest to date Bondi-to-gravitational radius scale separation of 3 orders of magnitude.The disk cooling in model M5C causes the jet to start out conical and transition to a cylindrical shape with  R r 100 200 jet g ( -) well inside the Bondi radius, already at a distance of r ; 300r g .The M5C jet radius approximately matches the observed values for the 3C 84 jet, where  R r 250 jet 3C 84g at r 3C 84  350r g(Giovannini et al. 2018).This figure shows the collimation profile of the jet as a function of the distance along the jet at t = 20k when the measured  M values are similar and the P jet values are within a factor of 2 between model M5C and the rest of the models (M7-M5).Jets in models M7, M6, and M5 have a similar parabolic collimation profile extending out to r ; 200r g , beyond which they expand conically.Outside of the sphere of the BH gravitational influence, r  r B , the jets in models M7, M6, and M5 become cylindrical; as we increase the ambient medium density, the transition point moves inward.The jet in model M5C initially expands conically up to around 300r g and then becomes nearly cylindrical.

Figure 5 .
Figure5.As time progresses, disk winds grow stronger, and the jet in model M5C transitions from a quasi-conical (R ∝ r 1 ) to a quasi-parabolic (R ∝ r 0.5 ) shape near the BH but remains quasi-cylindrical (R ∝ r 0 ) at r  700r g , well inside the Bondi radius.Panels (a)-(d) show density slices in the xz-plane at 20k t 50k every Δt = 10k.During this time period, the disk progressively shrinks as it thermally collapses toward the midplane due to cooling.Panel (a) shows that this opens up the space for jet sideways expansion; this also allows the backflows to come closer to the jet base and turn the jet into a cylinder.However, panels (b)-(d) show that the wind launched from the disk eventually becomes more prominent, opposes the backflows, and strongly collimates the jet.Panel (e) shows the jet radial profile at different times (see legend; we define the jet using the μ 3 condition; see text for details).The solid lines show the profile of the jet radius, averaged over Δt = 2k, at the times indicated in the legend; the shaded regions indicate one standard deviation.Panel (a) shows that at t = 20k, the jet expands conically before turning into a cylinder at r  300r g .However, panels (b) and (c) show that at t = 30k and 40k, the jet becomes quasi-conical near the BH because the jet power remains approximately constant, but the wind power increases in time.Panel (d)  shows that t = 50k, the jet assumes a quasi-parabolic shape.Interestingly, throughout the simulation, the jet profile remains robustly quasi-cylindrical at r  700r g , i.e., well inside the Bondi radius.At even later times, t  60k, the jet becomes unstable (not shown).