Jets with a Twist: The Emergence of FR0 Jets in a 3D GRMHD Simulation of Zero-angular-momentum Black Hole Accretion

Spinning supermassive black holes (BHs) in active galactic nuclei magnetically launch relativistic collimated outflows, or jets. Without angular momentum supply, such jets are thought to perish within 3 orders of magnitude in distance from the BH, well before reaching kiloparsec scales. We study the survival of such jets at the largest scale separation to date, via 3D general relativistic magnetohydrodynamic simulations of rapidly spinning BHs immersed into uniform zero-angular-momentum gas threaded by a weak vertical magnetic field. We place the gas outside the BH sphere of influence, or the Bondi radius, chosen to be much larger than the BH gravitational radius, R B = 103 R g. The BH develops dynamically important large-scale magnetic fields, forms a magnetically arrested disk (MAD), and launches relativistic jets that propagate well outside R B and suppress BH accretion to 1.5% of the Bondi rate, ṀB . Thus, low-angular-momentum accretion in the MAD state can form large-scale jets in Fanaroff–Riley (FR) type I and II galaxies. Subsequently, the disk shrinks and exits the MAD state: barely a disk (BAD), it rapidly precesses, whips the jets around, globally destroys them, and lets 5%–10% of ṀB reach the BH. Thereafter, the disk starts rocking back and forth by angles 90°–180°: the rocking accretion disk (RAD) launches weak intermittent jets that spread their energy over a large area and suppress BH accretion to ≲2% ṀB . Because the BAD and RAD states tangle up the jets and destroy them well inside R B, they are promising candidates for the more abundant, but less luminous, class of FR0 galaxies.


INTRODUCTION
Relativistic collimated outflows, known as jets, are prevalent across many astrophysical systems of vastly different scales.The largest and most energetic ones emanate from the galaxy centers that harbor supermassive black holes (SMBHs).Active galactic nuclei (AGN), whose central SMBHs consume gas, release energy that couples to the galactic environment in a process called AGN feedback.The general consensus is that the feedback comes in two flavors: (i) Radio or kinetic mode, which occurs at lower accretion rates (L ≪ 0.01L Edd /c 2 , where L Edd is the Eddington luminosity).In this mode, powerful radio lalakos@u.northwestern.edujets dominate the feedback (see Fabian 2012; Morganti 2017, for reviews).(ii) Radiative or quasar mode, which takes place at high accretion rates in luminous AGN (L ≳ 0.01L Edd /c 2 ).Here, the radiation coupled with the surrounding gas and dust drives powerful outflows at wide angles into the galactic environment.About 1 in 10 luminous AGN can produce powerful radio jets (Sikora et al. 2007).In both cases, (i) and (ii), the radio jets can propagate into the interstellar medium (ISM) and intracluster medium (ICM), displace the gas, and inflate X-ray cavities of up to several Mpc in size (Oei et al. 2022; see Mc-Namara & Nulsen 2007, 2012, for reviews).Whereas these cavities appear empty in X-ray images, they are filled with relativistically-hot magnetized plasma that emits copious synchrotron emission in the radio band.The cavities result from AGN feedback that can offset the runaway cooling of the ICM via shocks, acoustic waves, and/or turbulent heating and thus can regulate star formation (Zhuravleva et al. 2016;Martizzi et al. 2019;Li et al. 2020).However, how SMBHs feed on the gas, power the jets, and exert feedback on their environment remains poorly understood.
The Fanaroff-Riley (FR) classification system (Fanaroff & Riley 1974) is widely used to categorize radio galaxies into two main types: FRIs and FRIIs, determined by the morphology of their radio emission.FRIs are bright near the core and grow gradually fainter with increasing distance.FRIIs are brighter near the edge of the radio lobes and grow fainter towards the core.Recently, a new classification has been established, the FR0s that morphologically lack radio emission at scales r ≳ 1 kpc (Ghisellini 2011).They are more abundant than FRIs and FRIIs (Stecker et al. 2019), hinting that they might be the dominant population in the local universe (Baldi et al. 2018; see Baldi 2023 for review).Additionally, the Fermi observatory recently detected bright γ-rays from 3 out of ≳ 100 FR0 galaxies, with 2 being an order of magnitude brighter in radio than the rest (Paliya 2021).This, along with their abundance, suggests that FR0s can significantly contribute to the isotropic γ−ray background (Stecker et al. 2019), the diffuse neutrino background (Tavecchio et al. 2018;Ackermann et al. 2022), as well as ultra-high-energy cosmic rays (Merten et al. 2021;Lundquist et al. 2021).It is not clear yet whether FR0s are part of the evolutionary stage of FRI/IIs.The SMBH mass and host environment appear similar to the FRIs, however, the jet size is smaller, r ≲ 1kpc, usually unresolved, and moves slower γ ≲ 2 (Balmaverde, B. & Capetti, A. 2006;Baldi et al. 2018).This might suggest that either the BH spin or the magnetic field strength are different enough that can affect the jet power and velocity (Baldi et al. 2018).
Spinning BHs surrounded by whirlpools of hot magnetized plasma, or accretion disks, can produce relativistic jets via the Blandford & Znajek (1977, BZ hereafter) process.The BZ process capitalizes on the fact that as the hot plasma in the disk accretes on a spinning BH, it brings with it poloidal (pointing in the Rand z-directions) magnetic flux.BH rotation drags the inertial frames via the Lense & Thirring (1918) effect, which twists the poloidal flux in the vicinity of the BH, causing it to build up a strong toroidal component, which in turn powers collimated Poynting-flux dominated outflowsthe twin polar jets (see Tchekhovskoy 2015).The rotation of the disk launches disk winds (via a Blandford & Payne 1982 like process) that collimate the jets into small opening angles and enable the jets to accelerate (Komissarov et al. 2008;Tchekhovskoy et al. 2009).This suggests that rotation is an important ingredient for accretion, jets and outflows, and their feedback on the environment.In this work, we aim to reveal the importance of ambient gas rotation on the jets and their feedback.
To understand how the jets escape out of the BH sphere of influence, one must follow them from their formation at the BH event horizon to their interaction with the ISM or ICM.For this, one needs to model the BH feeding by following the gas infall from the ISM to the BH event horizon.These are extremely arduous tasks due to the enormous scale separation of the problem.The gas originates at the edge of the BH sphere of influence, the Bondi radius (Bondi 1952), which is significantly larger than the BH gravitational radius, R g = GM BH /c 2 = 5 × 10 −5 pc × (M BH /10 9 M ⊙ ).Here, M BH is the SMBH mass and c ∞ is the ISM sound speed.For a typical ISM temperature, T ∞ ∼ 1 keV, the scale separation between R B and R g reaches 6 orders of magnitude (Russell et al. 2015(Russell et al. , 2018)).Adding to the challenge is that the ICM and dark matter halos reach scales of ∼ 100 kpc, marking the scale separation of the full problem ∼ 9 − 10 orders of magnitude.Dedicated galaxy simulations are capable of following gas flows from the dark matter halo down to ∼ 1 kpc scales (Anglés-Alcázar et al. 2015, 2017), and state-of-the-art hyper-refined Lagrangian simulations can even reach sub-pc scales (Anglés-Alcázar et al. 2021).However, to bridge the "last mile" of the scale separation and connect the smallest scales in galaxy simulations to the event horizon, general relativistic magnetohydrodynamic (GRMHD) simulations are essential.
A popular approach for bridging the scale separation is the Bondi model (Bondi 1952;Shapiro & Teukolsky 1986), which approximates BH accretion as a spherically symmetric hydrodynamic flow.Although an elegant and simple approximation, it does not allow for jets.In this work, we add the minimum ingredients to the Bondi model that enable the study of BH-powered jet formation and propagation: we retain the zero angular momentum of the ambient gas but add BH rotation and ambient vertical magnetic field.We also consider the system at sufficiently high resolution and in full 3D to resolve the jet propagation and allow for the development of nonaxisymmetric instabilities in the jets.
GRMHD simulations found that magnetized accretion can accumulate large-scale poloidal magnetic fields on the BH to the point that the fields become dynamically important, i.e., able to obstruct gas infall, at which point the system enters the magnetically arrested disk state (MAD, Bisnovatyi-Kogan & Ruzmaikin 1974, 1976;Narayan et al. 2003;Igumenshchev et al. 2003;Igumenshchev 2008) that can launch jets with energy efficiencies exceeding 100%, i.e., whose power exceeds the accretion power (Tchekhovskoy et al. 2011).In this state, the accretion continuously brings in the magnetic flux to the BH, thereby flooding both the BH and inner disk with the vertical magnetic flux.At the same time, the dynamicallyimportant magnetic flux periodically erupts from the BH, rips through the disk, and escapes.In the presence of sufficient angular momentum, the inward advection of the magnetic flux wins, and the BH is always flooded with the magnetic flux (Tchekhovskoy & McKinney 2012).
In contrast, zero angular momentum accretion was found to not be conducive to powerful, stable jet production.Analytic models suggested that zero angular momentum accretion produces weak jets (with energy efficiency ≲ 1%, Das & Czerny 2012).Kwan et al. (2023, hereafter K23) used GRMHD simulations to model the accretion of spherically-symmetric magnetized gas on a rapidly spinning BH and found that zeroangular accretion could not remain in a MAD state for a prolonged period of time.The jet efficiency only transiently surpassed 100% due to the BH eventually losing the large-scale BH magnetic flux powering the jets (see also Gottlieb et al. 2022a, for a related phenomenon in collapsars) and entering the SANE state ("standard and normal evolution", Narayan et al. 2012), in which the large-scale magnetic flux is subdominant and the jets are weaker or non-existent.Ressler et al. (2021, hereafter R21) performed 3D GRMHD simulations of spherically symmetric gas accretion onto a rapidly-spinning BH (a = 0.9375), with a Bondi-to-event horizon scale separation of R B /R g = 100.1They varied the angle between the ambient magnetic field and BH spin directions and found that outflow efficiencies, in most cases, did not exceed 100%: the accretion flow had a hard time reaching and staying MAD.Based on the results of their simulations, R21 predicted that for R B /R g ≳ 800 all jets powered by low angular momentum accretion will end up falling victim to the kink instability inside the Bondi radius, as is the case for realistic systems, R B /R g ≳ 10 5−6 (e.g., SgrA* or M87).
Here, we evaluate the ability of low angular momentum accretion to power large-scale jets for an order of magnitude larger scale separation than has been possible until now.Lalakos et al. (2022, hereafter L22) simulated 3D GRMHD accretion of rotating ISM for the Bondi-to-event horizon ratio of R B /R g = 10 3 .Here, we follow the L22 approach, but consider zero angular momentum accretion.As in L22, we include a rapidly spinning BH, a = 0.94, and large-scale vertical magnetic field, which is aligned with the BH spin vector.To study the self-consistent jet formation and propagation to distances well outside the Bondi radius, we use adaptive mesh refinement (AMR) to ensure sufficient resolution of the tightly collimated jets.
More generally, we aim to understand the basic ingredients necessary for the formation of stable jets or, conversely, what it takes to "break" them.Namely, we aim to answer questions such as: Is gas angular momentum (and the formation of an ac-cretion disk) needed to maintain jet stability?Is there a critical power above which the jets manage to escape the BH sphere of influence?How does the ambient medium trigger jet instability?What are the observational signatures of unstable jets?
Throughout the paper, we work in units of G = M = c = 1.For conciseness, we sometimes measure the time in the unit of 1000R g /c, which we denote as 1k: e.g., 7.5 × 10 4 R g /c ≡ 75k.In Sec. 2, we describe our simulation setup and choice of the initial physical parameters.In Sec. 3, we present the properties of the accretion flow and jets as measured near the BH.In Sec. 4, we show how the jets can escape out of the Bondi sphere and discuss their stability.In Sec. 5 we show how the changes in the behavior of accretion flow angular momentum can lead to the destruction of the jets.In Sec. 6, we compare internal and external kink instabilities.Finally, in Sec. 7, we summarize and discuss our results.

