Bridging Bondi and Event Horizon Scales: 3D GRMHD Simulations Reveal X-Shaped Radio Galaxy Morphology

X-shaped radio galaxies (XRGs) produce misaligned X-shaped jet pairs and make up $\lesssim10\%$ of radio galaxies. XRGs are thought to emerge in galaxies featuring a binary supermassive black hole ($\rm SMBH$), $\rm SMBH$ merger, or large-scale ambient medium asymmetry. We demonstrate that XRG morphology can naturally form without such special, preexisting conditions. Our 3D general-relativistic magnetohydrodynamic (GRMHD) simulation for the first time follows magnetized rotating gas from outside the $\rm SMBH$ sphere of influence of radius $R_{\rm B}$ to the $\rm SMBH$ of gravitational radius $R_{\rm g}$, at the largest scale separation $R_{\rm B}/R_{\rm g} = 10^3$ to date. Initially, our axisymmetric system of constant-density hot gas contains weak vertical magnetic field and rotates in an equatorial plane of a rapidly spinning $\rm SMBH$. We seed the gas with small-scale $2\%$-level pressure perturbations. Infalling gas forms an accretion disk, and the $\rm SMBH$ launches relativistically-magnetized collimated jets reaching well outside $R_{\rm B}$. Under the pressure of the infalling gas, the jets intermittently turn on and off, erratically wobble, and inflate pairs of cavities in different directions, resembling an X-shaped jet morphology. Synthetic X-ray images reveal multiple pairs of jet-powered shocks and cavities. Large-scale magnetic flux accumulates on the $\rm SMBH$, becomes dynamically important, and leads to a magnetically arrested disk state. The $\rm SMBH$ accretes at $2\%$ of the Bondi rate ($\dot{M}\simeq2.4\times10^{-3}M_{\odot}\,{\rm yr}^{-1}$ for M87*), and launches twin jets at $\eta=150\%$ efficiency. These jets are powerful enough ($P_{\rm jets}\simeq2\times10^{44}\,{\rm erg\,s}^{-1}$) to escape along the spin axis and end the short-lived jets state whose transient nature can account for the rarity of XRGs


