We carry out idealized threedimensional generalrelativistic magnetohydrodynamic simulations of prograde weakly magnetized and geometrically thick accretion flows where the gas distribution is misaligned from the black hole BH spin axis. We evolve the disk for three BH spins a  0.5 0.75 and 0.9375 and we contrast them with a standard aligned disk simulation with a  0.9375. The tilted disks achieve a warped and twisted steadystate structure with the outer disk misaligning further away from the BH and surpassing the initial 24 misalignment. However closer to the BH there is evidence of partial alignment as the inclination angle decreases with radius in this regime. Standing shocks also emerged in proximity to the BH roughly at 6 gravitational radii. We show that these shocks act to partially align the inner disk with the BH spin. The rate of alignment increases with increasing BH spin magnitude but in all cases is insufficient to fully align the gas before it accretes. Additionally we present a toy model of orbit crowding that can predict the location of the shocks in moderatetofast rotating BHs illustrating a potential physical origin for the behavior seen in simulationswith possible applications in determining the positions of shocks in real misaligned astrophysical systems.
['Black hole physics', 'Accretion', 'Magnetohydrodynamics', 'General relativity', 'Magnetohydrodynamical simulations', '159', '14', '1964', '641', '1966', 'Astrophysics - High Energy Astrophysical Phenomena']
Shockinduced Partial Alignment in Geometrically Thick Tilted Accretion Disks Around Black Holes
{'Shock-induced partial alignment in geometrically-thick tilted accretion disks around black holes': '<!-- image -->  \n1 JILA and Department of Astrophysical and Planetary Sciences, University of Colorado, Boulder, CO 80309, USA', 'ABSTRACT': 'We carry out idealized three-dimensional general-relativistic magnetohydrodynamic (GRMHD) simulations of prograde, weakly magnetized, and geometrically thick accretion flows where the gas distribution is misaligned from the black hole spin axis. We evolve the disk for three black hole spins: a = 0 . 5 , 0 . 75, and 0 . 9375, and we contrast them with a standard aligned disk simulation with a = 0 . 9375. The tilted disks achieve a warped and twisted steady-state structure, with the outer disk misaligning further away from the black hole and surpassing the initial 24 · misalignment. However, closer to the black hole, there is evidence of partial alignment, as the inclination angle decreases with radius in this regime. Standing shocks also emerged in proximity to the black hole, roughly at ∼ 6 gravitational radii. We show that these shocks act to partially align the inner disk with the black hole spin. The rate of alignment increases with increasing black hole spin magnitude, but in all cases is insufficient to fully align the gas before it accretes. Additionally, we present a toy model of orbit crowding that can predict the location of the shocks in moderate-to-fast rotating black holes, illustrating a potential physical origin for the behavior seen in simulations-with possible applications in determining the positions of shocks in real misaligned astrophysical systems.  \nKeywords: black hole physics - accretion discs - magnetic field - Magnetohydrodynamics(1964) General relativity(641)', '1. INTRODUCTION': "Hot accretion flows have long been of astrophysical interest. These systems, in contrast to their thin-disk counterparts (Shakura & Sunyaev 1973), exhibit significantly reduced luminosity and operate with accretion rates far below the Eddington limit ( ˙ M ≪ ˙ M Edd ) (Xie & Yuan 2012; Abramowicz & Fragile 2013). This makes them suitable to explain the low/hard state of black hole X-ray binaries (Esin et al. 1997; Dexter et al. 2021) and in low-luminosity Active Galactic Nuclei (AGN) including the supermassive Black Hole in our Galactic Centre, Sagittarius A* (Sgr A*) (Yuan et al. 2002; Narayan & McClintock 2008; Yuan & Narayan 2014).  \nTraditionally such systems are assumed to have their angular momentum aligned with the spin of the Black  \nCorresponding author: Sajal Gupta  \nSajal.Gupta@colorado.edu,jason.dexter@colorado.edu  \nHole (BH). This presumption holds valid for thin, viscous accretion disks with small tilt angles, allowing viscous torques to gradually bring the inner flow into alignment (Bardeen & Petterson 1975; Liska et al. 2019). However, in the case of large tilt, nonlinear behavior such as the tearing of a disk are observed (Nixon & King 2012; Nealon et al. 2015; Liska et al. 2021; Kaaz et al. 2023), and it remains uncertain whether BardeenPetterson alignment can occur under these conditions. Alternatively, if the system has a high accretion rate, the spin of the BH may grow and eventually coincide with the angular momentum of the inflowing material (Volonteri et al. 2005). Furthermore, magneto-spin alignment driven by powerful relativistic jets can also play a crucial role in aligning the strongly magnetized disk with the BH spin (McKinney et al. 2013; Chatterjee et al. 2023). However, the low accretion rates and geometrically thick structure of hot accretion flows such as those observed in Sgr A* and M87 ( ˙ M ∼ 10 -7 -10 -9 M ⊙ /yr for Sgr A ∗ (Dexter & Fragile 2013; White et al. 2020))  \nmake it unlikely that the Bardeen-Petterson effect, cumulative mass accretion, or magneto-spin alignment can induce alignment. Consequently, these systems are more likely to remain misaligned, impacting various physical processes. These include growth rates of supermassive BHs (King et al. 2005; Stella & Vietri 1997), their spin measurements, and observable accretion rates. Additionally, the disk misalignment also affects the structure and size of the BH shadow (Dexter & Fragile 2011; White et al. 2020; Chatterjee et al. 2020) along with the formation and direction of the relativistic jet (King & Nixon 2018; Liska et al. 2019). For instance, recently Cui et al. (2023) observed that the radio jet in M87 precesses with a ∼ 11-year period, attributing this behavior to the misalignment of hot accretion flows.  \nPapaloizou & Lin (1995); Nelson & Papaloizou (2000) were among the first to study the evolution of misaligned disks. They developed a linear theory focused on the bending wave regime, and conducted three-dimensional smoothed particle hydrodynamics (SPH) simulations, demonstrating that warps in the low-viscous flows propagate in a wave-like manner rather than diffusing. This causes the disk to remain misaligned even at inner radii, with the disk's plane oscillating as a function of radius (Ivanov & Illarionov 1997; Lubow et al. 2002). Later, Zhuravlev & Ivanov (2011), while developing an analytical model to describe the structure of twisted disk around slowly rotating Kerr BHs, proposed that these oscillations strongly depend on viscosity in the prograde flows and are notably absent for disks with viscosity parameter, α (Shakura & Sunyaev 1973), between 10 -3 and 10 -2 , and effective scale-height, H/r/ √ a ∼ 10 -2 , typically lower than those found in numerical studies of geometrically-thick disks (Fragile et al. 2007; White et al. 2019). Here, a ≡ cJ BH /GM 2 BH , ranging from 0 to 1, is the dimensionless BH spin parameter. Recently, White et al. (2019) applied the linear theory to their GRMHD simulations of tilted accretion flows and observed no such oscillations. Instead, the linear theory predicted an exponential increase in disk tilt with decreasing radius, whereas their study showed that disk modestly warps further away from alignment at large radii before a rapid decline in tilt is noticed, highlighting both the limitations of linear bending wave theory and the tendency of the inner flow to partially align. Throughout our study, we refer to this phenomenon as 'partial alignment', denoting the reduction in inclination angle of the inner disk in its steady-state configuration. Previously, Sorathia et al. (2013a,b); Hawley & Krolik (2018) question the relevance of linear theory as they observed that MHD turbulence in misaligned flows disrupts the radial communication of warps, stressing  \nthe importance of strong warps driving non-linear effects such as shocks in the system.  \nFragile et al. (2007) carried out simulations of a weakly magnetized hot accretion flow onto a BH of spin a = 0 . 9, and discovered a pair of standing spiral shocks forming close to the BH (Fragile & Blaes 2008). The location of shocks found in their study, along with other GRMHD simulations (White et al. 2019; Generozov et al. 2013), coincides with regions of rapid fall in inclination of the inner flow, raising speculation that the shocks might be responsible for the partial alignment of the inner disk (Mewes et al. 2016; Kawaguchi et al. 2015). Despite this spatial correlation, no conclusive evidence has yet been presented that shocks act to increase the alignment of the inner disk.  \nMoreover, the potential mechanism for the formation of shocks in thick, tilted accretion disks is not yet fully understood. Previous studies led by Fragile & Blaes (2008) and Generozov et al. (2013) integrated the fluid trajectories in post-processing and suggested that the convergence of eccentric fluid orbits near the apocenters may be responsible for the shocks formation. However, their analysis were not able to explain the unique geometry of the shocks. It is also unclear how the spatial structure of shocks is affected by the bending and twisting of the disk as well as the spin of the BH. In light of these uncertainties, in this study we attempt to address this speculation of shock-induced partial alignment, and we describe the role of standing shocks in shaping the structure of the disk as well as illustrating a potential mechanism of their origin.  \nTo this end we perform GRMHD simulations. Despite the fact that prior 3D SPH simulations have shown a warp in disk shape (Lubow et al. 2002; Nealon et al. 2015; Drewes & Nixon 2021), none have included a magnetic field in the system. In addition, none of them have observed shocks, which appear to be an important non-axisymmetric feature in thick tilted disks. Previous GRMHDsimulations from Fragile and collaborators utilize an artificial viscosity scheme to capture the effects of shocks, whereas the parameter survey led by White et al. (2019) was ran for insufficient time to attain inflow equilibrium across the inner disk.  \nWe carry out a series of idealized 3D GRMHD simulations of weakly magnetized, geometrically-thick, tilted accretion disks called standard accretion and normal evolution (SANE) disks, initially making an angle of 24 · with the BH equatorial plane. We evolve the disks for three different BH spins: a = 0.5, 0.75, and 0.9375, allowing us to study the effect of BH spin on the disk structure. To compare the characteristics of misaligned disks with the conventional picture, we simulated an ac-  \netion disk with no initial tilt onto a BH of spin a = 0.9375. We self-consistently evolve an entropy equation for the electrons using the Ressler et al. (2015) scheme. The latter enables us to calculate the irreversible heating rate on the fly, enabling us to resolve the spatial structure of the shocks. All the cases are evolved with the same initial conditions, thus making the comparison between different models free of this variable. A brief discussion on the simulation setup and initial conditions are given in section 2. We provide our findings in section 3, and we conclude and discuss our results in section 4.", '2. SIMULATION SETUP': "We performed idealized 3D GRMHD simulations using the publicly available HARMPI code 1 (a parallel, 3D version of HARM (Gammie et al. 2003; Noble et al. 2006)) that solves the GRMHD equations for a SANE disk around a rotating BH defined by the Kerr metric with metric signature ( -,+,+,+). We initialize our simulations with a steady-state hydrodynamic equilibrium gas torus (Fishbone & Moncrief 1976), with inner radius and pressure maximum at r in = 12 r g and r max = 25 r g . Here, r g ≡ GM BH /c 2 denotes one gravitational radii. The torus was seeded with a single loop of poloidal magnetic field, thus facilitating the magneto-rotational instability (MRI) to kickoff. The internal energy and the gas pressure are related by adiabatic index: P gas = (Γ -1) u g , where we chose Γ = 5 / 3. The magnetic field strength is normalized by fixing plasma β -parameter, β =max P gas / max P b = 100, where P b = b µ b µ / 2 is the magnetic pressure. To study the misaligned flows, we rotate the torus with 24 · about the BH equatorial plane, such that the initial angular momentum of the disk lies in the X-Z plane. We evolve the system in modified spherical-polar KerrSchild coordinates ( r, θ, ϕ ), with a grid resolution of 320 × 256 × 160 cells to adequately resolve the MRI turbulence. A higher concentration of cells is implemented in the BH equatorial plane to resolve the gas flow and at the pole to resolve the jet at larger radius. Due to the disk's tilt, orbital parcels traverse regions with varying polar resolution. Specifically, the change in polar resolution reaches up to ∆ θ ∼ 0 . 05 radians over the disk extent. In all the cases we studied, the numerical setup covers the region of (0.88 r H , 10 5 r g ) × (0, π ) × (0,2 π ) in radial, polar and azimuthal directions. Here, r H = r g (1 + √ 1 -a 2 ) is the BH event horizon radius. The outer radial boundary was extended to 10 5 r g using a superexponential radial coordinate. In our simu-  \ne boundary conditions (BCs) in the θ -direction. Unlike transmissive polar BCs, which allow fluid to pass through the boundaries with minimal interference (Liska et al. 2018), reflective BCs can induce additional dissipation near the polar axis and potentially enlarge the jet artificially, influencing disk alignment (McKinney et al. 2013). To address these potential issues, we specifically exclude the polar axis region from our averages and remove regions with a magnetization ( σ ≡ b 2 ρ ) greater than one.  \nWe utilized Ressler et al. (2015)'s approach and computed the total heating rate at each time step by measuring the difference between the entropy calculated through energy conservation equations and that from entropy conserving equation. Magnetic and/or kinetic energy dissipated at the grid scale due to the turbulence is recaptured as internal energy, which we use to identify shocks. Our procedure is similar to White et al. (2019), with the difference being that we calculated the heating rate at each simulation time-step, which allowed us to incorporate the effects of time-rate change of entropy.  \nTo determine the radial dependence of physical quantities, we follow Penna et al. (2010) averaging approach, and define a density-weighted, time-averaged of quantity X as:  \n⟨ X ⟩ ρ ( r ) = ∫ ∫ ∫ X ρ ( r, θ, ϕ ) dA θϕ dt ∫ ∫ ∫ ρ ( r, θ, ϕ ) dA θϕ dt (1)  \nwhere dA θϕ ≡ √ -gdθdϕ is an area-element in spherical polar Kerr-Schild coordinates, and the integral over dt is a time average in the steady-state. The verticallyintegrated quantities are computed as  \n⟨ X ⟩ ρ,θ ( r, ϕ ) = ∫ ∫ X ρ ( r, θ, ϕ ) √ -gdθdt ∫ ∫ ρ ( r, θ, ϕ ) √ -gdθdt (2)  \nWe measured the disk thickness using:  \nH/r = √ 〈 ( θ ' -θ 0 ) 2 〉 ρ,θ ' ,ϕ ' (3)  \nwhere θ 0 ≡ π/ 2 + ⟨ ( θ ' -π/ 2) ⟩ ρ,θ ' ,ϕ ' is the disk midplane location. Here, ( θ ' , ϕ ' ) are the angular coordinates aligned with the disk, which are related to the coordinates ( θ, ϕ ) via Equation B1. Qualitatively, θ ' = π/ 2 represents the local midplane of the disk, and ϕ ' = 0 is set to the local precession angle. The notation ⟨ X ⟩ ρ,θ ' ,ϕ ' denoted a density-weighted, timeaveraged value of quantity X in disk aligned (or tilted) coordinates. At later times, the above equation yields H/r ∼ 0 . 3 -0 . 4 inside r = 20 r g .  \nFigure 1. Snapshots from the equilibrium states contrasting the properties of a hot, aligned disk (top row) with tilted accretion flows at BH spin, a ≊ 0 . 94 (bottom row). The color scales are in arbitrary/code units. Bending of the disk is visible in the poloidal slice of the density ( ρ ) and the magnetization ( σ ) (leftmost column), showing the disk midplane (white line), on average, is tilted with respect to equatorial plane of the BH (Z BH = 0). Non-axisymmetric structures in density-weighted, verticallyaveraged snapshots of heating rate ( ˙ q , middle column) signify the standing shocks, which extend well beyond the r ISCO (white dashed circular line) and correlates with low specific angular momentum (( l = -u ϕ /u t , rightmost column) regions. Here r ISCO denotes the innermost stable circular orbit (ISCO) -the smallest possible radius for a stable circular orbit.  \n<!-- image -->", '3.1. General Characteristics: Consequences of an initial tilt': 'The focus of this study is to investigate the steadystate properties of tilted flows. We evolve the simulations until a quasi-steady state is achieved, which we determined by monitoring the density-weighted radial profiles of tilt and twist, as well as the time it takes for standing shocks to fully develop. By these criteria, the tilted disks attain the inflow equilibrium up to r ∼ 20 r g within t ∼ 12000 r g /c . Tilt reaches equilibrium around t ∼ 10000 r g /c , while twist profiles develop from the beginning due to frame-dragging precession. We run the  \nhigh spin case (a = 0.9375) for t ∼ 24000 r g /c , to be confident about the steady-state whereas the untilted disk is evolved only for t ∼ 10000 r g /c , as we notice negligible changes in relevant physical parameters.  \nIn Figure 1, we showcase and compare misaligned accretion disk profiles with the standard aligned case for spin a ≊ 0 . 94. The general trend of the disk bending remains similar across other spins, as demonstrated in Appendix A. The 2D snapshots of different physical quantities for both scenarios are averaged during their respective steady states: from 8000 ≤ t ≤ 10000 r g /c for the aligned case (top row) and from 20000 ≤ t ≤ 22000 r g /c for the misaligned case (bottom row). Since both sim-  \n<!-- image -->  \nFigure 2. Density-weighted, shell averaged radial profiles of tilt (left) and twist (right) from each simulation (colors). The twist increases as the flow progresses closer to the BH due to frame-dragging precession, which increases with the spin of the BH. The increasing inclination with decreasing radius indicates that the disk warps further away from alignment as it accretes before beginning to align close to the BH (inside ∼ 6 r g ). This rate of partial-alignment also increases as the spin of the BH increases.  \n<!-- image -->  \nchieved steady states, averaging at different times does not represent distinct evolutionary stages.  \nThe left column shows the poloidal slice of the restmass density ( ρ ) and the magnetization ( σ ≡ b 2 ρ ) for aligned and misaligned flows. Comparing these profiles reveals a pronounced disk bending in the tilted flow, indicating persistent misalignment at the steady-state. Perpendicular to the disk midplane (white arrow) is its angular momentum vector ( L disk ), pointing toward the low-density, high-magnetization region called the funnel region. This suggests that jets may emit in the direction of disk angular momentum, and may precess with the disk as remarked by Liska et al. (2018). Near the BH, the tilted flow shows two separate density arms, which are the regions of high compression formed due to standing shocks (Fragile & Blaes 2008; Generozov et al. 2013).  \nTo identify standing shocks and their effects, we plot density-weighted, vertically averaged snapshots of the total heating rate ( ≡ ˙ q ), and the specific angular momentum ( l ≡ -u ϕ /u t ) in the central and right columns of Figure 1. In the aligned case (top row), both ⟨ ˙ q ⟩ ρ,θ and ⟨ l ⟩ ρ,θ are axisymmetric as expected. However, in the tilted flow (bottom row), the structure completely changes, and we find significant non-axisymmetric features in both quantities. Furthermore, the deviations from axisymmetry increase with BH spin, as illustrated in Figure 8.  \nThe m = 2 azimuthal structure in heating rate ( ˙ q ) plots (cf. figure 1) signifies the pair of standing, spiral shocks. Here the m = 2 mode signifies the Fourier decomposition of the azimuthal structure in the heating rate. These shocks are localized, stronger, and have a more well-defined structure at higher spins. Comparing the heating rate and specific angular momentum in the tilted flow reveals that regions of depleted angular momentum correspond to the standing shocks, suggesting that they may enhance the transport of angular momentum. The geometry of the shocks in the disk-aligned frame is illustrated in Appendix E.', '3.2. Inclination and twist': "We measure the tilt as a function of disk radius by defining inclination as:  \nI ( r ) ≡ cos -1 ( J BH · L disk ( r ) | J BH || L disk ( r ) | ) (4)  \nwhere J BH = a ˆ z is the angular momentum of BH in natural units. The angular momentum of the disk: L disk is defined as:  \n( L disk ) ˆ k = ϵ ijk r ˆ j ( T t ˆ k ) MA (5)  \nwhere ϵ ijk is the three-dimensional anti-symmetric Levi-Civeta symbol, and the hats indicate the cartesian coordinates related to spherical coordinates in the conventional way. We call ( L x , L y , L z ) ≡  \nFigure 3. Density-weighted, shell averaged radial profiles of tilt (top row), dimensionless warp ( ˆ ψ ) (middle row) and warp amplitude (bottom row) for tilted simulations (colors). At small radii (shaded areas), the simulations exhibit increased warping ( ˆ ψ ( r ) ≥ 1), indicating the presence of non-linear effects, which coincide with regions of decreasing tilt (top row).  \n<!-- image -->  \n( L disk (ˆ x ) , L disk (ˆ y ) , L disk (ˆ z ) ). Using the above information, we can write density-weighted, average inclination as:  \nI ( r ) = cos -1   ⟨ L z ⟩ ρ √ ⟨ L x ⟩ 2 ρ + ⟨ L y ⟩ 2 ρ + ⟨ L z ⟩ 2 ρ   (6)  \nand the twist as:  \nT ( r ) = tan -1 ( ⟨ L y ⟩ ρ ⟨ L x ⟩ ρ ) (7)  \nThe inclination I ( r ) is the measure of how elevated the disk midplane is from the BH equatorial plane, whereas the twist T ( r ) is the measure of the precession of the disk's angular momentum about the spin axis of BH. If the inclination and/or twist is a function of radius, the disk acquires a warped structure and is referred to as a warped disk as mentioned in section 1.  \nFigure 2 shows the radial, steady-state profiles of averaged inclination ( I ( r )) [left panel] and twist ( T ( r )) [right panel]. The radial profiles of angles are averaged in their steady-state, particularly, 8000 ≤ t ≤ 10000 for the untilted case, 12000 ≤ t ≤ 14000 for a = 0 . 5, and 0 . 75, and 20000 ≤ t ≤ 22000 for the high spin case ( a ≊ 0 . 94). To get a sense of the accretion flow, it is best to read the plot from right to left. The twist sets up very early in the simulation due to frame-dragging precession and later relaxes to a smaller value due to the transport  \nof angular momentum (White et al. 2019). The stronger twist in the high spin case is because frame-dragging precession is proportional to the spin of the BH (Lense & Thirring 1918). Overall, the shape of our twist profile indicates that the disk becomes extremely twisted as one approaches the BH.  \nIn the left panel of Figure 2, a tiny tilt of approximately 3 · (grey line) is noticeable in the aligned case. Ideally, we expect a flat line of 0 · . This discrepancy arises due to numerical errors associated with taking the projections of the angular momentum of the disk on the BH equatorial plane, as the out-of-plane angular momentum projections, namely L x and L y , are not conserved. Nevertheless, the associated error is small, suggesting that our routine for determining the inclination is accurate to within about O (3 · ). Our definition also yields the initial tilt and twist far away from the BH, as we expect the tilt and twist to asymptote to an initial value of 24 · and 0 · , which we can see is true for all cases.  \nAnalyzing the tilted profiles reveals that inclination increases as one approaches the BH, implying that the disk progressively bends further away from the BH's equatorial plane. Given the simulations' duration, which exceeds the inflow timescale, we conclude that the disk's warping is a stable characteristic, associated with bending waves. However, the predicted high-frequency oscillations from linear theory are virtually absent in our  \ntilt profile. This could be due to the small viscosity in our flows, which may suppress the large oscillations and mitigate the radial communication of warps (Zhuravlev & Ivanov 2011).  \nThe rapid decline in inclination, especially at very small radii ( r ≤ 6 r g ), causes the disk to partially align with the BH. This descent occurs at a significantly faster rate when the BH has a high spin, around a ≊ 0 . 94, indicating a direct influence of the BH's spin on the inner disk's physical structure. In our observations, we also confirm the presence of shocks in the disk. However, before quantifying their impact on the tilt profile, if any, it is essential to determine their radial extent.  \nWe achieve this by identifying the region where the dimensionless warp parameter, represented as ˆ ψ = ψ/ ( H/r ), exceeds 1. Here ψ ( r ) ≡ ∣ ∣ ∣ rd l ( r ) dr ∣ ∣ ∣ signifies the measure of local warp and l ( r ) = L disk ( r ) / | L disk ( r ) | is the unit vector in the direction of disk's angular momentum at radius r. In the presence of small warp ( ˆ ψ > 0, but ≪ 1), strong oscillatory motions are generated in a tilted disk due to warp-induced radial pressure gradients (Ogilvie & Latter 2013). However, when ˆ ψ ( r ) > 1, nonlinear features can emerge. This occurs because, in this regime, local warp surpasses the disk's height, causing these oscillatory motions to couple with the warp. As a result, the radial velocity approaches the sound speed (Sorathia et al. 2013a,b), leading to compression.  \nFigure 3 shows that ˆ ψ ( r ) > 1 occurs in the inner flows for spin a = 0 . 75 and 0 . 94 (shaded regions). This aligns neatly with the region of partial-alignment and coincides with the area where shocks form, as observed in the Figure 1. Notably, we do not identify any region with ˆ ψ ( r ) > 1 for the a = 0 . 5 case, as the warping in the lowest spin scenario is less pronounced due to minimal twisting and bending of the disk. In this case where ˆ ψ ( r ) < 1 we do not presume the absence of non-linear effects, but rather that they are much weaker than in the higher spin cases. It is essential to note that the condition ˆ ψ ( r ) > 1 is based on simple physical arguments and should not be interpreted rigidly, i.e. Sorathia et al. (2013b) chose an upper limit of 0.8 instead of 1.  \nFigure 4 visualizes the averaged disk plane, represented as a series of rings observed from the perspective of the BH frame. Each ring is characterized by its individual tilt and twist, in accordance with the radial profiles displayed in Figure 2. Notably, the figure distinctly showcases the pronounced warping of the disk, particularly evident in the vicinity of the black hole where the twisting is more intense.  \n3.3. Origin of the shocks  \nFigure 4. Representation of the averaged disk plane, depicted as tilted and twisted rings for the a ≊ 0 . 94 case. The visualization emphasizes pronounced warping, particularly near the black hole (black sphere). The BH spin axis is vertical in the figure, where the viewpoint camera is located at an elevation of -130 · and an azimuth of -90 · with respect to the BH equatorial plane.  \n<!-- image -->  \nTo understand the formation of shocks and how they are affected by radial variations of tilt, and twist, we employed two different approaches. First, we draw insights from the work of Ivanov & Illarionov (1997) and developed a toy model following the prescription below:  \n- 1. We define a series of centered elliptical orbits with an eccentricity of e = 0 . 1 2 in a disk-aligned frame, with the BH at the center. The diskaligned frame 3 is oriented such that its z-axis ( ≡ Z ' ) is parallel to the disk's angular momentum. The x-axis ( ≡ X ' ) is aligned with the line of nodes, and the y-axis ( ≡ Y ' ) completes the righthanded coordinate system. Initially, these orbits are untwisted, centered and are constrained to the X ' -Y ' plane.\n- 2. We then represent these orbits in the BH reference frame by calculating the angular coordinates ( θ, ϕ ) using Equation B1. This process manifest the geometry of the averaged tilted disk as seen from the BH frame. Figure 4 shows the orientation of such  \nFigure 5. Density-weighted, vertically integrated snapshots of the heating rate ( ˙ q ) with centroid positions (black) for all tilted flows, overlaid with locations of maximum crowding of orbits from the toy model (blue) and the predicted shock locations from the basis vector model (green). The close agreement between the centroid locations and predicted shocks location from both methods at BH spins, a = 0 . 75 and ≊ 0 . 94, suggests that the shocks form due to the crowing of elliptical orbits in the flow. The significant discrepancy between the black and blue/green points for radii r > 6 r g is attributed to the large radial variation in the tilt profile compared to the twist.  \n<!-- image -->  \norbits viewing at a specific angle for the a ≊ 0 . 94 case.  \n- 3. To locate the regions of maximum crowding, we identify the points where the distance between neighboring rings is minimized. These points are labeled as the sites of maximum crowding, which we predict to be the locations where shocks form.\n- 4. Finally, we superimpose the predicted sites of maximum crowding on the plot of verticallyintegrated heating rate ( ⟨ ˙ q ⟩ ρ,θ ) and compare them with observed shock-forming regions.  \nSecond, based on the work of Ogilvie (1999); Ogilvie & Latter (2013), we defined local basis vectors for our warped disk with a few adjustments:  \n- 1. m = rd l ( r ) dr , a unit vector orthogonal to l and pointing where adjacent warped annuli are most separated.\n- 2. ˆ n = l × m | l × m | , a unit vector pointing where adjacent warped annuli are least separated, and thus should point to the location of the shocks.\n- 3. Given the m = 2 nature of shocks, we also considered the vector -ˆ n , as both ˆ n and -ˆ n should point towards the least separated annuli.\n- 4. Finally, we superimposed the x and y components of both vectors, r ˆ n and -r ˆ n on the ( ⟨ ˙ q ⟩ ρ,θ ) plot,  \nwhere we multiply by ' r ' to account for the shocks differing locations at different radii.  \nThe findings from both our toy model and basis vector model are illustrated in Figure 5, where the black dots represent the centroid of vertically-averaged heating rate, and blue and green dots represent the predicted sites of shock location using the toy model and basis vector model, respectively. Both methods predict the shock locations quite well, especially for higher BH spins ( a ≊ 0 . 94 and 0 . 75). However, mapping the shock structure in the heating rate plot for the low spin case requires a comprehensive fluid trajectory analysis, as done by Generozov et al. (2013), for a meaningful comparison. Because shocks are weaker at lower spins, their m = 2 structure is less distinct. Consequently, using the centroid positions of the vertically-integrated heating rate as a proxy for shock locations introduces errors, particularly at larger radii and for the lowest spin.  \nSince both methods are majorly governed by tilt and twist values, and their radial variation, we found that in the very vicinity of the BH, ( r ≤ 6 r g ), twist is the dominant factor in determining the spatial structure of the shocks (see Appendix E) Whereas at larger radii ( r > 6 r g ), the observed disparities between our sites of maximum crowding and the centroid positions are primarily due to the radial dependence of the tilt profile, as at these radii the twist in the disk changes more gradually with radius compared to the tilt. Thus, for radii r > 6 r g , the predicted shock locations are significantly  \nFigure 6. Density-weighted, vertically integrated snapshots of change in the rate of inclination for BH spin, a ≊ 0 . 94. The sound agreement between the shock-forming locations depicted by centroid of ( ⟨ ˙ q ⟩ ρ,θ ) (black) and the regions where the tilt is decreasing demonstrates the possible role of shocks in partially aligning the inner disk with the BH.  \n<!-- image -->  \ninfluenced by the radial variation of the tilt profile. Our observations also indicate that while eccentric orbits are essential in our toy model for identifying regions of local density enhancement, the geometry of these shocks remains largely unaffected as long as the eccentricity remains small ( e < 0 . 2).", '3.4. Shock-induced partial alignment': "To understand the decrease in inclination within the inner flow, we investigate the rate of change of inclination along the path of a fluid parcel. This involves calculating the material derivative of inclination ( I ) using:  \nd I dt = ∂ I ∂t + V · ∇I (8)  \nwhere V is the coordinate three-velocity, the components of which are given by: V i = u i /u t .  \nGiven our focus on the steady-state behavior of tilted flows, the first term on the right-hand side of the equation becomes zero. Expanding the second term yields three additional terms that describe the change in inclination of a fluid parcel subjected to a space-dependent velocity field. To assess the influence of shocks, we compute the vertically-integrated rate of change of inclination, represented as:  \n⟨ V · ∇I⟩ ρ,θ = ⟨ V r ⟩ ρ,θ ∂ I ∂r + ⟨ V ϕ ⟩ ρ,θ ∂ I ∂ϕ (9)  \nAll the quantities in this equation are density-weighted and vertically averaged within the σ ≤ 1 region to only consider the dense area of the flow and reduce errors  \nassociated with high radial velocities in polar regions. In the above equation I is calculated using:  \nI ( r, ϕ ) = cos -1   ⟨ L z ⟩ ρ,θ √ ⟨ L x ⟩ 2 ρ,θ + ⟨ L y ⟩ 2 ρ,θ + ⟨ L z ⟩ 2 ρ,θ   (10)  \nFigure 6 illustrates the evolution of inclination for a fluid parcel moving in the Eulerian frame for the highspin case ( a ≊ 0 . 94). We chose the high spin case because the shock structure is well-defined, and the decline in inclination is much more pronounced compared to lower spins. The plot shows that the region with a negative rate overlaps with the shock-forming region, identified by black dots from the centroid positions of vertically-averaged heating rate. Additionally, the rate exhibits a non-axisymmetric structure similar to the heating rate [cf. Figure 1]. This indicates that the tilt of a fluid parcel decreases as it passes through a shock front, demonstrating the critical role of shocks in reducing the tilt of inner flows. However, fluid parcel also experiences an increase in tilt as it moves out of shock forming regions. This feature, where d I dt is positive, can be attributed to the angular momentum fluxes along with the gravitational torque.  \nHowever, it is important to note that the centroid positions do not mark the exact location of the shocks front. Moreover, calculating the rate of change of inclination in the Eulerian frame involves numerical derivatives, which can introduce small errors. Despite these limitations, the regions with decreasing inclination are consistent with the shock-forming regions. The rate of  \nFigure 7. Plot comparing the alignment timescale (blue) to the accretion timescale (orange) within the disk. Both timescales are normalized by the local orbital period ( t dyn ). The fact that the gas is accreting much faster than it is aligning indicates the inability of the inner disk to fully align with the spin of the BH.  \n<!-- image -->  \nchange of inclination for a = 0 . 5 and 0 . 75 are shown in Figure 9.  \nFurthermore, we observe that the second term on the right-hand side of equation 9-which represents the flux transport of the change in inclination of a gas parcel through azimuthal motion, is dominant and largely determines the V ·∇I . This causes the sign of d I /dt to flip at the shock position. To ensure the robustness of this result, we performed additional analyses where we defined the rate of change of inclination in three different ways, as detailed in Appendix D. Although these definitions showed minor differences, the overall features of the rate of change of inclination variable remained consistent.  \nA natural question arises: if shocks are promoting alignment within the inner disk, why do we not observe complete, mutual alignment between the disk and the black hole, akin to the Bardeen-Petterson effect? Why is partial alignment the prevailing outcome? To address this question, we introduce an alignment timescale: t align ≡ I ( r, t ) / ( -d I ( r, t ) /dt ). The inclusion of a minus sign ensures a positive timescale, considering our expectation that d I ( r, t ) /dt < 0 during the alignment process.  \nWe compared this timescale with the inflow timescale: t inflow = r/ ( -⟨ V r ⟩ ρ ). Our hypothesis is that the disk does not completely align because the gas is forced to flow radially inward at a much faster rate than the time it requires to dissipate its out-of-plane angular momentum entirely. Thus, we examine the relationship between the alignment timescale and the inflow/accretion timescale. In the above definitions, ⟨ V r ⟩ ρ represents density-weighted, time averaged radial velocity. The alignment timescale is calculated by taking the ratio of  \nradial inclination profile (cf. Figure 2) to the azimuthal average of d I /dt that we calculated using Equation (8).  \nFigure 7 illustrates these timescales as a function of radius for the misaligned flow around the BH of spin, a ≊ 0 . 94. The y-axis is normalized by the dynamical timescale ( t dyn = 1 / Ω, where Ω = dϕ dt = ⟨ V ϕ ⟩ ρ ), representing time in periods of flow rotation. The figure depicts that in the inner regions where the tilt is decreasing, t align exceeds t inflow by a factor of ∼ 2 -100. This indicates that gas lacks the required time to attain mutual alignment at a specific radius since it's compelled to migrate inward at a considerably faster pace, preventing it from fully shedding its out-of-plane angular momentum, and hence, leaving the inner disk to partial-alignment only.", '4. CONCLUSIONS & DISCUSSION': "We have presented a series of global, idealized 3D GRMHD simulations of misaligned accretion onto a black hole with different spins, ranging from a = 0 . 5 to 0 . 9375, with an initial inclination of 24 · . Comparing the misaligned simulations to an aligned accretion flow reveals that in all tilted flows, the disks evolve to a highly warped and twisted steady-state structure showing no signs of the Bardeen-Petterson alignment (Bardeen & Petterson 1975). We do, however, discover a potential aligning process for weakly-magnetized, thick accretion disks ( H/r > α ) that brings the inner disk into partial alignment with the BH. The standing shocks in the tilted flows are torquing the inner disk into the BH equatorial plane, likely through the redistribution of angular momentum via the redirection of fluid elements. For instance, while analyzing the evolution of inclination of a comoving fluid [cf. Equations (8) and (9)], we observe that the decrease in inclination at the shock location is far greater due to azimuthal motion than radial motion, suggesting the orbital mixing of fluid elements of different orientations of angular momentum. This importance of azimuthal mixing is also hinted at Liska et al. (2019) and Sorathia et al. (2013a), as Sorathia et al. (2013a) postulated that for alignment to occur azimuthal mixing of angular momentum is required-the second term in equation (9) may quantify such a mechanism. Overall, the presence of spiral shocks results in an even angular momentum distribution (or, less non-axisymmetric) in the azimuthal direction.  \nHowever, it is also possible that the partial alignment we observe is a result of the oscillatory nature of bending waves, as suggested by Ivanov & Illarionov (1997) and Zhuravlev et al. (2014). Given that the nonlinear features arise close to the BH, this decrease in inclination  \ncould result from oscillatory behavior associated with standing, nonlinear bending waves.  \nOur study reveals that the strength of shocks increases with the BH spin, in line with findings of White et al. (2019). This relationship offers an explanation for the absence of shocks observed by Teixeira et al. (2014) for a BH spin of a = 0 . 1. Our findings that shocks partially align the inner flow is consistent with the steeper alignment rate observed in high-spin cases. White et al. (2019) demonstrated that shocks become more pronounced as the initial misalignment of the disk increases for the same BH spin. Combining this finding with our shock-induced partial alignment, it is conceivable that in severely misaligned flows, the inner flow may undergo complete alignment with the spin of the BH. In future research, we intend to develop models to quantify the strength of shocks as a function of varying initial tilt and BH spins, enabling us to better understand the process of alignment between a disk and its black hole.  \nWhile our simulations exhibit general characteristics consistent with previous GRMHD studies of hot, tilted flows (Fragile et al. 2007; Fragile & Blaes 2008; Generozov et al. 2013; Teixeira et al. 2014), and to a certain extent in warped thin disks (Kaaz et al. 2023), we did not observe definitive evidence of significant global precession in our system. In the case of high BH spin ( a ≊ 0 . 94), simulated over an extensive period of t ≈ 25 , 000 r g /c , our warped disk exhibited only a modest precession of ≈ 10 · for the radial shell of 15 to 50 r g , suggesting a global precession timescale of ≈ 10 6 r g /c ≈ 0 . 5( M/M ⊙ ) s . This timescale is consistent with the predictions derived from the formulation provided by (Fragile et al. 2007, Equation (43)), where we took inner and outer radius of the evolved disk to be, r i ∼ 15 r g , r 0 ∼ 50 r g as fiducial parameters 4 , and ζ ∼ -0 . 83 (derived from Σ = Σ i ( r/r i ) -ζ ) as estimated in our study. This observation suggests the possibility of the tilted disk precessing as a rigid body under the influence of Lense-Thirring torque (Lense & Thirring 1918). For lower spins, however, the simulation duration was insufficient to capture significant precession, thus limiting our ability to draw meaningful conclusions.  \nThe results of our toy model confirm the arguments of Ivanov & Illarionov (1997) and their close agreement with the local basis vector method suggests that nodes of the shocks form predominantly at ± r ˆ n , i.e. at the location where the adjacent warped annuli are least sepa-  \nted. Additionally, the simplicity of our model enables us to quantify the effects of tilt and twist. where we discover that the twist predominantly governs the geometry of the shocks, but their location can change if the disk is significantly misaligned. Moreover, while global precession was not prominently observed in our simulations, it is conceivable that, akin to the findings of Fragile & Blaes (2008), spiral shocks might have precessed on a similar timescale if global precession had been present. Considering the pivotal role of twist in determining the locations of maximum crowding, it is plausible that our toy model could predict the location of shocks in scenarios involving global precession of the disk.  \nOur model does, however, have limitations. First, it cannot predict the geometry of the shock from the initial misalignment of the torus. Complete knowledge of tilt and twist is needed to predict the spatial structure of shocks. One can avoid performing the simulations by predicting the tilt and twist profile using the linear theory, but as shown by White et al. (2019), the problem is highly non-linear. Secondly, because we use the radial profiles of tilt and twist, we could not predict the vertical structure of the shocks. To determine this we would need to completely track the full parameter variations, such as density, pressure, and velocity across the shock front. Despite having these limitations, our model shows great potential for locating shocks and can be applied to observations of warped disks. While it is true that finding the tilt and especially twist of real warped disks is a highly complex task, one could take the jet inclination as an average tilt and use this as an initial condition in simulations to potentially predict where shocks will develop.  \nThe increase in inclination at the outer radii is consistent with the physics of bending waves. However, as shown by White et al. (2019), the straightforward application of the linear theory does not properly yield the observed profile of tilt. The high sensitivity of linear theory to the surface density profiles, Shakura-Sunyaev α viscosity, and the number of assumptions that are violated in GRMHD simulations make the predictions of linear theory inapplicable. An uncharted territory being explored by Zhuravlev & Ivanov (2011) is deriving the governing equations for each component of angular momentum and assessing the contribution of the torque produced due to frame-dragging of the spacetime, pressure forces, Reynolds stress, and Maxwell stress. Approaching the problem in the said manner may also help in determining the source of shocks-induced partial alignment of the inner disk and the dominance of V ϕ ∂ I ∂ϕ in governing d I /dt .  \nThe timescale analysis shows why the inner disk never fully aligns. We define an alignment timescale and observed that t align > t inflow , suggesting that the gas does not have enough time to dissipate its out-of-plane angular momentum. While this is a simple, dimensionalitybased estimate, it explains our results quite well and we believe it to be a good order of magnitude check. Also, as one would expect from the 'alignment timescale', our definition is moot in the outer regions of the disk where the flow is significantly more misaligned.  \nFinally, standing shocks can accelerate the plasma, and can be the sites of particle acceleration and nonthermal electron distribution in collisionless accretion flows like for Sgr A ∗ (Fragile & Blaes 2008). We suggest that a radiative 3D GRMHD simulation of a misaligned accretion disk rotating around a fast-spinning BH can help us to identify the signatures of standing shocks in spectra. It may also help us to quantify the difference between the spectrum of aligned and misaligned flows. Furthermore, it will allow us to investigate the argu-  \nnts of Zhuravlev & Ivanov (2011), who claimed that strong irradiation of the outer parts of the disk by the inner parts could change the spectra of the disk. To effectively probe the effects of shocks, we also propose a parameter survey similar to this, incorporating passive scalars as tracers to investigate the fluid trajectories across the shocks front.", '5. ACKNOWLEDGMENTS': "We thank M. C. Begelman, C. J. White, K. Long, C. Echibur'u -Trujillo, and N. Scepi for stimulating discussions related to this work. We thank the anonymous reviewers for their insightful comments and suggestions, which have significantly improved our research. This work was supported in part by National Science Foundation award AST-2034306, the NASA Astrophysics Theory Program grant 80NSSC20K0527, and by an Alfred P. Sloan Research Fellowship (JD).  \nSoftware: Matplotlib (Hunter 2007), NumPy (Harris et al. 2020)", 'A. PROPERTIES OF a = 0 . 5 and 0 . 75 TILTED FLOWS': "The purpose of this appendix is to describe the properties of tilted flows for spins a = 0 . 5 and 0 . 75. In the section 3.1, we describe the warping of the disk and the emergence of non-axisymmetric features in vertically-integrated profiles of ⟨ ˙ q ⟩ ρ, θ and ⟨ l ⟩ ρ, θ for spin a ≊ 0 . 94. These profiles reveal the presence of standing shocks. Similar characteristics are observed in tilted flows with spins a = 0 . 5 and 0 . 75, as depicted in Figure 8.  \nThe warping of the disk is evident from poloidal slices of density-magnetization profiles, ρ -σ . As expected, the low-density funnel region, aligned with areas of high magnetization ( σ > 1), stands perpendicular to the disk midplane. Comparing σ profiles of tilted flows with spins a = 0 . 5, 0 . 75, and ≊ 0 . 94 (Figure 1), we note a reduction in the regions with strong magnetization as the BH's spin increases. This phenomenon is attributed to more substantial vertical motions, notably noticeable in the high-spin case, demonstrated in Figure 2, which exhibits pronounced warping.  \nThe vertically averaged snapshots of the heating rate ( ˙ q ) reveal the presence of standing shocks. Notably, the degree of non-axisymmetry intensifies as the BH's spin increases, suggesting that the strength of the shocks and their impact on disk partial-alignment and angular momentum transport escalates with higher BH spin.", 'B. TILTED FRAME ANGULAR COORDINATES': "To assess parameters within the tilted frame, which aligns with the disk, we utilize angular coordinates ( θ ' , ϕ ' ) linked to grid coordinates ( θ , ϕ ) through Equation (7) in White et al. (2019). For thoroughness, we present the relationship here:  \nθ ' = cos -1 (cos I cos θ +sin I cos T sin θ cos ϕ +sin I sin T sin θ sin ϕ ) , ϕ ' = tan -1 ( -sin T sin θ cos ϕ +cos T sin θ sin ϕ, -sin I cos θ +cos I cos T sin θ cos ϕ +cos I sin T sin θ sin ϕ ) , θ = cos -1 ( -sin I sin θ ' cos ϕ ' +cos I cos θ ' ) , ϕ = tan -1 (sin I sin T cos θ ' +cos I sin T sin θ ' cos ϕ ' +cos T sin θ ' sin ϕ ' , sin I cos T cos θ ' +cos I cos T sin θ ' cos ϕ ' -sin T sin θ ' sin ϕ ' ) (B1)  \nFigure 8. Similar to Figure 1, the above plot shows the properties of the tilted flows for BH spin, a = 0 . 5 and 0 . 75. Warping of the disk and standing shocks are apparent from poloidal slices of ρ -σ and density-weighted, vertically averaged ⟨ ˙ q ⟩ ρ,θ and ⟨ l ⟩ ρ,θ profiles. The strength of the shocks is shown to be proportional to the BH spin. However, unlike shocks strength, the strongly magnetized region shrinks with increasing spin due to strong vertical motions that result from enhanced warping.  \n<!-- image -->", 'C. RATE OF CHANGE OF INCLINATION FOR a = 0 . 5 and 0 . 75 TILTED FLOWS': 'This section delves into the rate of change of inclination for BH spin values of a = 0 . 5 and 0 . 75. Figure 9 shows that the rate is largely negative in the inner flow, especially for a = 0 . 75. This is consistent with the radial drop in the tilt profile seen in Figure 2. The scattered black points in the plot correspond to centroid positions of vertically-integrated heating rate, which we used as proxy for the determining the location of the shocks front. For the a = 0 . 75 scenario, we observe a sound agreement between the positions of shocks and regions exhibiting negative inclination rates. This suggests that the shocks affect the inclination of disk, particularly compelling the inner flow to align with the spin of the BH.  \nHowever, for the a = 0 . 5 scenario, the situation is less straightforward. This is because the shocks are relatively weak, suggesting that their influence is subtler and less discernible through our simple analysis of d I /dt . Additionally, the observed tilt reduction for the low spin is less than 2.5 · , possibly due to errors related to angular momentum projections.  \nFigure 10 presents the resulting calculations using the said methods for the high spin case. As evident from the figure, there are noticeable discrepancies in the measured values of d I dt , particularly in the vicinity of the BH. The differences between Methods 2 and 3 are straightforward to interpret: Method 2 considers only the change in the tilt of a fluid parcel as it moves with the mean fluid velocity, whereas Method 3 incorporates both mean and turbulent fluid  \n<!-- image -->  \nFigure 9. Similar to Figure 6, the above plot shows the change in rate of inclination for the spins of BH a = 0 . 5 (left) and 0 . 75 (right). The good agreement between the black points describing the location of shocks and the regions of negative rate for a = 0 . 75 case elucidates the impact of shocks in aligning the inner disk. However, because the shocks are weak in the low spin case, the decline we observe in the tilt of the inner disk may be due to numerical errors.  \n<!-- image -->', 'D. DIFFERENT METHODS TO COMPUTE RATE OF CHANGE OF INCLINATION': 'In this section, we evaluate whether the rate of change of inclination measured in the Eulerian frame at steady-state varies depending on the definition employed. The following table outlines three distinct methods, all based on the same fundamental principle, d I dt = V · ∇I .  \nTable 1. Comparison of different methods for calculating d I dt . Method 1 uses vertically-averaged velocities and inclination, Method 2 uses mean fluid velocity and three-dimensional inclination gradients, and Method 3 incorporates both mean and turbulent fluid velocities and three-dimensional inclination gradients.  \nFigure 10. The above plot shows the rate of change inclination ( d I dt ) for the BH spin a ≊ 0 . 94, calculated using three different methods. While these methods showed minor differences, the overall features of the plot remained consistent, supporting the reliability of our findings despite the potential numerical errors.  \n<!-- image -->  \nvelocities. The deviation between Method 1 and other two methods arises primarily due to the different definitions of inclination we chose. Despite these variations, the rates measured by all methods display similar non-axisymmetric structures and regions of decreasing inclination, aligning consistently with the non-axisymmetric structure in the vertically-integrated heating rate. We selected Method 1 for our primary analysis because it defines inclination as the angle between the angular momentum vector of the vertically-averaged disk and the BH spin axis.', 'E. GEOMETRY OF THE SPIRAL SHOCK IN THE DISK-ALIGNED FRAME': "To substantiate the robustness of the m = 2 standing shocks, we quantified the heating rate within the disk-aligned frame. Initially, we transformed the metric tensor to the tilted frame using the relation g µ ' ν ' = ∂x α ∂x µ ' ∂x ζ ∂x ν ' g αζ , where g αζ represents the metric defined in Kerr-Schild coordinates. The Jacobian ( ∂x α ∂x µ ' ) for this transformation is detailed in White et al. (2019). Subsequently, we computed the area element in the new coordinates as √ -g ' , and calculated the angular coordinates θ ' and ϕ ' following the equations outlined in Appendix B. These calculations revealed that at specific radii, the angular coordinates were neither uniformly distributed nor covered the complete polar ([0 , π ]) and azimuthal ([0 , 2 π ]) ranges.  \nTo address this, we constructed a uniform angular grid and interpolated the relevant physical quantities onto it. The density-weighted vertical integration of the heating rate was then performed using the equation:  \n⟨ ˙ q ⟩ ρ,θ ' ( r, ϕ ' ) = ∫ ∫ ˙ q ( r, θ ' , ϕ ' ) ρ ( r, θ ' , ϕ ' ) √ -g ' dθ ' dt ∫ ∫ ρ ( r, θ ' , ϕ ' ) √ -g ' dθ ' dt (E2)  \nFigure 11 illustrates these calculations for the high-spin tilted case and corroborates the persistence of the standing shocks. Notably, the nodes of the shocks predominantly align with the Y Disk = 0 plane, corresponding to a local precession angle of ϕ ' = 0 and π . This alignment explains the observed influence of the twist angle on the spatial configuration of the shocks in the original BH frame.", 'REFERENCES': "Abramowicz, M. A., & Fragile, P. C. 2013, Living Reviews  \nin Relativity, 16, 1, doi: 10.12942/lrr-2013-1  \nBardeen, J. M., & Petterson, J. A. 1975, ApJL, 195, L65,  \ndoi: 10.1086/181711  \nChatterjee, K., Liska, M., Tchekhovskoy, A., & Markoff, S. 2023, Misaligned magnetized accretion flows onto spinning black holes: magneto-spin alignment, outflow power and intermittent jets. https://arxiv.org/abs/2311.00432  \nChatterjee, K., Younsi, Z., Liska, M., et al. 2020, Monthly Notices of the Royal Astronomical Society, 499, 362, doi: 10.1093/mnras/staa2718  \nFigure 11. The above plot shows the density-weighted, vertically averaged heating rate ( < ˙ q > ρ,θ ' ) in the tilted frame. Standing shocks predominantly form near ϕ ' = 0 and π , highlighting the influence of twist in determining the shock locations in the un-tilted BH frame.  \n<!-- image -->  \nCui, Y., Hada, K., Kawashima, T., et al. 2023, Nature, 621,  \n711, doi: 10.1038/s41586-023-06479-6 Dexter, J., & Fragile, P. C. 2011, The Astrophysical Journal, 730 -. 2013, Monthly Notices of the Royal Astronomical Society, 432, 2252, doi: 10.1093/mnras/stt583 Dexter, J., Scepi, N., & Begelman, M. C. 2021, The Astrophysical Journal Letters, 919, L20, doi: 10.3847/2041-8213/ac2608 Drewes, N. C., & Nixon, C. J. 2021, The Astrophysical Journal, 922, 243, doi: 10.3847/1538-4357/ac2609 Esin, A. A., McClintock, J. E., & Narayan, R. 1997, The Astrophysical Journal, 489, 865, doi: 10.1086/304829 Fishbone, L. G., & Moncrief, V. 1976, ApJ, 207, 962, doi: 10.1086/154565 Fragile, P. C., & Blaes, O. M. 2008, ApJ, 687, 757, doi: 10.1086/591936 Fragile, P. C., Blaes, O. M., Anninos, P., & Salmonson, J. D. 2007, The Astrophysical Journal, 668, 417, doi: 10.1086/521092 Gammie, C. F., McKinney, J. C., & Toth, G. 2003, The Astrophysical Journal, 589, 444, doi: 10.1086/374594 Generozov, A., Blaes, O., Fragile, P. C., & Henisey, K. B. 2013, The Astrophysical Journal, 780, 81, doi: 10.1088/0004-637x/780/1/81 Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2 Hawley, J. F., & Krolik, J. H. 2018, The Astrophysical Journal, 866, 5, doi: 10.3847/1538-4357/aadf90 Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55  \nIvanov, P. B., & Illarionov, A. F. 1997, MNRAS, 285, 394,  \ndoi: 10.1093/mnras/285.2.394 Kaaz, N., Liska, M. T. P., Jacquemin-Ide, J., et al. 2023, The Astrophysical Journal, 955, 72, doi: 10.3847/1538-4357/ace051 Kawaguchi, K., Kyutoku, K., Nakano, H., et al. 2015, Phys. Rev. D, 92, 024014, doi: 10.1103/PhysRevD.92.024014 King, A., & Nixon, C. 2018, The Astrophysical Journal, 857, L7, doi: 10.3847/2041-8213/aab8f9 King, A. R., Lubow, S. H., Ogilvie, G. I., & Pringle, J. E. 2005, Monthly Notices of the Royal Astronomical Society, 363, 49, doi: 10.1111/j.1365-2966.2005.09378.x Lense, J., & Thirring, H. 1918, Physikalische Zeitschrift, 19, 156 Liska, M., Hesp, C., Tchekhovskoy, A., et al. 2018, MNRAS, 474, L81, doi: 10.1093/mnrasl/slx174 Liska, M., Hesp, C., Tchekhovskoy, A., et al. 2021, Monthly Notices of the Royal Astronomical Society, 507, 983, doi: 10.1093/mnras/staa099 Liska, M., Tchekhovskoy, A., Ingram, A., & van der Klis, M. 2019, Monthly Notices of the Royal Astronomical Society, 487, 550, doi: 10.1093/mnras/stz834 Lubow, S. H., Ogilvie, G. I., & Pringle, J. E. 2002, Monthly Notices of the Royal Astronomical Society, 337, 706, doi: 10.1046/j.1365-8711.2002.05949.x McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2013, Science, 339, 49, doi: 10.1126/science.1230811 Mewes, V., Font, J. A., Galeazzi, F., Montero, P. J., & Stergioulas, N. 2016, Physical Review D, 93, 064055 Narayan, R., & McClintock, J. E. 2008, NewAR, 51, 733, doi: 10.1016/j.newar.2008.03.002"}