NUMERICAL METHOD AND SETUP
We carry out our simulations using the general relativistic magnetohydrodynamic (GRMHD) code H-AMR that includes advanced features such as graphics processing unit (GPU) acceleration, AMR, and local adaptive timestepping (Liska et al. 2022).We initialize the simulation by placing a BH inside a uniform static ambient medium.We consider a rapidly spinning BH, with dimensionless spin a = 0.94, to favor jet launching.To allow the system to evolve naturally, we avoid prescribing the conditions inside R B : we carve out an empty cavity inside the Bondi sphere (r < R B ) and place the ambient gas of uniform rest-mass density, ρ = ρ ∞ , outside the sphere (r ≥ R B ).
We choose the ambient sound speed, c ∞ , to achieve the desired scale separation, R B = 10 3 R g (eq.1).This is similar to L22, but here we set the ambient medium angular momentum to zero.In this work, we focus on low-luminosity BH accretion systems.Because of this, we do not include any radiation effects (e.g., radiative cooling), which are important at higher mass accretion rates.Thus, our simulations are non-radiative and scale-free: the simulation results trivially rescale to any value of ρ ∞ .As appropriate for AGN environments at the Bondi radius, we adopt a monatomic non-relativistic (Γ = 5/3) gas with an ideal gas equation of state, p g = (Γ − 1)u g , where p g and u g are the pressure and internal energy density of the gas.
Outside the Bondi radius, we include a large-scale vertical magnetic field in the direction parallel to the BH spin.Asymptotically far away, we set ⃗ B = B 0 ⃗ a. Closer to R B , we deform the field such that the radial component of the magnetic field, B r , smoothly vanishes towards the edge of the cavity and no field enters the cavity, r ≤ R B .To achieve this, we adopt the magnetic vector potential, A φ ∝ (r 2 − R 2 B ) sin 2 θ.We normalize the strength of the magnetic field such that the thermal to magnetic pressure ratio is β = p g /p m = 100 asymptotically far away (this also ensures that β ≥ 100 everywhere).Here, p m = b 2 /8π is the magnetic pressure, where b 2 = b µ b µ , and b µ is the comoving contravariant magnetic field 4-vector (defined in Appendix A) In order to break axisymmetry and provide the seeds for the growth of the magneto-rotational instability (MRI, Balbus & Hawley 1991), we include random pressure perturbations at a 2% level in the initial conditions.
We note that the cavity is, technically, not entirely empty because GRMHD codes cannot handle vacuum.To prevent ρ and u g from becoming extremely low or negative in highly magnetized regions (i.e., the jet launching regions), our numerical scheme adopts the following density floors.If at any point in the simulation the densities drop below the floor values, ρ < ρ fl = max[b 2 /(60π), 10 −7 r −2 , 10 −20 ] and/or u g < u g,fl = max[b 2 /(3000π), 10 −9 r −2Γ , 10 −20 ], then we add mass and/or internal energy, respectively, in the drift frame until the floor values are reached (see Appendix B3 of Ressler et al. 2017).
We construct our grid in spherical polar coordinates, r, θ, and φ.The radial grid is uniform in log r, and r spans 0.83R H ≤ r ≤ 10 6 R g .There are 6 cells inside the event horizon, , and the radius of the outer boundary is larger than the light travel distance in a simulation duration: these ensure that both the inner and outer radial boundaries are causally disconnected from and cannot influence the solution.The polar and azimuthal grids are uniform in the θand φdirections, and span 0 ≤ θ ≤ π and 0 ≤ φ ≤ 2π, respectively.We use outflow, transmissive, and periodic boundary conditions in the radial, polar, and azimuthal directions, respectively (Liska et al. 2022).The base-grid resolution is N r × N θ × N φ = 448 × 96 × 192 cells in the r-, θ-, and φ-directions.At r ≥ 6.5R g , we activate 1 level of static mesh refinement (SMR): this doubles the resolution in each dimension, leading to an increased effective resolution of 896 × 192 × 384.On top of the SMR, we also include 2 additional adaptive mesh refinement (AMR) levels that we dynamically activate during the run to ensure sufficient resolution to resolve the collimated jets and cocoons at large radii as we describe in Appendix B. The maximum effective resolution in the jets can therefore reach 3, 584 × 768 × 1, 536 cells.
In order to avoid the polar singularity interfering with highly magnetized jetted regions, we tilt the entire BH-gas system by 90 • relative to our computational grid: as a result, the BH rotational axis is perpendicular to the polar singularity of our computational grid.However, for the sake of narrative simplicity, below we present our results as if we did not perform the tilt: for presentation purposes, we direct the z-axis along the BH spin vector and count off the polar angle, θ, from the BH spin direction.

MAD PROLOGUE
The simulation starts with the constant-density ambient gas at rest, located outside of the empty cavity, r ≥ R B (Sec. 2).The pressure gradient at the edge of the cavity, along with the gravitational attraction from the BH, push the gas inward.To study the mass flow in our simulation, we define the rest-mass (RM) component of the energy flux density, whose angular integral gives us the net energy rate, or power, where dA = √ −gdθdφ is the differential surface element, g = |g µν | is the determinant of the metric, and u µ is the coordinateframe contravariant proper velocity vector.In a steady state, ĖRM is conserved and independent of radius.We evaluate BH mass accretion rate as the negative of the rest-mass power, Ṁ ≡ − ĖRM (r = 8R g )/c 2 , which we measure at r = 8R g to avoid potential contamination by the density floors near the event horizon.
Figure 1(a) shows the time-dependence of Ṁ, which peaks at approximately the analytic Bondi prediction, at about a free-fall time, , after the beginning of the simulation.Here, for our choice, Γ = 5/3, we have λ s = 1/4 (Shapiro et al. 1976;Di Matteo et al. 2003).That Ṁ reaches the simple analytic Bondi prediction, ṀB , is not entirely surprising because, in the absence of any external angular momentum supply, the gas never encounters a centrifugal barrier, which would inhibit the accretion relative to the Bondi expectation.The infalling gas drags inward the large-scale vertical magnetic flux, and some of it makes it all the way to the event horizon, resulting in the increase of the absolute BH magnetic flux, Here, the integral is over the entire event horizon of the BH, and the factor of 0.5 converts it to a single hemisphere.
which measures the strength of the magnetic flux relative to the onslaught of the infalling gas; here, to clarify the units, we explicitly include the length and velocity scaling factors, and Ṁ τ is the rolling average of the mass accretion rate over the time interval of τ = 3k.This time interval is sufficiently long to average over the strong Ṁ oscillations in the MAD state.We also consider a signed magnetic flux through the Northern hemisphere, as well as its normalized form, ϕ N , defined analogous to eq. ( 6).
As the infalling gas continuously drags more of the ambient magnetic flux inwards, the BH magnetic flux grows in strength and ϕ BH steadily increases.Similar to eq. ( 3), we compute the various components of the energy flux using the stress-energy tensor: We are interested in the radial energy flux and set κ = r and λ = t.Thus evaluated, eq. ( 8) gives us the total radial energy flux, f TOT .The electromagnetic (EM), thermal (TH), and kinetic energy (KE) flux density components of the total energy flux, respectively, are: where u r is the radial contravariant 4-velocity component and u t is the temporal covariant velocity component, which is a conserved quantity for point masses.We also define the components of power (integral versions of the flux densities), eqs.( 9)-( 11), following eq.( 3): where # = EM, TH, KE, or RM.The integral version of eq. ( 9) gives us the jet power, L j = ĖEM (r = 8R g ), which we measure through a sphere of radius r = 8R g .Note that L j includes the power of both jets.Figure 1(c) shows that L j abruptly increases at t ≃ 24k once ϕ BH exceeds a critical value, ϕ BH ≃ 15, and the jets form.Figure 1 (d) shows that the outflow energy efficiency, η ≡ L j /⟨ Ṁc 2 ⟩ τ , which is jet power measured in units of BH accretion power: at t ≃ 24k, it is not very high yet, but it is still significant, η = 0.1 ≡ 10%, where from now on we will express all fractions in terms of percent.Although at this time the jets have not yet reached the maximum efficiency, they already exert significant feedback on the accretion flow and suppress BH accretion rate relative to its peak by nearly an order of magnitude, Ṁ ∼ 0.2 ṀB , via injecting the energy into the accretion flow and partially unbinding it.This marks the end of the initial transient phase, during which the system settles into a quasi-steady state.
Figure 1(b) shows that the dimensionless BH magnetic flux grows until it saturates around ϕ BH ≳ 50 at t ≳ 31k.At this time, the BH magnetic flux becomes dynamically important, i.e., the magnetic pressure can withstand the onslaught of the total momentum flux of the infalling gas.This signals the formation of a MAD (Tchekhovskoy et al. 2011).The jets produced during the MAD state attain the maximum power for a given Ṁ, increasing their chance to reach the Bondi scale and produce feedback.This suppresses the mass accretion rate to Ṁ / ṀB = 1.5%, a staggering reduction by nearly 2 orders of magnitude from the peak, and the analytical Bondi (1952) expectation, ṀB .The average jet power, ⟨L j ⟩ / ṀB c 2 ≃ 2.9%, translates into extremely high jet efficiency, ⟨η⟩ ≃ 190%: this implies that the BH energy output in the form of jets exceeds the energy input in the form of the rest-mass energy.This is typical for rapidly spinning (a ≥ 0.9) MAD BHs, whose rotational energy is extracted -by the continuous winding of magnetic fields on the event horizon -faster than the accretion can replenish it (Tchekhovskoy et al. 2011;Tchekhovskoy & McKinney 2012;McKinney et al. 2012).Figure 1(b) shows that at t ≳ 65k the normalized magnetic flux drops below the characteristic MAD value, ϕ BH = 50, the flow exits the MAD state, and the jets become progressively weaker (Fig. 1c) until their complete destruction at the event horizon at t = 78k, when η ≤ 0.1% (Fig. 1d).With no strong jets to obstruct the accretion, the mass accretion rate (Fig. 1a) increases from Ṁ/ ṀB ∼ 1.5% in the MAD state to Ṁ/ ṀB ≲ 10% in the SANE state.