INTRODUCTION
Modeling how the hot interstellar medium (ISM) or intercluster medium (ICM) reaches from the outskirts of a galaxy or cluster of galaxies down to the central super-massive black hole (SMBH) remains an unsolved, grand challenge problem. Observationally, it is clear that some poorly constrained fraction of the gas makes it down all the way to the SMBH, forms an accretion disk, and powers an active galactic nucleus (AGN). Such systems can produce copious radiation and mechanical outflows, broad winds and tightly collimated relativistic jets. The jets can propagate through the galaxy, generate shocks, inflate bubbles, and heat up the the surrounding medium (Zhuravleva et al. 2016;Li et al. 2020;Martizzi et al. 2019). It's commonly accepted that jet launching near the event horizon is powered by extraction of black hole (BH) aretaioslalakos2022@u.northwestern.edu rotational energy (Blandford & Znajek 1977). The jets initially align with the BH spin axis (McKinney et al. 2013), but interactions with the disk outflows can change their direction and they follow the angular momentum of the accretion disk (Liska et al. 2018). Thus, the jet direction is determined by the complex interplay between the properties of the BH and its feeding.
A small but significant, ∼ 5 − 10%, fraction of observed radio galaxies displays an intriguing X-shaped morphology, suggesting the presence of two pairs of jets at different angles (Leahy & Williams 1984;Yang et al. 2019). Such X-shaped radio galaxies (XRGs) have been hypothesized to emerge due to: (i) restarted activity of the central BH producing jets with different orientation (Bruni et al. 2019(Bruni et al. , 2020, (ii) jet deflection by the hot and oblique halo/ISM of the galaxy (Cotton et al. 2020), (iii) the presence of a SMBH binary, each component producing its own pair of jets (Lal & Rao 2007), (iv) reorientation of the SMBH spin axis due to SMBH coalescence (Merritt & Ekers 2002;Roberts et al. 2015), and (v) precession of the jet axis (Parma et al. 1985) due to e.g. Lense-Thirring precession of a tilted accretion disk (Liska et al. 2018). These might not be the only scenarios for forming the X-shaped morphology, and combinations of these are also plausible.
The goal of this work is to investigate whether X-shaped jet morphology can emerge spontaneously, without the above special, preexisting conditions. To reach the SMBH, the gas needs to make its way from the ISM/ICM, located outside the sphere of influence of the SMBH, or the Bondi radius R B = GM BH /c 2 ∞ ∼ 20 pc × (M BH /10 9 M ), at which the thermal pressure is high enough to counter the free-fall of the gas (Bondi 1952), down to the scale of the 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 sound speed in the ISM/ICM. Galaxy simulations have been able to follow the gas flow down to ∼ 1 kpc (e.g., Anglés-Alcázar et al. 2015 and even sub-pc scales in hyper-refined Lagrangian simulations (Anglés-Alcázar et al. 2021), approaching and even reaching inside of the Bondi sphere. However, the event horizon is still orders of magnitude below the parsec scale, demanding dedicated numerical effort to bridge this gap in scale separation.
The simplest description of accretion is the sphericallysymmetric analytic Bondi (1952) model, which does not include any rotation. Bondi accretion with nonzero angular momentum has been simulated by multiple groups (e.g., Proga & Begelman 2003a,b;Cunningham et al. 2012;Li et al. 2013;Suková & Janiuk 2015;Suková et al. 2017;Xu & Stone 2019;Kaaz et al. 2019;Palit et al. 2019;Waters et al. 2020), often utilizing axisymmetric, non-relativistic, hydrodynamic simulations. Accretion from a realistic ISM will be necessarily magnetized and have at least some angular momentum support. Thus, it is important to model both rotation and magnetic fields. Rotation breaks the symmetry of the problem, provides rotational support, and tends to enhance mass outflows, reducing the mass reaching the BH. Magnetic fields are important, as they form magnetically-powered outflows that can inject energy into the ambient gas. It is also crucial to include general relativistic (GR) effects, in order to accurately represent the energy and momentum feedback by the central SMBH via launching BH-powered jets and disk-powered winds and properly account for the poorly understood effects the BH has on the inner boundary condition of the flow. It is also critical to extend the studies to 3D, to model the essentially non-axisymmetric magnetized turbulence, the associated angular momentum transport in the accretion disk, and to account for the development of 3D magnetic kink instabilities in the jets (see also Ressler et al. 2021;Kaaz et al. 2022;Jia et al. 2022). Tchekhovskoy & Bromberg (2016) initiated relativistic jets by the magnetized rotation at the inner grid boundary of radius R B = 0.1 kpc and followed their propagation through the ISM over distances of tens of kpc for long durations, suf-ficient for the instabilities to develop. They found that the jets inflated cavities, morphologically similar to M87* and other low-luminosity AGN as seen, e.g., by Chandra (Forman et al. 2017). Barniol Duran et al. (2017 showed that the jets can undergo internal 3D kink instabilities and dissipation triggered by the change in the radial density profile at the Bondi radius (Russell et al. 2015).
Incorporating all of these effects, 3D GRMHD simulations are unique tools enabling self-consistent studies of gas accretion from the Bondi scale, including the formation of a turbulent accretion disk, the launching of relativistic jets, and the development of 3D magnetic kink instabilities that can dramatically affect the jet morphology and the state of the ambient gas. Despite the success of GRMHD simulations, direct modeling of 3D accretion spanning the full range of log 10 (R B /R g ) 5−6 orders of magnitude in distance (and 8−9 orders of magnitude in time) between the BH and Bondi scales still remains computationally prohibitive. Our approach is to reduce the scale separation to the maximum that is computationally feasible. While the scale separation is smaller than in reality, the more orders of magnitude it spans, the closer the flows can approach a self-similar regime that approximates nature. By varying the degree of scale separation, we can evaluate the effects our particular choices for the R B /R g ratio have on the structure of the flow, and analytically extend the results to the values of R B /R g ratio found in nature.
In this Letter, we present the results of the first 3D GRMHD simulation that follows the accretion of an initially magnetizedrotating gas from a separation scale up to R B /R g = 10 3 , which in 3D GRMHD is the largest to date and close to the maximum possible value attainable with current computational resources when evolved over astrophysically-interesting times. Starting with the simplest, axisymmetric state of constant density, rotational profile of the gas, and vertical magnetic fields outside R B , we evolve the system until the formation of an accretion disk, a pair of relativistic jets, and their interaction with the ambient gas. In Sec. 2 we present the numerical setup and the physical parameters of our model, and in Sec. 3 we present the results: the emergence of the X-shaped morphology (Sec. 3.1), the jet shape (Sec. 3.2), and synthetic X-ray images (Sec. 3.3). In Sec. 4 we conclude. We use units of G = M = c = 1.

SETUP
We carry out a 3D GRMHD simulation using our GPUaccelerated GRMHD code H-AMR , see Porth et al. 2019 for comparisons with other current GRMHD codes). We employ spherical polar coordinates, r, θ, ϕ, and choose a spherical polar grid that is uniform in log r, θ, and ϕ variables. The grid spans the range of (0.97R g , 10 5 R g ) × (0, π) × (0, 2π) and has the resolution of N r × N θ × N ϕ = 192 × 256 × 128 cells in the r-, θ-, and ϕ-directions, respectively. We choose it such that the first 5 radial cells are inside the event horizon: this en-sures that the BH exterior is causally disconnected from the inner radial grid boundary. We use outflow boundary conditions in the radial, transmissive conditions in the polar, and periodic conditions in the azimuthal directions. We adopt axisymmetric initial conditions described by the following physical parameters: (i) Density profile: Outside the Bondi radius, we set the ambient density to a uniform value of ρ ∞ = 1. Since the simulation includes neither the cooling (radiation) nor self-gravity effects, the results are scale-free in density: the simulation results can be freely rescaled for any value of ρ ∞ . Within the Bondi radius, we place an empty cavity. The reason for this is we do not want to impose an initial density profile that may influence the subsequent steady-state solution.
(ii) Bondi radius: As we discuss in Sec. 1, in order to make the simulations affordable computationally, we choose a smaller R B than inferred from observations (R B ∼ 10 5−6 R g ). Namely, we adopt a value of R B = 10 3 R g : this is the largest scale separation between R B and R g that has ever been achieved in 3D GRMHD that we are aware of, yet it still allows us to evolve the system for astrophysically-interesting times. When discussing large-scale jet propagation outside R B , we convert the simulated lengths and times to physical lengths and times by associating the Bondi scale and its light crossing time in simulation units, R B and R B /c = 10 3 R g /c, with those in the physical units for M87*, R B = 0.15 kpc and R B /c = 0.15 kpc/c = 475 yr, respectively (Russell et al. 2015). Our choice of R B /R g naturally sets the sound speed in the ambient gas, c s = (R B /R g ) −1/2 c = 10 −3/2 c ≈ 0.03c. We adopt a non-relativistic ideal gas equation of state, p g = (γ − 1)u g , where γ = 5/3 is the polytropic index for a monatomic gas, and u g and p g are the gas internal energy density and pressure as measured in the comoving frame of the gas.
(iii) Circularization radius: We assign a rotational profile to the gas, such that on each spherical shell of radius r, the gas undergoes solid-body rotation at the angular velocity, ω = l 0 /r 2 , around the z-axis. Here, l 0 characterizes the specific angular momentum of the gas in the equatorial plane. We choose l 0 = constant, such that gas specific angular momentum, l l 0 sin 2 θ, reaches its maximum value at the equatorial plane and smoothly drops down to zero near the poles (see also Palit et al. 2019). 1 The value of l 0 depends on the circularization radius, R circ , at which the angular momentum in the equatorial plane equals the local Keplerian value. We adopt which divides approximately equally the scale separation between R g and R B and results in sub-Keplerian rotation at the Bondi radius, l 0 /l K (R B ) = (R circ /R B ) 1/2 ≈ 0.17.
(iv) Initial gas magnetization: We adopt the initial magnetic field threading the gas to be asymptotically a homogeneous vertical magnetic field in the z-direction and initialize it by setting the covariant magnetic vector potential, which has only one non-zero component: A ϕ = (r 2 − R 2 B ) sin 2 θ at r ≥ R B and 0 otherwise. This ensures that A ϕ and the magnetic field vanishes inside R B and does not bias the formation and evolution of the accretion flow there. Physically, the covariant ϕ-component of the vector potential, A ϕ , gives the poloidal (pointing in the rand θ-directions) magnetic flux enclosed by an axisymmetric ring, (r, θ), divided by 2π. We characterize the strength of the magnetic field via the plasma-β parameter, defined as the ratio of thermal to magnetic pressure, β = p g /p m , where p m is the magnetic pressure. We normalize the magnetic field strength by choosing the characteristic value of plasma-β to be high enough so that the initial accretion stage is gas pressuredominated, min β = 100.
(v) BH spin: We consider a rapidly spinning BH, with high dimensionless spin, a = 0.9375, to give the jets the best chance to form and survive in their fight against the onslaught of the infalling gas: if the jets form, their power will be close to the maximum possible value for a given magnetic flux (P jets ∝ a 2 , Blandford & Znajek 1977), enabling us to study the physics of their interaction with the infalling gas.
We add to the initial pressure random perturbations (independent for each numerical cell) at the level of 2%: without the perturbations, the system would maintain the exact axisymmetry and its evolution would be identical to a 2D model. We then let the system evolve out to t = 2.3×10 5 R g /c = 230R B /c = 0.1 Myr and report the simulation results. 3. RESULTS

Natural Development of X-shaped Jet Morphology
The simulation starts 2 with a uniform gas distribution outside the Bondi radius, and a vacuum hole at r < R B . The gravitational forces along with the pressure gradient at the interface of the hot gas and the empty cavity, at r = R B , push the gas inward. Because the gas possesses a non-zero angular momentum (Sec. 2), it undergoes what is nearly a freefall until it hits the centrifugal barrier at R circ = 30 R g . At this point, the radial infall slows down, and inside R circ an accretion disk forms. Because on the way to the BH the gas spends most of the time at the largest distances, we can estimate the time it takes for it to travel from R B to the BH as the free-fall time, t ff = 2 −1/2 (R B /R g ) 3/2 R g /c ≈ 2.2 × 10 4 R g /c. As seen in Fig. 1(a), around t ff the mass accretion rateṀ reaches its first peak and oscillates thereafter. Here,Ṁ = − ρu r dA, where dA = √ −gdθdϕ and g = |g µν | is the determinant of the metric. Analogously, we define the energy accretion rate, E = [(ρ + u g + p g + 2p m )u r u t − b r b t /4π]dA, where b µ is the fluid frame magnetic field 4-vector. We evaluate bothṀ anḋ E at r = 5R g , to avoid potential contamination by the density floors near the horizon. Because in a steady stateṀ andĖ are conserved and independent of radius, this does not affect time-average values and only slightly shifts the dependencies in time by ∆t 5R g /c, i.e., shorter than the sampling time of the simulation.
We use the sign convention such that positiveṀ andĖ imply mass and energy entering the BH. We also define the accretion efficiency as η = (Ė −Ṁ)/Ṁ. The gas drags with it the magnetic flux, and as Fig. 1(b) shows, this leads to an increase in the absolute magnetic flux on the BH, Φ BH = 0.5 |B r | dA, where the integration is over the area of the event horizon. We further normalize the magnetic flux by the time-smoothed mass accretion rate, φ BH ≡ Φ BH / Ṁ R 2 g c 1/2 . We smooth all quantities in Fig. 1 over a timescale of 1000R g /c using a 0th-order Savitzky-Golay filter, to make the plots more readable.
As the accretion disk forms, magnetized rotation of the BH works to form highly-magnetized jets. The increasing magnetic flux on the BH would ordinarily translate into magnetically-powered jets via the Blandford & Znajek (1977, BZ, hereafter) effect. However, the BH is engulfed in the infalling gas from all directions, which suppresses the event horizon electromagnetic (EM) luminosity of the jets, L, well below the analytic BZ expectation, L BZ : Fig. 1(c) shows that L/L BZ 0.25 at t 5 × 10 4 R g /c. In order for the jets to successfully launch and avoid falling apart due to gas-jet interaction and the development of magnetic kink instabilities, their total pressure needs to be higher than approximately the ram pressure of the infalling gas (Gottlieb et al. 2022). This happens only when the BH horizon accumulates enough magnetic flux, which powers the jets along with the spin of the BH, via the BZ process. The high-density gas circulation near the BH gives rise to magnetic flux polarity flips (Christie et al. 2019), and or this reason, at early times, t 5 × 10 4 R g /c, the jets continuously work intermittently, and get disrupted soon after launch. Figure 1(b) shows that, apart from a few bumps, the magnetic flux steadily increases until t 10 5 R g /c, and L/L BZ displays several prominent peaks reaching near unity, Figure 1(a) shows that the dips inṀ correlate with the spikes in the magnetic flux, luminosity and efficiency in Fig. 1(b)-(d). These are the moments when the jets are successfully launched with η ∼ 50% before the turbulent polarity flips disrupts them.
To better understand the behavior of the jets, we determine their tilt T and precession P angles, by computing the polar and azimuthal centroid positions, respectively, for both the northern and southern jets. For this, we introduce the total- Figure 1. To survive the onslaught of the infalling gas, jets deflect towards the equatorial plane (T ∼ 90 • , panel e), before aligning with the z-axis (T 10 • ). This jet reorientation naturally results in an Xshaped jet morphology (Fig. 2). The early period, t ≤ 5 × 10 4 Rg/c, of high mass accretion rateṀ/ṀB 0.06 (panel a) features low values of the dimensionless BH magnetic flux φBH 20 (panel b), total-EM to BZ jet-power ratio L/LBZ 0.25 (panel c), and jet energy efficiency η 10% (panel c). The infalling dense gas easily deflects the weak jets sideways, towards the equatorial plane, resulting in high jet tilt angle T ∼ 90 • (panel e) and large variations, ∼ 180 • , in precession P angle (panel f). At later times, t 10 5 Rg/c, the accumulation of largescale vertical magnetic flux on the BH leads to strong φBH ∼ 50 (panel b) and a MAD state (highlighted in green; steady state averages shown with horizontal dashed orange lines): the powerful jets, L/LBZ ∼ 1 (panel c) and η ∼ 150% (panel d), suppressṀ (panel a) and stably propagate along the z-axis (T 10 • , panel e), about which both the BH and the ambient gas rotate, out to large distances (Figs. 2 and 3). During the MAD state,Ṁ saturates at 2% ofṀB (panel a), implying that only 2% of the Bondi accretion rate reaches the BH and the rest 98% leaves as outflows. Figure 2. The first demonstration that X-shaped radio galaxy morphology naturally emerges from initially axisymmetric conditions: BH and ambient gas rotating around the vertical z-axis, with the gas threaded with a vertical magnetic field. Panels (a)-(i) show a time sequence of vertical slices through the simulated density (see color bar). The times and lengths shown in each panel have been scaled to M87*: 0.15 kpc = RB = 10 3 Rg, 475 yr = RB/c = 10 3 Rg/c. The early-time t 10 5 Rg/c ∼ 0.05 Myr evolution results in intermittent low-density jets (seen in green) that frequently disrupt due to the onslaught of infalling gas and transient reductions in jet power (see also Fig. 1). In spite of this, the jets and cavities they inflate manage to reach outside of the Bondi radius, which is shown with a black circle. Throughout the simulation, jets and remnant cavities point in different directions and resemble the X-shaped jet morphology in XRGs. The jets stabilize around t ∼ 10 5 Rg/c ∼ 0.05 Myr, once the BH saturates with the vertical magnetic flux and the disk enters the MAD state: the jets become strong enough to avoid getting deflected by the infalling gas and propagate along the vertical z-direction. energy to mass-energy flux ratio, µ = −u t (h + σ + 1) (Chatterjee et al. 2019), which gives the maximum Lorentz factor the flow could attain if all of its energy converted into kinetic energy. Here u t = g tν u ν is the covariant time-component of the 4-velocity, h = (u g + p g )/(ρc 2 ) is the specific enthalpy, and σ ≈ 2p m /(ρc 2 ) is the magnetization. We identify the jets on a sphere of radius 50R g using the condition µ ≥ 2, which selects the relativistic regions (jets) and eliminates non-and mildlyrelativistic regions (accretion flow and mildly relativistic outflows). Surprisingly, Fig. 1(e)-(f) shows that jet tilt and precession angles vary strongly as they erratically launch. In fact, at early times, t 5 × 10 4 R g /c, the jets launch nearly perpendicular to the polar axis, essentially in the equatorial plane. We see a clear anti-correlation between the inclination and precession angles (Fig. 1(e)-(f)) and the energy efficiency of the jets (Fig. 1(d)): when η peaks at t ∼ 7.5 × 10 4 R g /c, both T and P drop. Conversely, as η drops at t ∼ 10 5 R g /c, both T and P in-crease. This suggests that weaker jets are more easily deflected sideways by the pressure of the infalling gas. Figure 2(a)-(e) shows that the jets, which appear as underdense green regions, indeed propagate nearly horizontally away from the center. Their launch direction strongly deviates from that of the BH spin and ambient gas angular momentum vector directions, both of which are vertical and point along the z-axis. Shortly after launch, the jets subsequently quench. The timescale between the launching and quenching of the jets is approximately t ∼ 10 4 R g /c, or ∼ 0.01 Myr when rescaled to M87*. Furthermore, in a surprising fashion, the intermittent jets substantially wobble in time, both in tilt and precession angles, as the jets struggle to pierce through the dense surrounding gas.
Later on, at t 10 5 R g /c, Fig. 1(a)-(d) shows that the jet efficiency increases, and all quantities settle around their respective asymptotic values. The accretion rate reaches the low value of Ṁ /Ṁ B = 0.02 in units of the analytic spherical Bondi accretion rate while the magnetic flux asymptotes at the value φ BH = 50. The jet power levels off at L = 0.9 L BZ with an outflow efficiency reaching η = 150%. Figure 1 shows these steady state average values with horizontal dashed orange lines. Figure 1(e) shows that the jets manage to avoid getting deflected by the infalling gas and launch with a near-zero tilt, along the vertical z-axis. We can also see this in Fig. 2(c)-(i) that shows the jets propagating vertically outwards, displaying only slight bends, producing powerful backflows, and maintaining their overall stability.
The above high values of dimensionless BH magnetic flux, jet luminosity and efficiency, as well as the ability of the jets to avoid getting deflected by the infalling gas, make it clear that the jets are now strong enough to overcome the destructive effects of the infalling gas and operate at full strength. In fact, such high jet efficiency values, η 100%, are characteristic of the magnetically arrested disk (MAD) state in which the largescale vertical magnetic flux saturates the BH and becomes strong enough to periodically overcome the gravity pulling in the accreting gas and escape from the BH (Tchekhovskoy et al. 2011). These repeating magnetic flux eruptions manifest themselves through the fluctuations of the magnetic flux and jet efficiency in Fig. 1(b) and (d). For the rest of the simulation, we do not see any activity that indicates any deviation from this steady-state, although the rest of the duration is as long as the first highly-variable state.
The formation of wobbling jets is in contrast with the general expectation that when the BH spin and accreting gas angular momentum vectors are aligned, the accretion disk will produce jets that propagate along their direction. The reason for this discrepancy lies in the interaction of the jets with the surrounding gas, which deflects the jet head. This implies that even though this initial oscillating behavior of the jets might be a transient phenomenon, its effects can be long lasting and potentially probed as lower density cavities in the surrounding medium. The resulting morphology remarkably resembles the strongly asymmetric X-shaped jet morphology of XRGs (Sec. 1), even though the initial set-up is axisymmetric about the z-axis, apart from 2%-level pressure perturbations. Because realistic ISM has much stronger pressure fluctuations, any effect due to wobbling seen here will be much more pronounced for realistic ambient media. The cancellation of gas angular momentum near the BH can shut-off a jet outburst and power newly-formed pair of jets in a different direction. Eventually, enough magnetic flux accumulation on the BH leads to a magnetically arrested disk (MAD, Tchekhovskoy et al. 2011) and the power of the jet reaches its maximum value and manages to overcome the ambient gas density. In Fig. 2(e)-(i), the system is in its MAD state and not only do the newly formed jets remain stable near the highest density regions around the BH, but also retain their stability all the Figure 3. Northern jet shape (Rjet vs z) transitions cylindrical shape at early time to parabolic shape at late time (southern jet shows qualitatively similar behavior; not shown). Different colors and thickness indicate different times (see legend). At t 0.03 Myr, the jets exhibit a cylindrical shape at z 10Rg. At early times, the ram-pressure is high (largeṀ in Fig. 1a) and BH magnetic flux and jet power are low (small φBH and η in Fig. 1b,d), leading to strong collimation of the jets. At later times, after enough magnetic flux has accumulated on the BH, the jets become stronger, and take on a parabolic shape. We compute the cylindrical radius of the jet via the solid angle subtended by the jet (which we identify here as the region of µ > 1.5). way out to ∼ 7000R g 1 kpc, albeit the jet half-opening angle at these large scales is marginally resolved by about 7 cells.

Cylindrical and Parabolic Jets
Jets propagating through a dense medium can have rather different shape and stability compared to jets in vacuum. Observations indicate that AGN jets often show transition from the parabolic shape inside to conical shape outside the Bondi radius (Nakamura & Asada 2013;Kovalev et al. 2020;Boccardi et al. 2021). The origin of such a transition could either be due to the internal jet evolution (e.g., Nokhrina et al. 2020) or due to changes in the ambient density profile (e.g., Barniol Duran et al. 2017). Figure 3 shows the effective cylindrical radius of the northern jet, implied by the solid angle it subtends, versus the distance along the jet (the shape of the southern jet is similar and not shown). Since we use a logarithmic grid, at r 4000R g numerical diffusion in the jets increases the mixing, and thus we set the condition to isolate the jet to µ ≥ 1.5. In the very early intermittent jet phase, t 7 × 10 4 R g /c = 0.03 Myr, the jet is approximately cylindrical, R jet ∼ constant. At later times, while still in the early-phase and substantially tilted relative the z-axis, the jet transitions to a parabolic shape, R jet ∝ z 1/2 . The parabolic shape persists to late times, and the jets stably . Synthetic bolometric X-ray images (emissivity PX ∝ ρ 2 T 1/2 , neglecting absorption and scattering) demonstrate the formation of X-shaped jet morphology. Shock compressed and heated regions dominate the Bremsstrahlung emission, whereas empty cavities show jets or jet-inflated cavities. These low-emission regions present a morphology similar to X-shaped jets in XRGs. The color depicts the emissivity normalized by its angle-averaged value (Sec. 3.3). Each row and column corresponds to a different inclination angle and time, respectively. Times and lengths shown in each panel have been scaled to M87* (see also Fig. 2).
propagate from the BH along the z-axis out to z ∼ 7000R g = 7R B 1 kpc, without showing any tell-tale signs of disruptive kink-instability. Although the jets leave the Bondi sphere, they do not appear to exhibit the transition from parabolic to conical geometry at the Bondi radius. It is possible that longer and higher resolution simulations will show both the jet disruption, due to the kink instability, and the parabolic-to-conical transition, when the backflows are no longer able to reach from the jet head back to R B and affect the jet shape there.

Shocks and Cavities in Synthetic X-ray Images
Each panel of Fig. 2 indicates the position of the Bondi radius with a black circle and the spatial extent of the system with a black scale bar. Figure 2 shows that gas infall launches an expanding spherical accretion shock that forms as the gas finds itself rushing to the BH too fast. The shock reduces the gas radial infall velocity and helps the system to reach a steady state. Jets drive into the ambient gas an additional pair of bow shocks that by the end of the simulation elongate the accretion shock in the polar directions and start outrunning it.
The accretion and jet bow shocks compress and heat the ambient medium, resulting in thermal Bremsstrahlung X-ray emission. We are presenting a crude estimate of what the X-ray emission would look like, by calculating the Bremsstrahlung bolometric emissivity which scales as P X ∼ ρ 2 T 1/2 for fullyionized hydrogen. We assume an optically-thin gas and do not include any absorption or scattering effects. To construct synthetic X-ray images, we compute the projected fluid variables along the desired line of sight and use them to calculate the emissivity.
Columns in Fig. 4 depict a time series of synthetic Xray bolometric emissivity images. To highlight structure, we normalize the emissivity by its angle-averaged radial profile. Rows show views at different inclination angles relative to the z-axis, i = [15,30,90] • . At early times, in both the left and middle columns, both accretion shock and bow shocks are visible, leading to emissivity increase of ≈ 5−20% compared to the angle average (see the color bar). The bow shock encompasses the bulk of excess X-ray emissivity and surrounds the newlyformed jets. The dark regions show both the old jet-inflated cavities and newly-formed jets. The jet-inflated low-density cavities buoyantly rise, slower than the jets, and are soon outrun by them. The resulting X-shaped morphology is especially prominent in Fig. 4(b1),(b2),(c1),(c2): here, the jets are ap-proximately at the same distance from the SMBH as the older, jet-inflated cavities, forming a distinct X-shape. In the right column, even though the bow shock has outrun the accretion shock and the jets have reached ∼ 7R B 1 kpc, we can still see the old jet-inflated remnant cavities closer to the BH.

SUMMARY AND DISCUSSION
We have presented the first-ever 3D GRMHD simulation of a SMBH accreting rotating magnetized ambient gas from outside the Bondi radius, at a scale separation of R B /R g = 10 3 that is unusually large in 3D GRMHD. Initially located outside an empty cavity of radius R B , the uniform-density ambient gas contains a weak vertical magnetic field in the z-direction, which is aligned with both the BH spin and gas angular momentum vectors. As the gas falls in, it self-consistently forms a thick accretion torus whose size is set by the circularization radius of the ambient gas (eq. 1). More generally, the properties of the accretion system are controlled by several manifestly physical dimensionless parameters (see Sec. 2): (i) dimensionless Bondi radius, R B /R g , (ii) dimensionless circularization radius, R circ /R g , (iii) ambient gas plasma-β, (iv) dimensionless BH spin, a. This offers certain advantages over the standard equilibrium torus initial conditions typically used in GRMHD simulations (Fishbone & Moncrief 1976;De Villiers & Hawley 2003), where the parameters of the torus are harder to relate to observables.
As the infalling gas approaches the sonic surface too fast an accretion shock is developed. Around the same time, t 2 × 10 4 R g /c, a low-density intermittent jet launches close to the equatorial plane and disrupts soon thereafter. The disrupted jet remnants forms low-density bubbles, which buoyantly rise beyond the Bondi radius. The repeated disruption and recreation of the jets can be caused by: (i) the highly variable largescale vertical flux on the BH horizon (Fig. 1b), which leads to variability in the jet power (Fig. 1d) and even intermittent destruction of the jets; (ii) high ram pressure of the infalling gas on the BH horizon overcomes the ram pressure of the jets at times of low jet power and deflects them towards the equatorial regions. Indeed, minima in the φ BH and η in Fig. 1(b,d) correlate with the times of jet destruction, L/L BZ 1, seen in Fig. 1(c). At early times, when the mass accretion rate is high (largeṀ/Ṁ B ∼ 0.06) and jet ram pressure and power are low (η 50%), the tilt and precession angles of our jets reaches large values, T ∼ 90 • and P ∼ 180 • (Fig. 1e,f), implying that the infalling gas has deflected the jets towards the equatorial plane, i.e., perpendicular to the z-axis. Figure 2 shows that these early-time intermittent equatorial jets inflate buoyant bubbles whose direction is misaligned relative to the late-time stable jets. The synthetic X-ray images in Fig. 4 show multiple sets of jet-inflated cavities at different orientations. This morphology strongly resembles X-shaped jets in XRGs. Thus, the deflections of jets by infalling material can provide a nat-ural and simple explanation of XRG formation. Furthermore, the short-lived intermittent jet phase can explain the rarity of XRGs, which make up less than 1 in 10 radio galaxies.
At t 10 5 R g /c, the BH accumulates so much magnetic flux that the system goes MAD ( φ BH = 50) and the jets attain the maximum power for a given mass accretion rate (Tchekhovskoy et al. 2011). Namely, the jet luminosity is comparable to the BZ power, L ∼ L BZ (Fig. 1c), and the outflow energy efficiency reaches η = 150% (Fig. 1d). The accretion rate on the BH saturates at a level of Ṁ = 0.02Ṁ B (Fig. 1a), implying that 98% of the infalling gas at the Bondi radius is ejected from the system, with only 2% of the gas reaching the BH.
The Bondi accretion rate for M87* iṡ where λ s = 1/4 for γ = 5/3 (Shapiro & Teukolsky 1986;Di Matteo et al. 2003), and we assumed ionized hydrogen. Using this and the steady-state values for the system given above, the simulated mass accretion rate on the BH scaled to M87* becomes, in agreement with the inferred range of jet power, 10 42−45 erg s −1 (e.g., Stawarz et al. 2006;Broderick et al. 2015;Prieto et al. 2016;Nemmen 2019). Our expectation is that as we increase the scale separation from R B /R g = 10 3 closer to more physically-motivated values, 10 5−6 , more gas will be ejected in outflows, less gas will reach the BH, and the agreement with the observations will improve. We measure the shape of the jets to study whether our jets exhibit the parabolic to conical jet shape transition, which is observed near the Bondi radius in some AGN (Nakamura & Asada 2013;Kovalev et al. 2020;Boccardi et al. 2021). Figure 3 shows that whereas early-time jets have a cylindrical shape, at later times, t 6 × 10 4 R g /c, the jets transition to a parabolic shape. We do not see any signs of transition from parabolic to conical. We will investigate this issue in the future, by carrying out higher-resolution and longer-duration simulations, both of which might help to reveal the transition in jet shape (Sec. 3.2).
We have only studied a limited parameter space. In the future we will deploy a set of multiple high-resolution simulations and explore the parameter space, with what we think are the most crucial parameters. Specifically: (i) Choose a physically motivated density profile of ρ ∝ r −1 , which is associated with the ISM and ICM. (ii) Vary the scale separation R B /R g to lower and higher values, and quantify the differences in accretion rate efficiencies. (iii) Vary the circularization radius R circ /R g , and study the effects the size of the disk has on the accretion rate and wind power. (iv) Reach durations of t = 10 6−7 R g /c sufficient for the jets to transition into a conical shape, and/or fall victims to instabilities and get disrupted. (v) Include radiation transport and the effects of cooling on the dynamics of the accretion and power of the outflows. n introduce feedback in directions the jets normally can't access.