A TWISTED INTERLUDE
In this section, we explore the large-scale stability and propagation of the jets, and how they reach and impact the BH feeding scales at r ≳ R B .Some jets can traverse vast distances, sometimes larger than the galaxy itself (e.g., Cygnus A, Perley et al. 1984), whereas others appear to perish within the galaxy (e.g., M87, Biretta et al. 1991).This involves a scale separation of ∼10 orders of magnitude from the BH to the outskirts of the galaxy.Along the way, the jets survive in the face of many obstacles, each of which can become a major dissipationtriggering event (e.g.internal, external, and/or recollimation shocks, MHD instabilities, magnetic reconnection, etc.), where a significant fraction of the jet energy can transform into radiation, and the jet can become locally or globally disrupted.Indeed, magnetized jets are prone to current-driven instabilities, with the kink instability being perhaps the most destructive (Bateman 1978;Lyubarskij 1992;Eichler 1993;Begelman 1998;Lyubarskii 1999;Narayan et al. 2009).Kink instability acts on the jet core, dislodging it from its original location and twisting it into a helical structure.In ambient media with a flat density profile, ρ ∝ r −α , where α < 2, the kink instability is inevitable (Bromberg & Tchekhovskoy 2016, hereafter BT16): in such a medium, a jet with a constant opening angle would displace progressively more and more gas as it propagates out.The jet responds by becoming progressively more collimated.As we will see below (Sec.4.1), this makes it more susceptible to the kink instability that can cause the jets to fall apart.
To visualize the large-scale morphology and evolution of the jets, we have created a 3D volume-rendered movie.We present a sequence of key snapshots in Fig. 2, taken at the times indi- 3).At t ≳ 65k, the system exits the MAD state, the jets become weaker, ⟨η⟩ ≲ 15%, and the mass accretion rate higher, Ṁ / ṀB ≲ 5% (red-shaded region, 69k ≤ t ≤ 95k).The jets sometimes completely shutoff and fall apart (e.g., at t = 78k and 96 − 102k; Sec. 4) when the normalized Northern hemisphere BH magnetic flux vanishes, ϕN = 0 (thin line in panel b), η drops (to ≲ 0.1%, panel d), and the reduced jet feedback allows higher Ṁ/ ṀB ≲ 10% (panel a).However, even such weaker jets can suppress Ṁ similar to the MAD state, Ṁ / ṀB ≃ 1.7%, for extended periods of time if they continuously reorient (orange-shade region 110k ≤ t ≤ 130k; Sec. 5 and Fig. 4).In all panels, horizontal dashed lines show time average values, which are also labeled with call-outs.cated by the vertical lines in Fig. 1.The jets are launched at t = 24k (Fig. 1c), and by t ≃ 28k they reach the Bondi sphere and start interacting with the ambient medium (not shown in Fig. 2). Figure 2(a) shows that by t = 58k, the high-power MAD jets (red) drill their way through the ambient medium out to r ≃ 4000R g = 4R B .As the jets ram into the ambient medium, they form strong twin bow shocks at the jet heads.
The shocks increase the internal energy at the heads, which typically renders these regions observable as hotspots.At the head, jet material spills sideways and flows back to create the inner cocoon (orange), and the strong bow shock, caused by the jet's relativistic motion, heats up the ambient gas to create the outer cocoon (blue).The jets are mostly straight apart from the bends in the outer half of the jets (r ≳ 2000R g ). Figure 2(b) shows that at later times, t = 65k, these bends become more pronounced, e.g., one of them twists the jet into a knot, as seen at r ≃ 2000R g in the left jet.
Soon after this, the MAD state comes to an end: by t = 71k, the jet power drops by a factor of ∼ 4 (Fig. 1c). Figure 2(c) shows that the jets become more helical: the kink instability affects a larger fraction of the jet length, with the first jet bend showing up already at r ≲ 1000R g .Figure 2(d) reveals that by t = 78k the jets get globally destroyed: at this time, the power of the jets essentially vanishes (Fig. 1c), and they no longer focus their energy into the twin cocoons.At later times, Figure 2(e)-(f) shows that short-lived jets still form close to the BH but fall apart around, or even inside, the Bondi radius, and starve the cocoons of energy.With the energy supply via the jets mostly shut off, the mildly magnetized cocoons become relics, with a combination of leftover momentum and buoyancy, and occasional mergers with smaller cavities inflated by the wobbling jets, driving them outwards.At large scales, although weakened, bow shocks are still present and outrun the cocoons while injecting energy and momentum into the ambient gas.

Jet Stability Criterion
To quantitatively analyze the development of the kink mode in the jets, we use eqs.( 3) and ( 9)-( 11) to express the ratio of total to rest-mass energy fluxes, which gives the maximum Lorentz factor a fluid element can attain if all of its EM and thermal energy fluxes were converted into the kinetic energy flux.To obtain eq. ( 13), we used the facts that f TOT = f KE + f TH + f EM and and introduced the relativistic gas enthalpy per unit mass, where h is the non-relativistic gas enthalpy per unit mass (i.e.without the rest-mass contribution), and magnetization, where the approximate equalities in eqs.( 14) and ( 16) are accurate at large distances from the BH (Chatterjee et al. 2019).Below, unless stated otherwise, we will use the right-hand sides of eqs.( 13), (15), and ( 16) as expressions for µ, h, and σ, respectively.
Jets consist of an inner core of cylindrical radius, R c , dominated by the poloidal comoving magnetic field component, b p , and an outer region, dominated by the toroidal comoving field component, b tor (BT16).Here, R c = rsinθ c corresponds to the opening angle of the jet core, θ c , at distance r, and we measure both b p and b tor in the fluid frame (see Appendix A).The kink instability growth timescale for a local fluid element is approximately the time it takes for an Alfvén wave to travel around the circumference of the jet at the Alfvén velocity, We define the kink instability growth time scale, as measured in the lab frame: where we used the jet Lorentz factor, γ, to convert the timescale from the fluid-frame to the lab-frame.The ratio, b p /b tor , is the magnetic winding, and R × b p /b tor is the pitch, where R = rsinθ is the cylindrical radius of the jet location in question.The factor, η kink ≃ 5 − 10, is the poorly constrained prefactor that enters the kink timescale for the jet to become considerably deformed (Mizuno et al. 2009(Mizuno et al. , 2012;;Bromberg et al. 2019).The approximate equality in eq. ( 18) evaluates the timescale at the edge of the jet core, where by definition, b p = b tor , and the value of the pitch is simply R c .The kink instability develops in the fluid frame, and the available time for its growth is the propagation time of the jet fluid from the BH, moving at velocity v.In the lab frame, this dynamical time is: where we assume constant jet velocity, v, ignoring the initial acceleration of the jet, and that the jet travels only in the radial direction (not true in a bent jet).The ratio of the kink to dynamical timescale (eqs.18 and 19) gives us the stability parameter: The jet is stable at Λ > 2, and unstable at Λ < 2 (BT16).Calculating the radius of the jet core is a non-trivial task since the jet core can displace itself from the z-axis in a non-axisymmetric fashion.To compute R c , we first evaluate the solid angle subtended by the jet core, Ω c , which we identify as satisfying both µ > 3 and b p ≥ b tor conditions.Then, we can calculate the opening angle from Ω c = 2π(1 − cosθ c ), and then R c = rsinθ c .Without loss of generality, we focus on one of the twin jets, the one pointing in the direction of the BH spin, ⃗ a.

Jet Bends Set the Speed Limit
Figure 3(a1) shows a transverse projection of the jet proper velocity, γv, at t = 58k (the same time as in Fig. 2a).Here, we [panels (b)]: Whereas at the earlier time (t = 58k, panels a1 and b1) only the outer jet, r ≳ 2000Rg, shows bends and satisfies the kink instability criterion, Λ ≲ 2, at the later time (t = 71k, panel a2 and b2) most of the jet, r ≳ 800Rg, shows bends and satisfies the criterion.Jet bends cause the core-average proper jet velocity γv (magenta) to significantly decrease from γv ≃ 3 (panel b1) to γv ≲ 2 (panel b2) and result in a comparable decrease in Λ (eq.20).Although the jet radius, Rj (dark red), also decreases, this noticeably affects neither the core radius, Rc (gold), nor core-average Alfvén velocity, vA (green).[panels (c)]: Angle-integrated fluxes, or power, through the jet (µ > 3) normalized to the time-average total MAD power, ⟨ ĖMAD⟩.The normalized total jet power, ĖTOT (black), remains approximately constant in stable jet regions.The kinetic energy component, ĖKE (light blue), gradually increases at the expense of the EM component, ĖEM (purple), which dominates the total energy power until the instability sets in at r ≳ 2000Rg in panel (c1) and r ≳ 800Rg in panel (c2).Panels (a1) and (a2) show that in the unstable regions, the jet develops one or more bends.Whereas just before each substantial bend the thermal component (orange-red) grows up to equipartition with the EM component (purple), ĖTH ∼ ĖEM, after each bend all energy power components drop, followed by a gradual increase until the next bend (e.g., for bends at r/Rg ≃ 800, 1300, 2000, 2500, and 3000 in panel c2): we attribute this to the "accordion" effect -the jet bends allow the jets to elongate and rarefy longitudinally.The rest-mass power ĖRM (orange) remains subdominant to the rest of the components.[panels (d)]: Jet-average proper velocity, γv (light blue, obtained via γ ≡ ĖKE/ ĖRM), increases as the jets accelerate at the expense of the decreasing jet-average magnetization, σ ≡ ĖEM/ ĖKE (purple), before decelerating at jet bends.Jetaverage enthalpy, h ≡ ĖTH/ ĖKE (orange-red), increases at jet bends due to energy dissipation.The µ ≡ ĖTOT/ ĖRM parameter (black) decreases due to the jet losing its energy to the ambient medium (via ambient gas displacements and shocks).
average γv along the line of sight (which is perpendicular to the direction of BH spin, ⃗ a), and weigh it with the magnetization, σ, to highlight the internal structure of the jet.At this time, most of the jet is free of significant bends, except towards the head, at r ≳ 2000R g = 2R B .As the jet bends, its proper velocity drops from γv ≃ 3 (dark blue) to γv ≃ 1.7 (green), which corresponds to γ ≃ 2: the slower the jet, the easier for it to make sharp turns (see below).The jet maintains its transrelativistic proper velocity until it reaches the head (hotspot), splashes sideways and backward, and creates the backflows (brown) at γv ≲ 1. Figure 3(a2) shows the jet after the MAD state has ended, at t = 71k (same time as in Fig. 2c).By this time, the average jet power has decreased by a factor of ∼ 4, compared to Fig. 3(a1).Weaker jets are less rigid, thus their interaction with the ambient medium can bend them more easily.Indeed, the jet develops bends that are much stronger and closer to the BH than in the MAD state.These bends force most of the jet to become mildly relativistic, γv ≤ 1.7.Even tiny, visually imperceptible, bends can have dramatic consequences for relativistic jets.
For relativistically-magnetized jets, the consequences can become particularly dire when the jets become super-fast magnetosonic, i.e., when they outrun their own fast magnetosonic waves, which for cold flows translates to a proper velocity, (Tchekhovskoy et al. 2010).If we introduce the fast Mach number, M F = γv/u F ≈ γv/σ 1/2 , we can write the fast wave Mach cone opening angle as where the first approximate equality applies when the jets are in the highly super-fast magnetosonic regime, M F ≫ 1, and the latter applies in the cold limit.Why do relativistic jets have a hard time making turns?This is because in order for a jet to navigate a turn smoothly, the jet must be able to anticipate that the turn is coming.Here, θ Mach gives the jet's "field of vision" angle, i.e., the maximum angular deviation within which fast waves can mediate the total pressure perturbations.Thus, the jets cannot bend by more than the fast Mach cone opening angle, or, equivalently, the jets cannot exceed the "speed limit", u max , where the latter approximate equality applies in the cold limit.
The stronger the bend (larger θ bend ), the smaller the speed limit, u max , which is the maximum allowed jet proper velocity, γv, for which the jets avoid the development of internal shocks.If γv exceeds u max , the jets develop oblique internal shocks that discontinuously reduce the jet velocity.Bends are not the only mechanism through which jets can develop shocks and decelerate.For instance, if the jets conically expand into a medium, they accelerate rapidly and soon become super-fast magnetosonic.At some point, the pressure of the hot cocoons engulfing the jets starts to dominate over the jet internal pressure and squeezes the jets.As a result, the jets develop collimation shocks behind which the jet material slows down.Magnetic stresses operating downstream of the shocks introduce additional compression forces that cause the jets to pinch and form narrow nozzles in which the kink instability can grow more efficiently and dissipate energy (e.g. Figure 3(b) shows the radial profiles of proper, γv (magenta), and Alvén, v A (green), velocities, which we have angleaveraged over the jet core of radius R c (shown in gold, see Sec. 4.1 for definition).For comparison, we also show the radius of the jet, R j (dark red).Using eq. ( 20), we can now compute the stability parameter, Λ (blue): here, the line thickness indicates the uncertainty range of η kink = 5 − 10 in the definition of the stability parameter.That Λ ≳ 2 at r ≲ 1500R g in Figure 3(b1) implies that according to the stability criterion the jet is expected to be stable: Fig. 3(a1) shows that at these distances, the jet appears to be mostly straight.At larger radii, however, we have Λ ≲ 2, suggesting an instability: consistent with this expectation, the jet develops bends at r ≳ 1500R g .At the later time, Fig. 3(b2) shows that Λ ≲ 2 throughout most of the jet, already at r ≳ 500R g , indicating that according to the stability criterion, we expect that most the jet is susceptible to the kink instability.Fig. 3(a2) shows multiple strong bends throughout the length of the jets.We note that strictly speaking, eq. ( 20) is only valid for small perturbations (i.e., in the limit of small jet bends) and not applicable once the bends become extreme (see Sec. 6).
Figure 3(b1) shows that in the stable region, at r ≲ 1500R g , the jet core is relativistic, γv ≃ 3, with an average radius of R c ≃ 20R g .The multiple dips in γv correspond to recollimation or oblique shocks in the jet, at r/R g ≃ 800, 1300, and 1500.At these locations, we see that the velocity drops to γv ≃ 2, and the jet total and core radii dip; Figure 3(a1) shows that the jet tends to develop bends at the same locations.At r ≳ 2000R g , the jet displays its largest bend: at the beginning of the bend, the jet velocity dips to γv ≃ 1.5 and the core radius dips to R c ≲ 10R g .As the jet bend proceeds, the core radius gradually increases to R c = 50R g likely reflecting the dissipation of the toroidal field; the velocity also recovers and levels off at γv ≳ 2. Neither the jet radius R j , which remains ∼ 3 − 4 times wider than the core, nor the Alfvén speed, which remains roughly constant at v A /c ∼ 0.5 − 0.7, appear to be affected by the bend.
At the later time, Fig. 3(b2) shows that the jet core moves at a transrelativistic speed, γv ≲ 2, and both R c and Λ show multiple dips, likely associated with oblique shocks (e.g., at r/R g ≃ 400 and 600), at which the stability parameter dips down to an unstable level, Λ < 2, indicating the presence of the instability.Figure 3(a2) shows that the jet disrupts at r = 800R g , and the velocity decreases to γv ≲ 1.At the same time, both R c and R j spike due to the r = constant cross-section slicing the jets at an angle and making jet features appear wider.Even at such high resolutions as ours, the simulation resolves the transverse extent of the jet core (∼ 30R g ) out to r ≃ 1500R g by 5 cells.The simulation resolves the transverse extent of the jet (∼ 100R g ) out to r ≃ 5000R g by 5 cells, although at such large distances, it might lose some details of its internal structure.

Energy Partition and Dissipation
Figure 3(c)-(d) shows angle-average quantities, over the entire jet cross-section, which we define as the relativistic region, µ > 3, to avoid the contamination caused by the mildly magnetized cocoon and weakly magnetized ambient gas.(This is in contrast to Fig. 3b that considered averages over the jet core.) Figure 3(c1) shows various components of the radial angleintegrated jet energy flux, essentially the energy rate, or power.We normalize all of the power components by the total power flowing through the entire jet cross-section (µ > 3) at r = 8R g , averaged over the duration of the MAD state, ĖMAD ≡ ⟨ ĖTOT (r = 8R g )⟩ MAD ; here ⟨. . .⟩ MAD denotes the time-average over the MAD state, which we highlighted in purple in Fig. 1 2 .The total jet power, ĖTOT (black line in Fig. 3c1), remains approximately constant at r ≲ 2000R g , and drops slightly thereafter, just as the jet starts to bend, as seen in Fig. 3(a1).The jet electromagnetic power, ĖEM (purple), dominates the jet power at r ≲ 2000R g , whereas the thermal component of the power, ĖTH (orange-red), steadily increases until, at r ∼ 2000R g , it reaches equipartition with ĖEM , a tell-tale sign of energy dissipation (BT16).The increase in ĖTH just as the jet starts to bend most likely comes from the dissipation of EM and kinetic energy into heat due to kink, mixing instabilities, and/or shocks caused by the bend (see Sec. 4.2). 3 The kinetic energy power component, ĖKE (light blue), steeply rises at r ≲ 100R g , while the rest-mass power component ĖRM (orange) is subdominant compared to the other components.
At r ≲ 2000R g , ĖKE remains approximately constant because the jet becomes cylindrical (due to external collimation by the constant density ambient medium).However, once the jet becomes unstable to the kink instability and develops strong bends, ĖKE drops.At the later time, t = 71k, Fig. 3(c2) shows that ĖTOT steeply drops at r ≃ 800R g , right where the jet makes a sharp 90−degree turn, as seen in projection on Fig. 3(a2).In fact, all of the energy power components experience similarly sharp drops, since we are only accounting for the radial component of the fluxes, which vanishes.At the location of such drops, which correspond to strong bends, the jet is disrupted and the energy rate does not increase, meaning it is lost.This is anticipated since (i) when the jets fall apart, they can mix with the ambient medium, lower the specific energy, and thus, they won't be picked out by our jet-criterion µ ≥ 3, and (ii) the jet head is spending energy into displacing the surrounding gas, while also losing a fraction of its energy into the formation of the cocoons.
In Fig. 3(d1) we plot the jet area-averaged radial profiles of µ parameter (black), the magnetization, σ (purple), the proper velocity, γv (light blue), and specific gas specific enthalpy, h (orange-red): To compute the jet-average values we expressed them through the ratios of jet power components (i.e. using Ė# in place of f # in the definitions 13 -16).Figure 3 (d1) shows that the jet accelerates to relativistic velocities, γv ≃ 2 at r ≲ 100R g and γ ≃ 3 at r ≲ 1000R g .At r ≲ 2000R g the jet is still magnetized, σ ≳ 2, the proper velocity remains (on average) relativistic, γv ≳ 3, and the enthalpy, h, slowly increases, similar to ĖTH .Once the jet becomes unstable and develops bends, the magnetic energy dissipates into heat, until σ ≲ h, i.e. until the magnetic energy comes into equipartition with the thermal energy.Kinetic energy is also reduced, with the jet decelerating down to the transrelativistic value γv ≃ 1.7, as the jet needs to slow down in order to be able to navigate jet bends (see Sec. 4.2).In Fig. 3(d2), we see a similar behavior as in Fig. 3(d1), but with the kink instability happening much closer to the BH, already at r = 800R g .Both the magnetization and proper velocity drop, σ ≲ h and γv ≃ 1.7, and enthalpy h increases.
Figures 4(a) show 3D volume-rendered images of the lowdensity (purple) jet and the high-density accretion flow (orange).Figure 4(a1) shows the system at t = 50k, when it is in the MAD state: the jets are powerful and mostly bend-free, and the accretion flow near the BH shows a disk-like structure.Figure 4(a2) shows the system after it exited the MAD state, where the jet power has decreased by a factor ∼ 4. The weakened jets appear to have bends originating from wobbles near their launching region.Figure 4(a3), shows jets with extreme bends, misaligned from the BH spin (z−)axis.The accretion flow is also significantly misaligned, compared to Figs. 4(a1, [panels (a)]: We show 3 distinct 3D volume-rendered images of density.Panel(a1) captures the system during the MAD state, at t = 50K, where the jets (purple) follow the direction of the BH's spin, and the accretion disk (orange) lies perpendicularly to the jets.Panel(a2) captures the system in the BAD state, at t = 71k, with weakened jets showing strong bends near the BH and getting significantly kinked.Panel(a3) captures the system in the RAD state, at t = 135k.The jets are on average 1 order of magnitude weaker and misaligned from the spin direction, The disk is also significantly tilted, at about π/2, and processing, causing significant bends and wobbles in the jets.[panel (b)]: We plot the density-weighted specific angular momentum normalized to the local Keplerian value in a space-time diagram.Angular momentum builds around the BH promptly after freefall (tff ≃ 22k) with an average disk of size 10 − 50Rg, chosen for l ≳ 0.3lK.Between, 65k ≤ t ≤ 100k the accreted angular momentum forms intermittent and tilted disks with l ≤ 0.3lK which results in intermittent and wobbly jets.After t ≃ 100k angular momentum starts building up, this time reaching larger radii ≃ 200Rg, however, near the BH the intermittent jets push most of the gas away while oscillating wildly, stopping the angular momentum from coherently adding up and stabilizing the jets.[panel (c)]: We show the tilt and precession angles of the disk and north jet, measured at a distance of r = 50Rg.During the MAD state (purple shaded), both the tilt of the disk Pd ≲ π/16 (orange) and jet Pj ≲ π/16 (purple), as they are aligned with the BH spin.The disk precession angle (orange circles), Td, is illdefined as the disk tilts about the z-axis, while the jet precession angle (purple crosses), Td ≃ π/2, which means it propagates at a preferential direction.Once we transition from the MAD to the BAD state (red shaded), at t ≥ 65k, the disk angular momentum has decreased, its tilt increases, reaching up to π/2, and rapidly precesses.The jet tries to follow and forms large bends (see Fig. 4a2) that lead to its destruction.Between 95k ≲ t ≲ 110k, the disk's tilt goes from π/2 to π, flipping to retrograde, while the jet is disrupted.During the RAD state, at t ≥ 110k, the disk flips again, which coincides with the flipping of the magnetic flux on the Northern hemisphere, in Fig. 1(b).After t ≳ 105k we show the tilt and precession of the counter jet, as it attempts to follow the rocking motion of the disk.a2).This shows that the evolution, stability, and destruction of the jets are correlated with the underlying state of the accretion flow near the launching region.
To understand the evolution of the system, we calculate the density-weighted specific angular momentum, ⃗ l = ⃗ r × ⃗ v = (l x , l y , l z ).In particular, we compute the magnitude of ⃗ l: The angular momentum computed this way is appropriate in flat spacetime, r ≳ 10R g .We normalize the resulting specific angular momentum with the specific Keplerian angular momentum l K around a BH of spin a (Shapiro & Teukolsky 1986): Figure 4(b) shows the magnitude of l/l K on a spacetime diagram.The three vertical black lines indicate the times of the 3D panels in Fig. 4(a).We choose l/l K ≃ 0.3 as a fiducial value to separate the accretion disk from the rest of the gas, which we overplot as a green contour in Fig. 4(b).Although initially the gas had zero angular momentum, within several R g from the BH the frame-dragging (Lense & Thirring 1918) due to the BH high spin (a = 0.94), combined with large-scale magnetic fields, to force the accretion flow to partially corotate with the BH, and the flow develops an azimuthal angular momentum component (see also R21).Additionally, the MRI excites disk turbulence and transports the angular momentum outwards, so the disk can grow in size.Indeed, by t ≃ 50k, the accretion disk size has grown up to 50R g .
Figure 4(a1) shows that once formed, the jet travels along the BH spin axis (perpendicular to the x − y disk plane, i.e., along the z− direction).Figure 4(c) shows the tilt and precession angles of both disk and one jet, at r = 50R g (dashed black line in Fig. 4b).We calculate the disk tilt T d (orange line) and precession P d (orange circles) angles using T d = cos −1 (l z /l) and P d = cos −1 [l y /(l 2 x + l 2 y ) 1/2 ], respectively.To compute the jet tilt and precession angles, we first find the Cartesian coordinates of the jet core centroid, at every radius r, using the magnetization-weighted average of x j , y j , z j , in regions where µ > 3 and b p ≥ b tor .The jet tilt is then given by T j = cos −1 (z j /r j ) (purple) and precession angle by P j = cos −1 [y j /(x 2 j + y 2 j ) 1/2 ] (purple crosses), where r 2 j = x 2 j + y 2 j + z 2 j .Figure 4(c) shows that during the MAD state (purple shaded) the tilt angle is small, P d ≤ π/16, at t ≲ 65k: in other words, the angular momentum vector is nearly parallel to the BH spin.For this reason, the disk precession angle, P d , is ill-defined, and we lower its opacity to minimize clutter.The jet tilt is, on average, similar to the disk tilt, with T j ≤ π/16.The jet precession angle shows highamplitude variability but on average clusters around P j ≃ π/2.This means that the jet on average has a preferred direction that is distinctly different than the rotational axis of the BH: indeed, Fig. 4 shows that both jets are skewed sideways, towards the yaxis.
At 50k ≲ t ≲ 65k, the radial extent of the disk shrinks down to ∼ 10R g , and by the end, it disappears, as the MAD state ends.This hardly seems like a coincidence, and we suggest that the system exits the MAD state when the accretion flow loses the rotational support.But what can cause the loss of angular momentum in the disk?It is possible that once launched, the jets perturb the infalling gas, and induce turbulence and vorticity, with angular momentum, on average, misaligned with the BH spin vector.Once the misaligned angular momentum reaches the disk, it can cancel out the disk's coherent angular momentum, which was relatively low, to begin with, due to the limited disk radial extent (R disk ≃ 50R g ≪ R B = 1000R g ).Moreover, strongly magnetized disks can launch winds that carry angular momentum outwards (Blandford & Payne 1982) and can further reduce the disk's angular momentum.Whatever the case, the disk of small size, ≲ 10R g , cannot support a jet with large-scale winds.The jet is left exposed to violent interactions with the ambient medium, for which the winds would normally act as a cushion.Additionally, the magnetic flux on the BH, once held by the disk, leaks out (Fig. 1b) and the power of the jet decreases (Fig. 1c).
From here on the system behavior changes drastically: at t ≳ 65k, our quasi-steady-state MAD transitions to a state dominated by the accretion of gas with a continuously varying direction of the angular momentum vector.For the duration, 69k ≲ t ≲ 95k, the disk tilt increases and the disk starts processing, while the angular momentum magnitude decreases.The flow with such low angular momentum is barely a disk, which we label as BAD state.The weakened and slower jet is now more affected by the presence of the ambient gas, and more specifically, the precessing and tilting disk start twisting the jet creating wiggles: this is similar to when one tags the end of a rope in a circular fashion creating circular waves that are transmitted along the rope.Around t = 78k, the jets transiently fall apart (due to the shutoff of jet power at the BH, Fig. 1c). Figure 4(b) shows that subsequently, a short-lived disk of size ∼ 30R g forms at 78k ≲ t ≲ 85k. Figure 1(c) shows that the jets remain active until t ≃ 95k.
At 95k ≤ t ≤ 110k, the disk gradually flips upside down and changes its sense of rotation from prograde (T d ∼ 0, rotating in the same sense as the BH) to retrograde (T d ∼ π, rotating in the opposite sense to BH). 4 This causes the jets to weaken, bend dramatically, and disrupt.We label this the rocking accretion disk state (RAD), where the accretion of gas with strongly misoriented angular momentum is capable of flipping the disk upside-down.We highlight a specific part of the RAD state, during 110k ≤ t ≤ 130k, where the disk gradually flips back to the prograde configuration, but stops at T d ≃ π/2.Interestingly, when looking at r = 15R g (not shown here) the tilt reaches all the way to T d ≃ π/4.The absolute dimensionless magnetic flux, ϕ BH , remains approximately constant during this transformation.In contrast, the magnetic flux through the northern hemisphere, ϕ N , smoothly changes sign in phase with the disk tilt angle: this is precisely the behavior we would expect if the entire disk-jet system underwent a flip.During this stage of the simulation, we witness a large angular momentum inflow (Fig. 4b), leading to an increase in the disk size to ∼ 200R g .The formation of an accretion disk, which grows in size, deprives the BH of the gas and suppresses the BH accretion rate.Although 15 times less powerful than in the MAD state, the weaker, intermittent and continuously reorienting jets fail to pierce through the infalling gas and are forced to bend around it, as seen in Fig. 4(a3), they efficiently couple their energy to the infalling gas and additionally reduce the BH accretion.This results in BH mass accretion rate suppression to Ṁ/ ṀB ∼ 1.7%, similar to the MAD state.

EXTERNAL VS INTERNAL KINK
Figure 5 shows a time series of jet snapshots, illustrating the propagation of the jets and the growth of the kink instability during the early time in our simulation, in the MAD state (31k ≤ t ≤ 65k): color shows the proper velocity, γv, projected along the y− and x−directions in left and right columns, respectively.We identify two early-time features, E1 and E2, that propagate along the jet and study how the kink instability grows in the jet.We pick these regions by either (i) visually searching for significant jet bends, or wiggles, or (ii) identifying jet regions with low values of the stability parameter, Λ ≲ 2 (e.g., in a movie of Fig. 3b).In both projections, the jet exhibits multiple wiggles along its length, that originate near the jet launching region and propagate out.Such wiggles can emerge due to the stochastic nature of BH accretion and wind-jet interactions (see, e.g., Fig. 4c).Both E1 and E2 features originate as small-scale wiggles and propagate out to r ≃ 1500R g .The jet remains mostly straight and maintains a relativistic velocity, ⟨γv⟩ ≃ 3.Such powerful jets as this one, in the MAD state, can accelerate to and sustain relativistic velocities (γv ≳ few): the relativistic motion stabilizes the jets against the kink instability via the relativistic time dilation (see eq. 18, which has the Lorentz factor is in the numerator).However, the jets propagate out to larger distances in a flat density distribution, they become more unstable (Sec.4.1): after reaching r ≃ 1500R g , the jets develop unstable regions of Λ ≲ 2 (Fig. 3), which implies that the kink instability growth timescale becomes comparable to the dynamical time of the jet, and the jets bends grow in amplitude.Indeed, both E1 and E2 features develop stronger jet bends that force the jet slows down, ⟨γv⟩ ≃ 1.7 (see eq. 22).These bends grow in amplitude and spatial scale, and culminate in essentially breaking up the jet into two segments, as seen in Fig. 5(d) for feature E1 and in Fig. 5(f) for feature E2.To calculate the propagation speed of E1 and E2 features, we fit straight lines through them, as we show in Fig. 5.Both features propagate along the jet at v/c ≃ 0.8 − 0.85, which is comparable but slightly lower than the jet velocity.
We have just discussed that in the MAD state the jets developed unstable features: jet bends that propagate outwards while slowly growing in amplitude.However, the jets managed to maintain their outward motion most of the way to the jet head.We also saw in Figs. 2 and 3 that the jets globally disintegrated shortly after the system left the MAD state at t ≃ 65k.Can we develop a criterion to predict when such catastrophic jet destruction is bound to happen?BT16 derived a simple analytic approximation for a global jet stability parameter that evaluates the ability of the kink mode to grow on the jet periphery leading to a global deformation of the jet body.For a jet of power, L j , that moves through an ambient medium of constant density, ρ a , at a Lorentz factor, γ = 1/ 1 − v 2 /c 2 , the stability parameter at the jet head distance, r = r h , is: where for simplicity we have omitted constant prefactors.Can eq. ( 26) predict the global jet disintegration we observe in our simulations?Upon exiting the MAD state, at t = 65k, the jet power drops four-fold (Fig. 1c).At the same time, the jet Lorentz factor decreases from ⟨γv⟩ ≃ 3 to ⟨γv⟩ ≃ 1.7, but v A remains roughly unchanged (compare Figs. 3b1 and 3b2).The distance to the jet head, r h does not change significantly, because the drop in power occurs over a short timescale (∆t ∼ 5k) and leaves the jet head little time to advance (compare Figs. 2b  and 2c).Thus, the stability parameter drops by an order unity factor Λ MAD /Λ tr ≃ 1.11, Why does the stability criterion, eq. ( 26), not capture the global jet instability we observe in the simulation? Figure 6 shows a time series of jet images at a later time, 68k ≤ t ≤ 70.5k, after the end of the MAD state and right up until the time when the jet becomes strongly kinked and globally unstable, a precursor to its global destruction.Figure 6(a)-(c) shows a relativistic jet (yellow), ⟨γv⟩ ≃ 3, which eventually slows down to transrelativistic speed (orange), ⟨γv⟩ ≃ 1.7, in Fig. 6(d)-(f).The dark purple regions are locations where ⟨γv⟩ ≲ 1 and where the jet is getting disrupted.We identify three bend features (L1, L2, and L3) that grow in time.The feature L1 is one of the last features to be ejected during the MAD state.In Figs. 6(a1,a2), the jet is mostly straight, apart from the large bend at r ≃ 2000R g : this is the bend associated with the L1 feature, and it grows as the feature approaches the jet head.The propagation speed of L1 is v L1 ≃ 0.75c, comparable to what we found for unstable regions at the earlier times, as seen in Figure 5. Early-time MAD jet exhibits multiple small-amplitude bends and a major bend located at r ≃ 1500Rg.A movie is available in Supplementary Information and on YouTube (link).The left and right columns show the jet's proper velocity, γv, projected along the y−, x− directions respectively.We fit a straight line through features along the jet, that show wiggles (E1, E2), and which we follow at different times to observe their growth.The slope of the white dashed line gives the propagation speed of the wobbles.For both E1 and E2 features, the velocity is approximately v/c ≃ 0.82 − 0.85.The wiggles, which can act as a seed for the kink instability, originate from the jet tilting and precessing about the BH spin.However, the fluid velocity is relativistic.γv ≃ 3, and the jet-fluid does not have enough time to become kink-unstable, not until it reaches the distance of r = 1500Rg, where it becomes globally-unstable.Fig. 5.The L2 feature was ejected after the system exited the MAD state after the jets weakened and became unstable, enabling the L2 bend to grow, as seen in Fig. 6(b2).Figure 6(a)-(d) shows that as the L2 feature propagates out and grows in amplitude, its propagation speed decreases from that comparable to L1 down to v L2 ≃ 0.6c.Ejected at an even later time, the L3 feature starts out as a bend of high amplitude already at r ≲ 500R g , in Fig. 6(c).Its propagation speed, v L3 ≃ 0.35c, is the lowest among all E and L features and is, therefore, the most unstable.Indeed, Fig. 6 shows that the L3 feature grows to become a 90-degree jet bend, at r ≃ 800R g , that eventually globally disrupts the jet.
Summing up, we are witnessing a runaway growth of the kink instability that generates large-scale, strong jet bends.This implies that eq. ( 26), for some reason does not apply to our weakened jets in the BAD state (red shaded region, 69k ≤ t ≤ 95k, in Fig. 1).Namely, once the system enters this state, the linear stability analysis (eq.26) under-predicts the growth of the global kink.What is missing from the linear stability criterion, eq. ( 26), are the non-linear effects that combine to cause a catastrophic run-away of the instability.Namely, the fluctuations in the disk tilt (Fig. 4c) cause the jets to wobble around the BH axis and propagate in directions offset from the jet pre-drilled path.This causes the jets to plow up the gas along new paths, through the shocked cocoon material.Be-Figure 6. Post-MAD jet gets progressively disrupted, with small bends fully developing into large bends.Similar to Fig. 5, the color shows the proper velocity, γv, and the two columns are the projections along y− and x−direction respectively.This jet has an average power of 4 times smaller than the average power during the MAD state.The drop in power, is also associated with the tilt of the disk slowly increasing, as shown in Fig. 4(c).The power drop leads to jets that cannot efficiently accelerate to previously relativistic velocities, and as the jets follow the tilt of the disk, they run into the ambient medium, shock, and decelerate even more (see Fig. 3b2).External kink starts affecting the jet bends.As the bends get progressively larger, the fluid can become internally kink-unstable leading to the "breaking" of the jet.We fit a straight line through regions that develop large bends, L1, L2, and L3, and we estimate the wiggle propagation velocity from the slope.The L1 spot follows a bend developed during the end of the MAD state when the fluid was mostly stable to the kink.At r ≃ 2000Rg, a bend develops and the wiggle grows without significantly bending the jet until the jet fluid smashes at the head.The estimated propagation speed is slightly smaller than the MAD equivalent, at v ≃ 0.75c.The L2 spot propagates slower than L1, at v ≃ 0.6c and the wiggle amplitude is noticeably larger.Finally, L3 is launched well into the weaker jet phase, and the jet bend becomes non-linear until the full jet disruption at r ≃ 800Rg.The L3 wiggle propagates at v ≃ 0.3c, significantly slower than either L1 or L2, due to the fact that the large amplitude bends launched near the BH interacts with more ambient gas which slows it down even more.
cause the cocoon was generated at earlier times when the jets were four times more powerful, the weakened jets have difficulty displacing the material, and get deflected: this causes the jets to bend and slow down (eq. 22 and Fig. 3a1,b1).Slower velocity accelerates the kink (and bend) growth and further slows down jet propagation.As a result, such interactions lead to a runaway slowdown not only of the jet fluid but also of jet bends (Fig. 6).In fact, the global kink instability we are witnessing here grows in the frame of the bends, which move slower than the jet fluid and hence more susceptible to the kink instability.This is an example of the non-linear development of an external kink instability, which deforms the jet interface with the ambient medium.This is different than the internal kink instability, which works on disrupting the jet core but leaves the jet-ambient medium interface intact.The two flavors of the kink instability can coexist, i.e., the development of the external kink instability can trigger the internal kink instability.
To more systematically study the nonlinear development of the global, external kink instability, Fig. 7 presents a spacetime diagram, where the color shows the logarithm of the dimensionless jet curvature, where κ is jet curvature and ⃗ s is a unit vector parallel to the jet.High values of the dimensionless curvature (yellow) indicate strong bends, whereas low values (black) correspond to mostly straight jets.We see that even at small distances, r ≲ 1000R g , the jet develops multiple bends along its length, due to the interactions with the ambient medium and disk winds.We plot the tracks of all features we identified in Figs.5-6: the slopes of the tracks give the bend propagation speeds.During the MAD state, t ≲ 65k, the jet on average is mostly straight, and the bends travel at nearly the same speed as the jet fluid.If we follow the E1 and E2 features, they develop significant (orange) bends, at r ≃ 1500R g , where dimensionless curvature reaches log 10 ( κ) ≳ 1.5.
Further away, at r ≳ 2000R g , the jets get even more twisted: they become globally unstable, break apart, and energize the cocoon.Inflated by the jets, the cocoon expands at a nonrelativistic speed of v ≃ 0.05c, as revealed by the slope of the jet-cocoon boundary, which is seen as the orange-black transition on the right of Fig. 7.In fact, Figure 7 reveals that the entire bright orange region of the strongest jet bends shifts in time towards larger radii.This suggests that the proximity to the cocoon exacerbates jet bends: cocoon convective motions can displace and bend the jets sideways.Indeed, one can visually identify such convective cocoon motions in the movie of Fig. 2. Consistent with this picture, once the cocoon expands out to larger radii, the jets manage to avoid significant bends out to larger radii as well: the radius beyond which the jets start to show significant bends slowly increases from r ≃ 1500R g at t = 50k to r ≃ 3000R g at t ≃ 68k, just before the MAD state ends and the jets become globally unstable.
Visually following the slope of the "late-time" track of feature L1 in Fig. 7, we estimate the propagation speed of jet bends prior to disruption, v L1 ≃ 0.8c, which is consistent with the best-fit values in Fig. 5.However, once the MAD state comes to an end, at t ≃ 65k, the jet power drops.This causes each subsequent feature to travel at a slower speed, as indicated by the steepening of the slopes, decreasing to v L2 ≃ 0.6c and eventually to v L3 ≃ 0.35c.This is consistent with the feature propagation speeds we measured in Fig. 6.At t ≳ 70k, the jet develops extreme bends and shortens to r ≲ 2000R g , before globally disrupting at t ≃ 78k.

DISCUSSION AND CONCLUSIONS
We study large-scale jet propagation and survival in global 3D GRMHD simulations of weakly magnetized, zero-angular momentum gas accretion onto rapidly-spinning BHs, a = 0.94.The gas accretes from the Bondi-to-event horizon radius, traversing the largest scale separation to date, R B = 10 3 R g , in a single 3D GRMHD simulation.Over time, the gas drags in from the ambient medium large-scale vertical magnetic flux, which readily accumulates on the BH and launches powerful relativistic collimated jets.
Although initially our accretion flow started out without any angular momentum, as the accretion flow goes MAD, it forms an accretion disk of size ∼ 50R g (see also R21; K23).Thus, no pre-existing accretion disk is required to form jets and collimate them via winds, as long as the BH is rapidly spinning and has accumulated dynamically-important large-scale poloidal magnetic flux on its event horizon.Interestingly, the MAD state achieved in our simulation is long-lived: its duration of ∆t = 34k is 3 times as long as the total simulation duration in R21 and about half as long as the simulation duration in K23, where in both cases the system only transiently enters the MAD state.Although our MAD state lasts comparatively longer, it still has a finite duration: the BH eventually loses the large-scale magnetic flux, and the jets fall apart.The MAD state ends in a fashion similar to that observed by R21 and K23.Once our jets transiently re-appear, they do not reach beyond the Bondi radius.
To explore the origin of these differences, we repeated our simulation at a reduced scale separation, R B /R g = 10 2 , to match that of R21, see Appendix C. We carried out the simulation for t ≃ 110k, a factor of ∼ 5 longer than R21.The system turns MAD (ϕ BH ≳ 50) at 70k ≲ t ≲ 100k, which suggests that the MAD state can establish itself, but it might take longer simulation durations in order to witness it.Why do MADs in the simulations of K23 survive for a rather short duration?One possibility is that our simulations contain more magnetic flux closer to the BH: although both simulations start with β = 100 gas, we use a constant vertical field, B z = constant, which results in the enclosed poloidal magnetic flux scaling as Φ ∝ B z r 2 ∝ r 2 .In contrast, K23 use a parabolic-like field, B z ∝ r −5/4 , whose flux has a flatter radial profile, Φ ∝ B z r 2 ∝ r 3/4 .Thus, out to the same distance, our simulations contain a larger flux reservoir.It is possible that this allows our simulations to transport the magnetic flux to the BH at a higher rate than the BH loses it due to the magnetic flux eruptions.Additionally, they have an adiabatic index Γ = 4/3 describing radiation-dominated gas.Because this makes the gas more compressible, this can reduce the disk scale height and make it harder for the disk to hold on to the magnetic flux on the BH.
Apart from the early-time transient, our MAD jets carry the maximum power and extract BH rotational energy at high efficiency, η ≃ 190%.The jets reach distances of r ≳ 4000R g , Figure 7. Space-time diagram shows that normalized jet-curvature, κ, grows due to the kink instability as it advects along the jet.We overplot the tracks of jet features in the MAD state (E1, E2) and BAD state (L1, L2, L3), also seen in Figs. 5, 6.The slopes of these tracks give us the propagation speeds of the features.During the MAD state, the large-scale bends tend to develop at r ≃ 1500Rg, and propagate to larger scales at roughly the speed of the jet fluid, v ∼ 0.8c, as the jet powers the cocoon.The cocoon outer boundary (right-most extent of the orange region boundary) expands at v ≃ 0.05c, i.e., much slower than the jet features.After the flow exits the MAD state (t = 65k) and enters the BAD state (t = 69k), the jet power drops, and strong jet bends form at much smaller distances, r/Rg ≃ 500 − 1000, and propagate at much slower speeds, v ∼ 0.3c (e.g., for feature L3).Eventually, the jet is globally destroyed, at t ≃ 78k, as its power vanishes (Fig. 1c).
or 3.5 orders of magnitude in scale separation.In the inner ∼ 50% of the jet length, r ≲ 1500R g , the jets remain mostly bend-free.However, further out, the jets develop twists and bends, as they become kink-unstable.Magnetic energy dissipates into thermal energy, and the jet magnetization drops until the magnetic fields reach equipartition with the gas thermal energy.Multiple locations along the jets show signs of recollimation and/or oblique shocks, associated with small-scale bends of the jets.To interpret this association, recall that at recollimation point the fluid has to decelerate, which can lead a marginally stable jet to become kink-unstable, and twist as a result.Figure 3(b) shows that eventually the instability grows non-linear and breaks up the jet into segments.Because each segment is no longer constrained length-wise by the formerly adjacent jet segment, it is free to expand longitudinally and inject part of its energy into the cocoon, resulting in a dramatic drop in jet energy flux at the beginning of each segment.
The fact that the jets can reach r ≃ 1500R g without getting globally disrupted is not trivial.R21 find that their jets become kink unstable at r ≃ 200 − 400R g .Using eq. ( 26) to extrapolate to larger Bondi radii and assuming a density profile of ρ ∝ r −1 within the Bondi radius, they predict that for R B ≳ 800R g , jets should become kink-unstable.Because observations show many jets survive much greater distances, this would suggest that such jets cannot be powered by low-angular momentum accretion.Surprisingly, our jets make it out to r ≃ 1500R g without significant global bends.One possible reason for this difference is that when our jets first form, they navigate a steeper density profile, ρ ∝ r −3/2 , than the one established at late times due to feedback processes, ρ ∝ r −1 .The steeper profile makes it easier for the jets to escape unscathed (BT16).Another difference is that soon after jet launching, the accretion flow goes MAD and maximizes the jet power.It is possible that, once our system exits the MAD state, and the jets transiently shut off, they would struggle to make it out to the same distances even if they went MAD again because the feedback causes the ambient density profile to flatten, ρ ∝ r −1 , at late times.
At 50k ≲ t ≲ 65k, the MAD disk shrinks in size down to ∼ 10R g , the MAD state terminates (ϕ BH < 50), and the mass accretion rate increases.Initially, the jet power stays roughly constant, before dropping.Without the disk wind to shield the jets, the infalling gas slams into and bends them.The bends propagate out along the jets in a spiral pattern.Once the jet power drops, the mass accretion rate increases even further, and the system enters the BAD state at 69k ≲ t ≲ 95k.During this state, the jet becomes even more susceptible to bends, which are amplified by the external kink instability: in fact, the bends grow non-linear, drag against the cocoon material, and slow down.This leads to a run-away growth of the bends that twists the jets into a helix.Thus, the low-angular momentum accretion system self-consistently forms precessing helical jets.We offer this as an observational alternative of periodical and curved jets, instead of binary companions or BH-torquing of a misaligned disk (Cui et al. 2023), especially since wobbles and precession can be hard to distinguish.
As the system went MAD, the jets reached the maximum energy efficiency and suppressed the accretion rate by a factor of ∼ 70: only a tiny fraction, Ṁ/ ṀB ≃ 1.5%, of the Bondi accretion rate reaches the BH, and the feedback (e.g., by BHpowered jets and disk-powered winds) ejects the remaining 98.5% of the gas.L22 used an identical setup, except for including non-zero angular momentum in the ambient medium (with circularization radius, R circ = 30R g ), and found similar levels of Ṁ suppression.This suggests that most of the feedback is done by the powerful jets as opposed to disk winds.Indeed, Fig. 1(a) shows that Ṁ remains roughly constant over the entire duration of the BAD state (69k ≤ t ≤ 95k), despite the formation of an accretion disk at 78k ≲ t ≲ 86k.In fact, we see that Ṁ transiently increases whenever the jet power vanishes, e.g., at t ≃ 70k and 95k ≲ t ≲ 105k.This indicates that it is the jets, especially in the MAD state, that dominate the Ṁ suppression, as opposed to the rotationally supported disk and its wind.At t = 71k, the power of the jet shuts off, and the jet gets globally destroyed.The subsequent jet formation has weak jets that mainly go through the pre-drilled path from the prior large-scale jets, with an average efficiency of η ≃ 15%, until t = 94k.For the rest of the simulation, t ≳ 95k, the system enters the RAD state, where the accreted material brings in gas with misaligned angular momentum, with respect to the BH spin.Newly formed intermittent and weaker jets don't surpass the Bondi radius intact, but promptly get disrupted, leaving behind weakly-magnetized buoyant bubbles.When the jets are not active, the mass accretion rate can reach all the way up to Ṁ/ ṀB ≃ 10%.Although the system cannot revert back to an apparent steady state, there are periods where the intermittent, and weaker, jets along with the formation of an ephemeral misaligned accretion disk, suppress the mass accretion rate down to values in the MAD state.For example, during 110k ≲ t ≲ 130k, the mass accretion rate averages to Ṁ/ ṀB = 1.7%, while the jet power is weaker compared to the MAD values by ∼ 15.The outflow efficiency averages at η = 11%, while the normalized magnetic flux at ϕ BH = 18.
During the MAD state, the powerful and relativistic jets (γ ≃ 4) showcase extensive lobes and hotspots similar to FRIIs.During the RAD state, the FRII jets have been disrupted at their base, and the new jets are intermittent, and propagate more slowly (γ ≃ 2) in varying directions without ever making it out of the Bondi radius (r ≲ 800R g ).This is supported by observations, with FR0 jets appearing mildly relativistic and with small sizes (r ≲ 1kpc) making them hard to observe.For the first time, in a GRMHD simulation, we are witnessing a failed FRII jet transitioning into an FR0 state.Thus, we think that FR0 should be abundant in environments with extremely low angular momentum.Even though such systems are capable of creating FRI or FRII jets, if the disk loses its angular momentum and shrinks, the system inevitably turns RAD.Additionally, the jets that are globally destroyed during the transition from MAD to BAD, might be the source of γ-rays observed by a handful of FR0s (Paliya 2021), as a significant amount of the jet energy is dissipated into heat.The jets in the BAD and RAD states might be significant contributors to the isotropic γ−ray background, the diffuse neutrino background (detectable by IceCube), and the ultra-high-energy cosmic rays.We will explore the radiation signatures of a globally kink-unstable jet, as well as the FR0 jets in a follow-up paper.

B. ADAPTIVE MESH REFINEMENT CRITERION
To identify the jet material, we use a cutoff in a proxy for entropy, S = u g /ρ Γ , where u g is the internal energy density, ρ is comoving frame density of the gas, and Γ = 5/3 is the gas polytropic index for a monatomic non-relativistic gas.The reason is that magnetic dissipation in the jets increases the value of S, making it appear as "hot".For our run, the cut-off is chosen at S ≥ 10.We additionally refine the cocoons, created when the jets mix with the ambient gas resulting in a shocked and mildly-magnetized gas.The cut-off for the cocoons is chosen at S ≥ 0.1.We also derefine blocks where the jets and/or cocoons left, when the entropy values dropped to [50% − 100%] of the cut-off.The AMR criterion is activated once the half-opening angle of either jets or cocoons becomes narrower than 48 cells (see Gottlieb et al. 2022b), which is the width of a single AMR block, whose resolution is N B r × N B θ × N B φ = 56 × 48 × 48.The maximum effective resolution in the jets can reach N eff r × N eff θ × N eff φ = 3, 584 × 768 × 1, 536.This ensures that both the jets and their cocoons are well-resolved as they propagate out to large radii: jet opening angles of 0.04 rad = 2.3 • are resolved with ≥ 10 cells out to distances of 5, 000R g .

C. REDUCED SCALE SEPARATION
We use the same physical parameters, as in the fiducial run, except we reduce the Bondi-to-event horizon scale separation by 1 order of magnitude, R B /R g = 10 2 .The grid setup is mostly identical to our fiducial run (see Sec. 2).In this case, we did not tilt the BH spin 90 • -degrees with respect to the polar axis.The base-grid resolution is the same: N r × N θ × N φ = 448 × 96 × 192 cells in the r-, θ-, and φ-directions.We activate 1 level of static mesh refinement (SMR) at r ≥ 6.5R g , doubling the resolution in each dimension for an effective resolution 896 × 192 × 384 cells.Figure 8 shows the normalized flux, ϕ BH , on the BH.The system turns MAD at 70k ≲ t ≲ 100k.
Figure 1(b) shows the time dependence of normalized absolute magnetic flux,

Figure 2 .
Figure2.Three-dimensional volume renderings of density at different times show that as jets propagate out of the Bondi sphere, they become kink unstable and globally fall apart.The BH spin is pointing in the direction of jet propagation to the right in panel (a).A movie is available in Supplementary Information and on YouTube (link).[panel a]: The highly powered MAD jets (red) easily escape the Bondi sphere without obvious signs of instability, except near the head.At the contact points (hotspots), at which the jets drill through the ambient gas, the shocked jet material spills sideways, creates backflows, and forms the cocoons (yellow).The relativistic motion of the jets and cocoons shocks the ambient gas via bow shocks (blue).[panel b]: As the jets propagate further out, this leaves enough time for the kink instability to grow, and it leads to helical bends near the hotspot.[panel c]: The jet power decreases by ∼ 4, and the kink unstable regions move in closer to the BH.The decrease in power propagates through the jets faster than the cocoons can react, and the weakened jets get squeezed by the cocoons.[panel d]: The central engine shuts off, and the jets break apart globally.[panel e]: Weakened jets wobble and lose focus: without consistent energy supply from the jets, mildly-magnetized cocoons become relics, which buoyantly rise outwards.Meanwhile, the wobbling jets transiently form and fall apart, without an opportunity to reach the Bondi sphere.[panel f]: Intermittent and weaker jets form, but get deflected sideways, as the accretion disk tilt continuously changes.The inset image shows that the jets develop extreme bends inside the Bondi radius.

Figure 3 .
Figure 3. Kink instability disrupts the jets as their power drops four-fold upon the exit out of the MAD state at t ≈ 65k.[panels (a)]: Transverse projection of σ-weighted proper velocity, γv, reveals large-scale jet bends, telltale signs of the kink instability (brown shows low and blue high values, see color bar).[panels(b)]: Whereas at the earlier time (t = 58k, panels a1 and b1) only the outer jet, r ≳ 2000Rg, shows bends and satisfies the kink instability criterion, Λ ≲ 2, at the later time (t = 71k, panel a2 and b2) most of the jet, r ≳ 800Rg, shows bends and satisfies the criterion.Jet bends cause the core-average proper jet velocity γv (magenta) to significantly decrease from γv ≃ 3 (panel b1) to γv ≲ 2 (panel b2) and result in a comparable decrease in Λ (eq.20).Although the jet radius, Rj (dark red), also decreases, this noticeably affects neither the core radius, Rc (gold), nor core-average Alfvén velocity, vA (green).[panels (c)]: Angle-integrated fluxes, or power, through the jet (µ > 3) normalized to the time-average total MAD power, ⟨ ĖMAD⟩.The normalized total jet power, ĖTOT (black), remains approximately constant in stable jet regions.The kinetic energy component, ĖKE (light blue), gradually increases at the expense of the EM component, ĖEM (purple), which dominates the total energy power until the instability sets in at r ≳ 2000Rg in panel (c1) and r ≳ 800Rg in panel (c2).Panels (a1) and (a2) show that in the unstable regions, the jet develops one or more bends.Whereas just before each substantial bend the thermal component (orange-red) grows up to equipartition with the EM component (purple), ĖTH ∼ ĖEM, after each bend all energy power components drop, followed by a gradual increase until the next bend (e.g., for bends at r/Rg ≃ 800, 1300, 2000, 2500, and 3000 in panel c2): we attribute this to the "accordion" effect -the jet bends allow the jets to elongate and rarefy longitudinally.The rest-mass power ĖRM (orange) remains subdominant to the rest of the components.[panels (d)]: Jet-average proper velocity, γv (light blue, obtained via γ ≡ ĖKE/ ĖRM), increases as the jets accelerate at the expense of the decreasing jet-average magnetization, σ ≡ ĖEM/ ĖKE (purple), before decelerating at jet bends.Jetaverage enthalpy, h ≡ ĖTH/ ĖKE (orange-red), increases at jet bends due to energy dissipation.The µ ≡ ĖTOT/ ĖRM parameter (black) decreases due to the jet losing its energy to the ambient medium (via ambient gas displacements and shocks).

Figure 4 .
Figure 4.The morphology, strength, and directions of the angular momentum of the accretion flow during the MAD, BAD, and RAD states.[panels(a)]: We show 3 distinct 3D volume-rendered images of density.Panel(a1) captures the system during the MAD state, at t = 50K, where the jets (purple) follow the direction of the BH's spin, and the accretion disk (orange) lies perpendicularly to the jets.Panel(a2) captures the system in the BAD state, at t = 71k, with weakened jets showing strong bends near the BH and getting significantly kinked.Panel(a3) captures the system in the RAD state, at t = 135k.The jets are on average 1 order of magnitude weaker and misaligned from the spin direction, The disk is also significantly tilted, at about π/2, and processing, causing significant bends and wobbles in the jets.[panel (b)]: We plot the density-weighted specific angular momentum normalized to the local Keplerian value in a space-time diagram.Angular momentum builds around the BH promptly after freefall (tff ≃ 22k) with an average disk of size 10 − 50Rg, chosen for l ≳ 0.3lK.Between, 65k ≤ t ≤ 100k the accreted angular momentum forms intermittent and tilted disks with l ≤ 0.3lK which results in intermittent and wobbly jets.After t ≃ 100k angular momentum starts building up, this time reaching larger radii ≃ 200Rg, however, near the BH the intermittent jets push most of the gas away while oscillating wildly, stopping the angular momentum from coherently adding up and stabilizing the jets.[panel (c)]: We show the tilt and precession angles of the disk and north jet, measured at a distance of r = 50Rg.During the MAD state (purple shaded), both the tilt of the disk Pd ≲ π/16 (orange) and jet Pj ≲ π/16 (purple), as they are aligned with the BH spin.The disk precession angle (orange circles), Td, is illdefined as the disk tilts about the z-axis, while the jet precession angle (purple crosses), Td ≃ π/2, which means it propagates at a preferential direction.Once we transition from the MAD to the BAD state (red shaded), at t ≥ 65k, the disk angular momentum has decreased, its tilt increases, reaching up to π/2, and rapidly precesses.The jet tries to follow and forms large bends (see Fig.4a2) that lead to its destruction.Between 95k ≲ t ≲ 110k, the disk's tilt goes from π/2 to π, flipping to retrograde, while the jet is disrupted.During the RAD state, at t ≥ 110k, the disk flips again, which coincides with the flipping of the magnetic flux on the Northern hemisphere, in Fig.1(b).After t ≳ 105k we show the tilt and precession of the counter jet, as it attempts to follow the rocking motion of the disk.