Radiation Transport Two-Temperature GRMHD Simulations of Warped Accretion Disks

In many black hole systems, the accretion disk is expected to be misaligned with respect to the black hole spin axis. If the scale height of the disk is much smaller than the misalignment angle, the spin of the black hole can tear the disk into multiple, independently precessing `sub-disks'. This is most likely to happen during outbursts in black hole X-Ray binaries (BHXRBs) and in active galactic nuclei (AGN) accreting above a few percent of the Eddington limit, because the disk becomes razor-thin. Disk tearing has the potential to explain variability phenomena including quasi-periodic oscillations (QPOs) in BHXRBs and changing-look phenomena in AGN. Here, we present the first radiative two-temperature GRMHD simulation of a strongly tilted ($65^{\circ}$) accretion disk around a $M_{BH}=10M_{\odot}$ black hole, which tears and precesses. This leads to luminosity swings between a few percent and $50 \%$ of the Eddington limit on sub-viscous timescales. Surprisingly, even where the disk is radiation pressure dominated, the accretion disk is thermally stable over $t \gtrsim 21,000 r_g/c$. This suggests warps play an important role in stabilizing the disk against thermal collapse. The disk forms two nozzle shocks perpendicular to the line of nodes where the scale height of the disk decreases $10$-fold and the electron temperature reaches $T_e \sim 10^8-10^9 K$. In addition, optically thin gas crossing the tear between the inner and outer disk gets heated to $T_e \sim 10^8 K$. This suggests that warped disks may emit a Comptonized spectrum that deviates substantially from idealized models.


INTRODUCTION
Accretion disks in black hole X-ray binaries (BHXRBs) and active galactic nuclei (AGN) are typically described by axisymmetric, time-independent models. A key factor that determines the most suitable model is the accretion rate with respect to the Eddington limit. At the Eddington limit, radiation pressure rivals gravity and plays a dynamically important role. The most well-known models are the geometrically thin standard accretion disk model (Shakura & Sunyaev 1973), which describes sources accreting at an appreciable fraction of their Eddington limit (L 0.01 − 0.1 L Edd ), and the geometrically thick advection dominated accretion flow model (Narayan & Yi 1994;Narayan et al. 2003), which describes highly sub-Eddington sources (L 0.01L Edd ). Slim disk models, in which the advection of trapped photons dominates, work well near and above the Eddington limit (e.g. Abramowicz et al. 1978). While these models can provide a reasonable description of the multi-wavelength emission, they -by construction -do not address the variability of their light curves and spectra. Such variability can encode important information about the structure of the accretion disk and corona (e.g. Kara et al. 2019;Uttley et al. 2014).
For example, recent observational evidence indicates that accretion disks in changing-look AGN can undergo drastic changes in luminosity and spectral shape over the course of less than a year (e.g. MacLeod et al. 2016;Yang et al. 2018a) or during quasi-periodic eruptions (QPEs) over the course of hours (e.g. Yang et al. 2018b). In these AGN, the luminosity can change by up to two orders of magnitude while the spectrum significantly hardens or softens. This suggest a dramatic change in the accretion flow on timescales which are too short to be explained by the viscous timescales inferred from the Maxwell stresses created by magnetorotational instability (MRI, Balbus & Hawley 1991 induced turbulence (e.g. Lawrence 2018;Dexter & Begelman 2019). Various physical mechanisms were proposed to reconcile theory with observations, including instabilities in the accretion disk (Sniegowska et al. 2020), magnetically elevated accretion disks (Begelman & Pringle 2007), reprocessing of pointsource radiation by the outer accretion disk (e.g. Clavel et al. 1992), disk tearing ) and an orbiting compact object (e.g. King 2020; Arcodia et al. 2021).
Similar to their changing-look AGN counterparts, the power spectra of BHXRBs display a wide range of mysterious features ranging from broad-spectrum variability to narrow spectral peaks known as quasi-periodic oscillations (QPOs). Such features can encode unique information about the structure and dynamics of the accretion disk in addition to the spin and mass of the black hole. QPOs are typically divided into low-and high-frequency QPOs, which can sometimes be observed together (e.g. Ingram & Motta 2019). Various physical mechanisms have been proposed to explain QPOs. These include geometric effects such as precession of a misaligned disk around the black hole spin axis (e.g. Stella & Vietri 1998;Ingram et al. 2009) and intrinsic effects such as parametric resonances (e.g. Kluzniak & Abramowicz 2002;Rezzolla et al. 2003;Abramowicz & Kluźniak 2004) and disco-seismic modes (e.g. Kato 2004;Dewberry et al. 2020a,b). All of these models make simplifying assumptions about the underlying physics. Most significantly, they typically do not include magnetized turbulence (though see Dewberry et al. 2020b;Wagoner & Tandon 2021), which can potentially dampen oscillatory modes. Thus, general relativistic magnetohydrodynamic (GRMHD) simulations, which can simulate accretion from first-principles, are very attractive for addressing the origin of QPOs.
Interestingly, numerous GRMHD simulations have demonstrated that a misaligned disk can precess (e.g. Fragile & Anninos 2005;Fragile et al. 2007;Morales Teixeira et al. 2014;Liska et al. 2018aLiska et al. , 2019aWhite et al. 2019). However, this only happens when the disk is unrealistically small and, in fact, Liska et al. (2018a) demonstrated that the precession period of more realistically sized disks (r f ew × 100r g ) becomes too long to explain any observed QPOs. This problem can be solved if a smaller precessing disk tears of from a larger non precessing disk due to the differential frame dragging (Lense & Thirring 1918) of space-time by a spinning black hole. This process, called disk tearing, was observed in various SPH simulations (e.g. Nixon et al. 2012a,b;Nealon et al. 2015; and recently confirmed by GRMHD simulations (Liska et al. 2019b(Liska et al. , 2022b. Disk tearing, interstingly, also leads to rapid burst of accretion possibly explaining large luminosity and spectral swings observed in accreting black holes (e.g. Nixon & Salvesen 2013;. Musoke et al. (2022) recently addressed the origin of both low and high-frequency QPOs using GRMHD simulations where a 65 • misaligned accretion disk tore apart. Based on these simulations, Musoke et al. (2022) argued that lowfrequency and high-frequency QPOs can be explained by disk tearing. In their model, a disk tears off at a nearconstant radius (r ∼ 10 1 r g ) and precesses for 10 periods before falling into the black hole. While precessing at ν = 2.5 Hz (for M BH = 10M ), the disk emits a periodically modulated light curve. Musoke et al. (2022) also found a prominent radial epicyclic oscillations at ν ∼ 56 Hz in the inner disk, which they argue can potentially explain some highfrequency QPOs. The geometric origin of low-frequency QPOs is supported by recent observations, which suggest that that iron line centroid frequency moves from blue-toredshifted with the same phase as the low-frequency QPO (Ingram & van der Klis 2014;Ingram et al. 2016). On the other hand, observations have been able to constrain neither an intrinsic nor geometric origin of the much rarer highfrequency QPOs (e.g. Ingram & Motta 2019).
Thus, the results of Musoke et al. (2022) suggest that disk tearing is a very promising mechanism to explain the multiwavelength variability in many accreting black holes. However, all GRMHD simulations that found tearing (Liska et al. 2019b(Liska et al. , 2022b have relied on a cooling function (Noble et al. 2009) to artificially set the temperature and scale height of the disk. In reality, the temperature of the disk is determined by complex physics involving viscous heating, radiative cooling and the advection of energy. Such physics can modify the radial and vertical structure of the disk, such as its scale height and vertical temperature profile. This makes it extremely challenging to benchmark such GRMHD models against multi-wavelength observations because they might not capture the important physics that sets the disk thermodynamic state. For example, in radiation pressure supported disks imbalance between the local dissipation rate and cooling rate can lead to thermal and viscous instabilities that are not captured by a cooling function (e.g. Shakura & Sunyaev 1976;Sadowski 2016;Jiang et al. 2019;Liska et al. 2022a). In addition, when the inner disk precesses, it will beam its radiation field periodically towards the outer disk. This radiation field will scatter off the outer disk, and the radiation back-reaction force can cause additional warping and dissipation (e.g. Pringle 1997;Wijers & Pringle 1999).
Here, we present the first radiation-transport twotemperature GRMHD simulation of an accretion disk tilted by 65 • relative to a rapidly spinning M BH = 10M black hole of dimensionless spin of a = 0.9375 accreting at 35% the Eddington luminosity. We show that the disk undergoes tearing. The accretion disk does not contain any large scale vertical magnetic flux to begin with and does not launch a jet: this makes the simulation applicable to the soft(-intermediate) states of BHXRBs (see Fender et al. 2004). In Sections 2 and 3 we describe the numerical setup and initial conditions, before presenting and discussing our results in Sections 4 and 5, and concluding in Section 6.

NUMERICAL SETUP
In this work we utilize the radiative version of our GPU accelerated GRMHD code H-AMR (Liska et al. 2018a(Liska et al. , 2022b. H-AMR evolves the ideal GRMHD equations and radiation transfer equations in addition to the electron and ion entropies. We use a (modified) spherical grid in Kerr-Shild foliation with coordinates r, θ and φ. Spatial reconstruction of primitive variables is performed using a 1-dimensional PPM method, which guarantees second order convergence in 3 dimensions. Magnetic fields are evolved on a staggered grid as described in Gardiner & Stone (2005); White et al. (2016). Inversion of conserved to primitive variables is performed using a 2-dimensional Newton-Raphson routine (Noble et al. 2006) or Aitken Acceleration scheme (Newman & Hamlin 2014) for the energy equation and a 1-dimensional Newton-Raphson routine for the entropy equation. The energy based inversion is attempted first, and, if it fails, the entropy based inversion is used as a backup. This dual energy formulation is now standard in many GRMHD codes .
The radiative transfer equations are closed with a twomoment M1-closure (Levermore 1984) whose specific implementation is described in Sądowski et al. (2013); McKinney et al. (2013). We additionally evolve (see e.g. Noble et al. 2009;Ressler et al. 2015) the electron and ion entropy tracers (κ e,i = p e,i /ρ γe,i ) and include Coulomb collisions (Stepney 1983). We use adiabatic indices γ e = γ i = 5/3 for the electrons and ions (and γ = 5/3 for the gas). This is appropriate for the ion and electron temperatures in the accretion disk. For a typical 'peak' temperature of T e ∼ 10 9 K and T i ∼ 10 10 K we find that γ e ∼ 1.57, γ i ∼ 1.66, γ g ∼ 1.66. We use a reconnection heating model (Rowan et al. 2017) to distribute the dissipative heating between ions and electrons. This leads to roughly 20 − 40% of the dissipation going into the electrons with the remainder heating the ions.
We include Planck-averaged bound-free, free-free, and cyclo-Synchrotron absorption (κ abs ), emission (κ em ) and electron scattering (κ es ) opacities as given in McKinney et al. (2017) for solar abundances X = 0.7, Y = 0.28, and Z = 0.02. In addition, we account for Comptonization through a local blackbody approximation (Sadowski & Narayan 2015). This works reasonably well when emission, absorption and scattering are localized, but becomes inaccurate when the radiation field is anisotropic and/or the energy spectrum deviates from a blackbody. Namely, the radiation temperature used to calculate the absorption opacity is approximated as a blackbody with T r = (Ê rad /a) 0.25 withÊ rad the fluid frame radiation energy density and a the radiation density constant. If the radiation deviates from a blackbody this approximation can severely under or overestimate the radiation temperature. We assume that only electrons can absorb and emit radiation and set the ion opacity to κ i = 0. This is a reasonable approximation in the accretion disk where bound-free processes, which involve energy exchange between ions/electrons/radiation, only become important when the density is high enough for Coulomb collisions to equilibrate the ion and electron tem-perature. The limitations of the M1 closure relevant to our work are further addressed in Section 5.
To resolve the accretion disk in our radiative GRMHD models we use a base grid resolution of N r × N θ × N φ = 840 × 432 × 288. We then use beyond r 4r g up to 3 layers of adaptive mesh refinement (AMR) to progressively increase this resolution in the disk from the base resolution at r ∼ 4r g to N r × N θ × N φ = 6720 × 2304 × 4096 at r ∼ 10r g . (In addition, we use static mesh refinement Liska et al. 2022b to reduce the φ-resolution to N φ = 16 at the polar axis.) This is exactly half of the resolution used in Liska et al. (2022b) and Musoke et al. (2022). We place the radial boundaries at R in = 1.13r g ≈ 0.84r H and R out = 10 5 r g , which ensures that both boundaries are causally disconnected from the flow; here r H = 1 + (1 − a 2 ) 1/2 r g is the event horizon radius. The local adaptive time-stepping routine in H-AMR increases the timestep of each mesh block independent of the refinement level in factors of 2 (up to a factor of 16) based on the local Courant condition (Courant et al. 1928). This increases the effective speed of the simulations several-fold while reducing the numerical noise .
Since ideal GRMHD is unable to describe the physical processes responsible for injecting gas in the jet funnel (such as pair creation) we use density floors in the drift frame (Ressler et al. 2017) to ensure that ρc 2 ≥ p mag /12.5. However, due to the absence of any poloidal magnetic flux, no jets are launched in this work and thus the corresponding density floors are never activated. Instead, the ambient medium is floored at the density ρ = 10 −7 × r −2 and internal energy u tot = 10 −9 × r −2.5 , where u tot = u gas + u rad is the sum of the gas, u gas , and radiation, u rad , internal energies.

PHYSICAL SETUP
To study the effect of disk tearing in a highly warped accretion disk we evolve our model in two stages. In both stages we assume a rapidly spinning black hole of mass M BH = 10M and dimensionless spin a = 0.9375. The initial radial (r) and vertical (z) density profile of the disk is ρ(r, z) ∝ r −1 exp (−z 2 /2h 2 ). The disk extends from an inner radius, r in = 6r g , to the outer radius, r out = 76r g at a constant scale height h/r = 0.02. The covariant magnetic vector potential is A θ ∝ (ρ − 0.0005)r 2 and normalized to give an approximately uniform β = p b /p tot ∼ 7 in the initial conditions. Here p b is the magnetic pressure and p tot is the sum of the electron/ion pressures, p e,i , and the radiation pressure, p rad . This rather low value of β ∼ 7 was chosen such that the grid is able to resolve both the toroidal and poloidal MRI wavelengths for a physically meaningful runtime. A parameter exploration considering higher β values will become possible in the coming 5 − 10 years after the advent of (post-)exascale GPU clusters. The whole setup is subsequently rotated by 65 degrees with respect to the black hole spin axis, which Figure 1. The first demonstration of disk tearing in a radiation pressure supported accretion disk, as seen through iso-countour 3D renderings of density at three different times. (Left panel) The disk at t = 28, 955rg/c forms a rigid body that is warped, but shows no differential precession. (Middle panel) The disk tears off a precessing inner disk at t = 31, 000rg that precesses for two periods between t = 35, 000rg/c and t = 37, 000rg/c. While precessing, the inner disk aligns with the BH, shrinking in size until it fully disappears. (Right panel) The disk forms a new tear at r ∼ 25rg/c around t ∼ 45, 000rg/c, and the tearing process repeats.
itself is aligned with polar axis of the grid. Physical equivalence between tilting the disk (e.g. Liska et al. 2018a;White et al. 2019) and tilting the black hole (e.g. Fragile & Anninos 2005;Fragile et al. 2007) is demonstrated in the Appendix of Liska et al. (2018a).
In the first stage (model T65), which is described in Liska et al. (2022b), we evolve the disk in ideal GRMHD for ∆t = 150, 000r g /c with a preset scale height of h/r = 0.02 set by a cooling function (Noble et al. 2009). The disk undergoes a large tearing event around t ∼ 45, 000r g /c, which is analysed in Musoke et al. (2022). During this tearing event a sub-disk of size ∆r ∼ 10 − 15r g tears off and precesses for ∼ 6 periods between t ∼ 45, 000r g /c and t ∼ 80, 000r g /c. While precessing, the inner disk exhibits a very prominent radial epicyclic oscillation at ν = 56 Hz for a 10 solar mass BH. Another tearing event between t ∼ 120, 000r g /c and t ∼ 145, 000r g /c produces a similar oscillation albeit at a slightly higher frequency of ν = 69 Hz. These frequencies are consistent with observed high frequency QPOs (e.g. Ingram & Motta 2019).
In the second stage (model RADT65), which is the focus of this work, we restart this simulation in full radiative GRMHD at t = 29, 889r g /c as described above. The resolution is downscaled by a factor of 2 to reduce the cost of the simulation and make it amenable to the pre-exascale GPU cluster OLCF Summit on which this simulation has been run. This is well before the occurrence of any large tearing event in stage 1. We assume thermodynamic equilibrium between gas and radiation by setting T r = T e = T i in the initial conditions and run this model (RADT65) for ∆t ∼ 21, 000r g /c until t ∼ 49, 787rg/c. The density scaling is set such that the average accretion rate corresponds to ap-proximatelyṀ/Ṁ Edd ∼ 0.35 with M Edd = 1 ηNT L Edd /c 2 the Eddington accretion rate, L Edd = 4πGMBH ckes the Eddington luminos-ity, and η NT = 0.179 the Novikov & Thorne (1973) efficiency. According to Piran et al. (2015) this choice of parameters is roughly consistent with a thin disk solution (Sadowski 2011) of scale height of h/r = 0.02 − 0.06 between r = r isco and r = 100r g . 4. RESULTS

Tearing Process
The disk in RADT65 undergoes a tearing event between t = 31, 000r g /c and t = 40, 000r g /c, as visualised using a 3D density isocontour rendering at three different times in Fig. 1). Figure  During this tearing event the disk precesses for 2 periods between 36, 000 < t < 38, 000r g /c (Fig. 2e) before it aligns with the spin axis of the black hole (Fig. 2f). Alignment happens through warp driven dissipation of misaligned angular momentum on the accretion time of the inner disk, which differs from the Bardeen-Petterson alignment (Bardeen & Petterson 1975) mechanism that manifests itself on much shorter timescales during which the structure and density of the disk remains in a steady state (see discussion in Liska et al. 2018b). In addition to the first tearing event, which will be the focus of this article, a much larger tearing event starts after t 47, 000r g /c that does not finish before the end of this simulation. Due to the large computational expense as- Figure 2. Panels (a, b) Space-time diagrams of density, ρ, and warp amplitude, ψ. The warp amplitude peaks where the disk is about to tear or has already torn. Since dissipation is enhanced in warps, maxima in warp amplitude correspond to minima in density. The main tearing event is denoted by black dotted lines. Panels (c, d) Space-time diagrams of the electron, Te, and ion, Ti, temperatures. The ion and electron temperature correlate with the warp amplitude and anti-correlate with the density. Close to the black hole and in the tear the ions are significantly hotter than the electrons since Coulomb collisions are unable to equilibrate the ion and eletron temperatures. Panels (e, f) Space-time diagrams of the tilt, T , and precession, P, angle of the disk. Around t ∼ 36, 500rg/c a disk tears off and precesses for 2 periods, during which it gradually aligns with the black hole spin axis. Another tearing event occurs after t ∼ 47, 000rg/c. sociated with extending the duration of RADT65, we leave the full analysis of this longer tearing event to future work.
Whether a disk is able to tear depends on how the internal torques in the disk react to differential Lense-Thirring precession (e.g. Nixon & King 2012). The Lense-Thirring precession rate of a point particle around a rapidly spinning black hole follows a steep radial dependence of ν LT ∝ a/r 3 . However, an accretion disk will precess with a single frequency set by the integrated LT torque in the disk if the viscous torques are able to counteract the differential precession rate (e.g. Fragile & Anninos 2005;Fragile et al. 2007). Naively, one expects that the warp amplitude decreases monotonically with the increasing distance from the black hole. This is (roughly) the case in the inner and outer disk of RADT65 for t 32, 000r g /c.
As Fig. 2(a),(b) shows, the warp amplitude and surface density in RADT65 are anti-correlated for the entire runtime. At r 20r g the warp amplitude is rather high and a cavity of low density gas forms at both early and late times. In this cavity the plasma thermally decouples into hot electrons with T e ∼ 10 8 −10 9 K and extremely hot ions with T i ∼ 10 9 −10 10 K.
The strong anti-correlation between warp amplitude and density suggests that accretion is driven by dissipation in warps (e.g. Papaloizou & Pringle 1983;Ogilvie 1999;Nelson & Papaloizou 2000). In a companion paper (Kaaz et al. 2022) we demonstrate that warp-driven dissipation is responsible for the bulk of the dissipation and that magneto-rotational instability driven turbulence (e.g. Balbus & Hawley 1991 plays only a very subdominant role. However, after t 32, 000r g /c the warp amplitude is not a monotonically decreasing function of radius. Instead, the warp amplitude forms local minima and maxima, and concentric rings of higher density gas centered at the (local) minima of the warp amplitude form. It is known that when the warp amplitude exceeds a critical value, it can become unstable (Doǧan et al. 2018;Dogan & Nixon 2020). When this happens, the viscous torques in the warp drop to zero and the disk tears. In RADT65 the tearing process starts at r ∼ 19r g around 33, 000r g /c and the disk fully detaches at r ∼ 10r g around t ∼ 36, 500r g /c. We see later that during this process warp-driven dissipation leads to a temporary burst in the BH accretion rate and luminosity around t ∼ 36, 500r g /c Figure 3. A volume rendering of the density during the main tearing event at t = 37, 026rg/c. The red plane cuts the disk at the φ = 0 surface. Streamers from the outer disk (green) deposit low angular momentum gas onto the inner disk (blue). This causes the inner accretion disk to shrink in size (Kaaz et al. 2022).
As the gas from the inner disk falls into the black hole it transports angular momentum outwards. Angular momentum conservation dictates that the radius of the inner disk should increase. This process, called viscous spreading, affects all finite size accretion disks that are not resupplied externally with gas (e.g. Liska et al. 2018a;Porth et al. 2019). However, the inner disk in RADT65 does not spread viscously, but, in fact, decreases in size. This discrepancy could potentially be explained by cancellation of misaligned angular momentum (e.g. Nixon & King 2012;Hawley & Krolik 2018). Namely, when gas from the outer disk falls onto the inner disk the misaligned components of angular momentum cancel out leading to a net decrease in the inner disk radius. In addition, it is possible that the inner disk transfers some of its excess angular momentum to the outer disk. In a companion paper (Kaaz et al. 2022) we quantify the relative contributions of the angular momentum cancellation and outward angular momentum transport to the evolution of the tearing radius.

Dissipation and Radiative Signatures
To better understand the structure of the disk during the tearing event at t = 37, 026r g /c we show a volume rendering of the density in Figure 3 and a vertical projection of several quantities in Figure 4. In this projection the line of nodes, where the disk crosses the equatorial plane, is aligned with the horizontal x-axis. The projected quantities include the surface density Σ = π 0 ρ √ g θθ dθ (panel a), the scattering optical depth τ es ∼ 0.34Σ (panel b), electron and ion temperature T e,i = π 0 pe,i √ g θθ dθ Σ (panels c and d), density scale height L edd (panel f), gas entropy per unit mass (panel e), and radiation entropy per Integration is performed in a spherical coordinate system aligned with the local angular momentum vector of the disk. An animation of Figure 3 and panels a, c and f of Figure 4 can be found on our YouTube channel (movie).
Interestingly, the electron and ion temperatures do not decline smoothly with distance from the black hole (Fig. 4b,c) as predicted by idealized models of thin accretion disks Novikov & Thorne (1973). At the tearing radius of 10r g , where the warp amplitude reaches a maximum, the electron temperature reaches T e 10 8 K and the plasma becomes optically thin. Here streamers of gas get thrown onto highly eccentric orbits and subsequently rain down on the inner disk (Fig. 3). The rapid rise in temperature in the tear suggests that gas crossing the tear is subject to additional dissipation as it undergoes a rapid orbital plane change (see also Nixon & Salvesen 2013;). If we assume the gas has a temperature of T i = 0K before crossing the tear at r ∼ 7r g , and the gas dissipates = 10 −2 of its orbital kinetic energy U kin ≈ 1 2r , we find that dissipation can easily heat the gas to T i ∼ 10 10 K. In reality, can be much higher if radiative cooling is efficient.
The scale height of the disk exhibits a prominent m = 2 azimuthal oscillation (Fig. 4e). This can also be observed in Figure 5 where we show a θ − φ slice at r = 25r g of the density, electron temperature, and gas entropy. At the line of nodes the scale height reaches a maximum of h/r ∼ 0.1 while perpendicular to the line of nodes the scale height drops to h/r ∼ 0.005. Around this nozzle-like 'compression', the gas reaches a temperature ranging from T e ∼ 2×10 7 K at r ∼ 20r g to T e 10 8 K within r 5r g . At select times T e can even reach T e ∼ 10 9 K closer to the black hole (see YouTube movie). The radial and azimuthal fluctuations in the temperature lead to a highly non-axisymmetric emission pattern (Fig. 4f). At this specific snapshot (t = 37, 026r g /c) the effective bolometric luminosity varies azimuthally between λ ∼ 0.15L edd and λ ∼ 0.5L edd at r = 200r g .
The specific gas entropy does not increase in the nozzle (Fig. 4g), suggesting heating of gas in the nozzle is adiabatic. This might seem surprising since nozzles in tidal disruption events (TDEs) are typically associated with shock heating and steep gradients in the gas entropy (e.g. Rees 1988;Kochanek 1994;Andalman et al. 2022). However, if the radiative cooling timescale is very short, shocks do not The gas entropy entropy, sgas, does not increase substantially when gas crosses the scaleheight compression along the vertical y-axis, but the radiation entropy, s rad , does. This suggests the shock heating is radiatively efficient. automatically lead to an increase of the gas entropy. Instead, they will lead to an increase in the radiation entropy. This increase in radiation entropy is visible along the vertical yaxis in Fig. 4(g). In a companion paper (Kaaz et al. 2022), we demonstrate that a shock forms in the nozzle that dissipates roughly ∼ 1.8% of the orbital energy each time the gas passes through it.

Outflows and Variability
To understand the variability associated with disk tearing and the energetics of the outflows we plot in Figure 6 the mass accretion rate (panel a), radiative emission rate (panel b) and (radiative) outflow efficiencies (panel c). The mass accretion rate measured at the event horizon (Ṁ = 2π 0 π 0 ρu r √ −gdθdφ with g the metric determinant and u µ the fluid 4-velocity) exhibits a factor ∼ 40 variation in time, ranging fromṀ = 2 × 10 −2Ṁ Edd toṀ =Ṁ Edd (defined in Sec. 3). However, the luminosity (L = 2π 0 π 0 R r t √ −gdθdφ with R r t the r-t component of the radiation stress energy tensor) only exhibits a factor ∼ 5 variation and peaks at L ∼ 0.5L edd .
The outflow efficiencies with respect to the BH mass accretion rate are defined at the event horizon as η wind =Ṁ −Ė |Ṁ|t , η rad = the energy accretion rate with T r t the stress energy tensor, and, |Ṁ| t is a running average ofṀ over a time interval of 500r g /c. Unless stated otherwise, we measureṀ andĖ at the event horizon, L at the event horizon for η adv and at r = 200r g for η rad . As Fig. 6(c) shows, when the accretion rate reaches the Eddington limit the advective efficiency, η adv , is of comparable magnitude to the radiative efficiency, η rad (Fig. 6d). This implies that a significant fraction of the emitted radiation will fall into the black hole before it reaches the observer suggesting the disk might have entered the slim disk regime at the peak of luminosity. Nevertheless, the timeaveraged radiative efficiency reaches |η rad | = 14.7%, which is only slightly below the canonical Novikov & Thorne (1973) efficiency of η NT = 17.8%.
We do not observe any evidence for a radiation or magnetic pressure driven wind. While the outflow efficiency measured at the event horizon exceeds η wind 15% this drops to η wind 0.1% when measured at r = 200r g , which demonstrates that both the magnetic pressure (e.g. Liska et al. 2018b) and radiation pressure (e.g. Kitaki et al. 2021) are insufficient to accelerate the wind to escape velocities. This suggests that poloidal magnetic fields are a key ingredient to accelerate sub-relativistic outflows, even when the accretion rate approaches the Eddington limit. Figure 5. A θ − φ slice through the disk at r = 15rg in tilted coordinates with the ascending node at φ = 1/2π and the descending node at φ = 3/2π. The pink and black lines correspond to the τes = 1 and τ b f = 1 surfaces respectively. (Upper panel) The gas density, ρ, peaks in the 'nozzle', which is located perpendicular to the local line of nodes at φ = 0 and φ = π. (Middle panel) The electron temperature, Te, increases by a factor ∼ 5 in the nozzle, suggesting the radiative emission will be hardened. (Lower panel) The gas entropy per unit mass, κgas, does not increase in the nozzle, suggesting the rise in gas temperature is adiabatic (Sec. 4.2)

Radial Structure
To better understand the structure of the accretion disk we show in Figure 7 time-averaged (between t = [43, 000, 44, 000]r g /c) radial profiles of the density,ρ (panel a), luminosity, L (panel b), and radiation, ion, electron, and magnetic pressure, p rad , p i , p e , and p b , respectively (panel c). During this time interval a new tear starts to form around r ∼ 30r g , which manifests itself as a drop in density and increase in the radiative flux. While the magnetic and radiation pressure are in approximate equipartition within r 20r g , the The wind efficiency, η wind , radiative efficiency, η rad , advective radiative efficiency, η adv , and Novikov & Thorne (1973) efficiency, ηNT , and total efficiency, ηtot , as defined in Sec. 4.3. When the accretion rate peaks, advection of radiation becomes important and suppresses the total luminosity of the system. (Panel e) The thermal scale height of the disk, θt , stays stable or increases throughout the disk, which suggests the disk is thermally stable (Sec. 4.5). accretion disk becomes fully radiation pressure dominated between r = [20, 40]r g . The sum of the ion and electron pressure is a factor ∼ 10 − 100 lower than the radiation pressure, suggesting that gas pressure does not play a dynamically important role.
We have already mentioned in Section 4.1 that warp driven dissipation plays an important role in RADT65. To illustrate how much warp-driven dissipation accelerates the inflow of gas in RADT65, we translate the radial infall speed of the gas into an effective viscosity, α eff = gives the y−weighted average of quantity x. In Fig. 7(c), we find α eff ∼ 1 − 5 × 10 1 throughout the disk. The Maxwell stresses α m = b r b φ ρ ptot ρ ∼ 10 −2 − 10 −1 induced by magnetized turbulence cannot account for this effective viscosity, which suggests that other physical processes such as shocks (e.g. Fragile & Blaes 2008) drive accretion in warped accretion disks.

Thermal Stability
Interestingly, Liska et al. (2022a) found that an aligned version of the accretion disk in RADT65 collapsed within t 5, 000r g /c to an exceedingly thin slab that could not be resolved numerically. This runaway cooling (see also Sadowski 2016;Fragile et al. 2018;Jiang et al. 2019;Mishra et al. 2022) is a manifestation of the thermal instability (Shakura & Sunyaev 1976). However, we do not observe any evidence for the thermal instability in RADT65. In fact, as Fig. 6(e) demonstrates, the thermal scale height, θ t = cs ρ v φ ρ , remains stable or even increases depending on the radius.
This unexpected result could potentially be explained within r 20r g by magnetic fields, which can stabilize a thermally unstable accretion disk if they reach equipartition with p b ∼ p rad (e.g. Begelman & Pringle 2007;Sadowski 2016;Jiang et al. 2019;Liska et al. 2022a;Mishra et al. 2022). However, magnetic fields will not be able to stabilize the disk between r = [20, 40]r g where the disk is strongly radiation pressure dominated with p rad /p b ∼ 10 1 − 10 2 .
We found in Sec. 4.3 that advection of radiation becomes significant at the event horizon. If the radiation diffusion timescale becomes longer than the accretion timescale, photon trapping might stabilize the disk against runaway cooling (e.g. Abramowicz et al. 1978). To test if photon trapping can explain the apparent thermal stability of the disk we plot the radiation diffusion time, t diff ∼ h v rad ∼ 3hτ es and the accretion time, t acc = r vr ρ , in Fig. 7(e). Here we use for v rad the optically thick radiation diffusion speed v rad ∼ 1 3τes c (e.g. Fig. 7(e), while the radiation diffusion and accretion timescales are comparable in the inner disk, the radiation diffusion timescale becomes much shorter than the accretion timescale in the outer disk. This suggests that photon trapping will not be able to stabilize the outer disk against runaway cooling or heating.

McKinney et al. 2013). As is evident from
Instead, we propose in Section 5 that nozzle shock driven accretion disks are inherently thermally stable.

Applicability to Spectral States
The accretion disk in RADT65 does not launch any relativistic jet. This suggests RADT65 is applicable to the highsoft and potentially soft-intermediate spectral states. In these spectral states the spectrum is dominated by thermal blackbody radiation and no radio jets are detected (though see Russell et al. 2020 for detection of compact jets in the infrared). The cumulative luminosity of the disk increases until 35rg, partially driven by enhanced dissipation in the forming tear between r = 20rg and r = 30rg. (Panel c) The disk is radiation pressure, p rad , dominated between r = [15, 40]rg and magnetic pressure, p b , dominated everywhere else. The ion, pi), and electron, pe, pressures are very weak compared to either the radiation or magnetic pressure and thus do not play an dynamically important role. (Panel d) The effective viscosity, αeff, is much larger than the Maxwell, αm, stress. This suggests that dissipation is not driven by magnetized turbulence induced by magnetic stresses. Instead, we demonstrate in a companion paper that shocks are driving accretion (Kaaz et al. 2022). (Panel e) The photon diffusion timescale, tdiff, is shorter than the accretion timescale, tacc, except very close to the black hole. This suggest advection of radiation internal energy is only dynamically important very close to the black hole.
Nevertheless, the emission profile of RADT65 presents a radical departure from the Novikov & Thorne (1973) model of a thin accretion disk. Instead of an axisymmetric powerlaw emission profile, emission in warped accretion disks can become hardened around the nozzle shock or at the tearing radius (see also Nixon & Salvesen 2013;. Here the temperature of plasma can reach T e ∼ 10 8 − 10 9 K, which might lead to Comptonization of cold accretion disk photons. This could potentially contribute to the high-energy power law emissions tail observed in the spectra of soft state BHXRBS (e.g. Remillard & McClintock 2006) in addition to the soft X-Ray excess in the spectra of high-luminosity AGN (e.g. Gierliński & Done 2004. Note that part of this high-energy emission can come from within the ISCO as suggested by recent GRMHD simulations (Zhu et al. 2012) and semi-analytical models (Hankla et al. 2022). Future general relativistic ray-tracing (GRRT) calculations will need to quantify how much the spectrum deviates from a blackbody spectrum and in which wavelength bands this optically thin shock-heated gas can be detected. The energy spectrum of the accelerated electrons will depend on the microphysics of the shocks, which will need to be addressed by particle-incell simulations.
In the hard-intermediate state (HIMS) jetted ejections imply that the accretion disks might be saturated by vertical magnetic flux (e.g. McKinney et al. 2012;Tchekhovskoy et al. 2011;Chatterjee et al. 2020). When an accretion disk is saturated by magnetic flux it enters the magnetically arrested disk (MAD) regime (Narayan et al. 2003) in which magnetic pressure overcomes gravity and accretion proceeds through non-axisymmetric instabilities. Recent two-temperature radiative GRMHD simulations (Liska et al. 2022a) have demonstrated that magnetic flux saturation can lead to 'magnetic' truncation during which the inner disk decouples into a magnetic pressure supported 'corona' of hot electrons (T e 5 × 10 8 K) and extremely hot ions (T i 10 10 K). This happens within the magnetic truncation radius, which is determined by the amount and location of the excess poloidal magnetic flux in the system. It is conceivable that the combined effects of magnetic reconnection in current sheets and shocks in warps can heat the plasma to temperatures far above what was observed in this work or in Liska et al. (2022a). Future work will need to address if, where and how such misaligned truncated disks tear, and, if their spectral and timing signatures are consistent with sources in the HIMS.

Quasi-periodic Variability
This work and Musoke et al. (2022) have demonstrated that disk tearing can lead to short luminosity outbursts of length ∆t ∼ 1000 − 3000r g /c every ∆t tear ∼ 10, 000 − 50, 000r g /c. Coincident with a tearing event the disk precesses on a timescale of ∆t prec ∼ 1000 − 6000r g /c suggesting tearing events might be accompanied by a sinusoidial modulation of the light curve. In parallel, the mass-weighted radius of the precessing inner disk oscillates on a timescale of ∆t osc ∼ 150 − 400r g /c suggesting another sinusoidal modulation of the lightcurve albeit at a higher frequency (Fig. 8 and Musoke et al. 2022). We summarize the timescales of these vari- −gdθdφ , of the inner disk (after tearing) in RADT65 (blue) features a radial epicyclic oscillation consistent with the radial epicyclic frequency of an oscillating ring of gas at r = 6.9rg (orange). Such oscillations may be associated with HFQPOs detected in BHXRB lightcurves. ability phenomena in Table 1 for a BHXRB (M BH = 10M ), a low-mass AGN (M BH = 10 6 M ), and a high-mass AGN (M BH = 10 9 M ). Please note that these timescales might vary by at least an order of magnitude depending the location of the tearing radius, which we expect will depend on e.g. the accretion rate, black hole spin, misalignment angle, and amount of large scale magnetic flux. Table 1 suggests that the timescales associated with disk tearing are consistent with variety of astrophysical phenomena (see also Nixon & Salvesen 2013;. The tearing events themselves might be associated with the heartbeat modes (on 1 − 100 s timescales) detected in the BHXRBs GRS 1915+105 and IGR J17091-3624 (e.g. Belloni et al. 2000;Neilsen et al. 2011;Weng et al. 2018;Altamirano et al. 2011) in addition to QPEs detected in various AGN (e.g. Miniutti et al. 2019;Giustini et al. 2020;Arcodia et al. 2021;Miniutti et al. 2022). If the tearing events in some spectral states are spaced further apart, they could potentially also explain some changing-look phenomena observed in AGN with masses of 10 7 − 10 9 M that occur on timescales of years to decades (e.g. Sniegowska et al. 2020). In all of these cases the luminosity fluctuates by a factor 10 0 − 10 2 suggesting some significant flaring event in the accretion disk reminiscent of a tearing event. The precession of the disk that follows each tearing event can potentially explain low frequency QPOs detected in numerous BHXRB with ν ∼ 10 −1 − 10 1 Hz (e.g. Ingram & Motta 2019). Precession is also consistent with the 44 day LFQPO detected in the 2 × 10 8 M AGN KIC 9650712 (Smith et al. 2018) and the 3.8H LFQPO detected in the 10 5 − 10 6 M AGN GJ1231+1106 (Lin et al. 2014). In addition, the oscillation of the mass-weighted radius might be associated with HFQ-POs of ν ∼ 10 1 − 10 2 Hz in BHXRBs and some of the short duration QPOs detected in various AGN. These include the ∼ 1H QPO in the lightcurve of the 10 6 − 10 7 M AGN RE J1034+396 (Gierlinski et al. 2008;Czerny et al. 2016), the Table 1. Time period between large tearing events, ∆ttear, disk precession period, ∆tprec, and period of radial epicyclic oscillations, ∆tosc, for various black hole masses. These timescales can be associated with various variability phenomena in BHXRBs and AGN such as QPEs, QPOs, and heartbeat oscillations (Sec. 5).
The heated gas at the tearing radius suggests the QPO emission will be substantially hardened through Comptonization. This is consistent with frequency-resolved spectroscopy of BHXRBs which has demonstrated that the Type-C QPO waveform is dominated by hardened emission (Sobolewska &Życki 2006) and recent observations of AGN that suggest an increase in the soft-X-ray excess during QPEs (Giustini et al. 2020). Interestingly, our simulations suggest that tearing events are accompanied by a low-frequency sinusoidal signal caused by precession of the inner disk and a high-frequency sinusoidial signal caused by a radial epicyclic oscillation at the tearing radius. This can potentially be tested using observational data for both BHXRBs and AGN.

Comparison to Semi-analytical Models
Semi-analytical models typically divide warped accretion disks into two categories (e.g. Papaloizou & Pringle 1983) depending on the scale height of the disk relative to the α−viscosity parameter. In the first category bending waves excited by oscillating radial and vertical pressure gradients in the warp are damped and the warp is propagated through viscous diffusion only. This happens when the scaleheight is smaller than the viscosity parameter, e.g. h/r < α. In this regime instability criteria were derived that predict the tearing of accretion disks (e.g. Ogilvie & Dubus 2001;Doǧan et al. 2018;Dogan & Nixon 2020). In the second category the warp is propagated through bending waves that travel at a fraction of the sound speed (Papaloizou & Lin 1995). These bending waves are excited by pressure gradients in a warp. Past work has made substantial progress in describing warp propagation in the bending wave regime (e.g. Ogilvie & Latter 2013;Paardekooper & Ogilvie 2018;Deng et al. 2020;Ogilvie 2022), and recently also addressed disk tearing in the bending wave regime using SPH simulations and linearised fluid equations (Drewes & Nixon 2021).
Interestingly, Deng & Ogilvie (2022) predicted using the affine model of Ogilvie (2018) that oscillations in the verti- Figure 9. A scatterplot of the effective viscosity, αeff, with respect to the thermal scaleheight, θt at r = 5rg in green, r = 10rg in black, and r = 20rg in blue. We show power-law fits to the data as solid lines with fitting function αeff ∝ (h/r) −γ . The effective viscosity increases as the disk becomes thinner, which suggests warp-driven dissipation becomes more important for thinner disks. We argue in Section 5.4 that warp-driven dissipation can stabilize the disk against runaway thermal collapse (e.g. Shakura & Sunyaev 1976). cal pressure gradient lead to an azimuthal m = 2 oscillation in the scale height of the disk similar to what is observed in our work. These authors also demonstrated that a parametric instability can induce turbulence and create ring-like structures in warped accretion disks. Since h/r ∼ α m (Fig. 7d) and the warp is highly non-linear it is unclear in which regime of warp propagation our disk falls (see also Sorathia et al. 2013;Hawley & Krolik 2018). We do note though that there is some evidence for wave-like structures in the disk (Fig. 10), which suggests the disk might be in the bending wave regime of warp propagation. Nevertheless, future work will need to isolate wave-like from diffusive (misaligned) angular momentum transport to make more direct comparison to semianalytical models possible.

Thermal and Viscous Stability of Warped Disks
In Section 4.5 we demonstrated that warped, radiation pressure supported, accretion disks remain thermally stable for much longer than similar non-warped, radiation pressure supported, accretion disks considered in previous work (Liska et al. 2022a). We concluded that neither magnetic pressure nor advection of radiation internal energy can explain this thermal stability. Here we propose an alternative explanation.
For a disk to remain thermally stable the derivative of the dissipation rate, Q vis , with respect to the temperature needs to be smaller than the derivative of the cooling rate, Q rad , with respect to the temperature such that Qvis dT < Q rad dT . For a radiation pressure supported α−disk Q vis ∝ T r φ H ∝ αp tot H ∝ αp 2 tot ∝ αT 8 and Q rad ∝ T 4 (Shakura & Sunyaev 1973;Meier 2012). This suggests that with a constant α-viscosity radiation pressure supported accretion disks are thermally unstable. However, if the dissipation rate follows Q vis ∝ T ζ with ζ 4 due to warping, the disk will remain thermally stable. Figure 10. A transverse slice through the density along the line of nodes at 50rg. This snapshot is taken beyond the 'official' duration of our simulation (Sec. 3), just after the onset of the next tearing event. Outwards propagating wave-like structures (visible as vertical 'stripes') form in the outer disk suggesting that wave-like angular momentum transport might be present. The vertical extent of the outer disk is artificially inflated because the frame is sliced at an angle with respect to the disk.
To determine ζ in our simulation we present in Figure 9 a scatter plot of the effective viscosity, α eff , with respect to the thermal scale height, θ t . We only use data after 32, 000r g /c to allow the disk to adjust after seeding it with radiation and exclude data where the tilt angle is smaller than T < 60 • . We visually fit the data to find that α eff ∝ H −γ with γ ∼ 1.25 − 1.75. From the relation H ∝ p tot ∝ T 4 we find that ζ = 8 − 4γ ∼ 1 − 3. This suggests that even though our accretion disk is radiation pressure dominated at r = 20r g , warp-driven dissipation can potentially prevent the disk from becoming thermally unstable.
Following a similar argument, one can show that the warpdriven dissipation also avoids the viscous  instability of radiation pressure dominated disks. If a disk is subject to the viscous inability, it will disintegrate into rings of high density gas (e.g. Mishra et al. 2016). Even though high-density rings do form in our simulation (Figs. 2 and 4), we conclude that this cannot be a manifestation of the viscous instability since similar rings form in the nonradiative analogue to RADT65 presented by Musoke et al. (2022). This suggests that the ring formation in RADT65 is warp-induced.

Radiation Warping
When the inner disk precesses it will periodically beam radiation towards the outer disk. This can potentially warp the outer disk and induce additional dissipation (similar to Pringle 1997;Wijers & Pringle 1999). In a similar context Liska et al. (2019b) has demonstrated that when a precessing jet collides with the outer accretion disk it can inject energy and angular momentum into the disk, leading to an increase of particle orbits in the outer disk. Since the energy efficiency of the radiative outflows in RADT65 is non-negligible (|η rad | ∼ 14.7%) we expect that radiation feedback can have a profound effect on the structure and dynamics of an accretion disk.
Radiative feedback is an interesting avenue of future research that is hard to address in this article due to limitations of our M1 radiation closure (Levermore 1984). M1 works well in modeling radiative cooling and energy transport in optically thick regions, and in optically thin regions for a single radiation bundle. However, M1 is unable to properly model crossing radiation bundles. The reason for this is that M1 treats radiation as a highly collisional fluid and is thus unable to model crossing radiation bundles in optically thin media. Instead the energy and momentum of the two radiation bundles will be summed together. While the M1 closure conserves the total amount of energy and momentum in the photon-fluid, it can lead to a nonphysical redistribution of energy and angular momentum.
Despite these limitations we demonstrate in Figure 11 that the structure of the inner accretion disk in RADT65 differs significantly after seeding it with radiation. In Figure 11 the radiation streamlines get deflected by the inner disk presumably leading to an increase in the warp amplitude compared to the start of the simulation. In this snapshot the tearing process has not started yet, which suggests that radiation warping can be important even in the absence of disk tearing. While these preliminary results are encouraging, more advanced numerical simulations will be necessary to reliably address radiation warping.

Role of Magnetic Fields
In this work we have seeded the accretion disk with a relatively strong toroidal magnetic field of β = pgas+p rad p b ∼ 7 (Sec. 3) in order to being able to resolve the poloidal MRI components and reach numerical convergence. Alternatively, we could have threaded the accretion disk with poloidal magnetic flux loops and chosen a higher plasma-β in the initial conditions as we have done in the past (e.g. Liska et al. 2019b). However, this leads to the formation of powerful jets and would make these simulations more applicable to the hard-intermediate spectral state (e.g. Liska et al. 2022a).
In an accompanying paper (Kaaz et al. 2022) we demonstrate that the turbulent Maxwell stresses seeded by MRI turbulence fall short by 1-3 orders of magnitude in explaining the measured accretion rate. Instead, we demonstrate that accretion is driven by shock dissipation in the nozzle. This naively suggests that magnetic fields do not play an important Figure 11. Radiation field warps the inner disk, as seen in a transverse slice of density through the x − z plane of our initial conditions (left) and after ∆t = 2, 000rg/c (right) with radiation streamlines shown in red. Even though our M1 radiation scheme is unable to properly model multi-beam radiation (Sec. 5.5), these results suggest that radiation warping might be an important effect. role in the disk dynamics and thus that the presented results would be insensitive to the chosen plasma-β.
However, magnetic fields might still play an important role in communicating the differential Lense-Thirring torque throughout the accretion disk. This can have important consequences for the tearing process and the formation of the nozzle shock. In fact, SPH models (e.g. Nixon et al. 2012a,b;Nealon et al. 2015;), which do not include magnetic fields, find that a thin accretion disk tears apart into multiple deferentially precessing rings instead of a radially extended subdisk observed in our GRMHD simulations. In addition, SPH models do not find a nozzle shock. These arguments suggest that magnetic fields might play an important, though indirect role, in both the tearing process and the formation of the nozzle shock, which needs to be addressed in future work.

Numerical Convergence
RADT65 is the largest radiative 2T GRMHD simulation to date, featuring on average over 3 × 10 9 cells, and requiring over 6000 V100 GPUs on OLCF Summit. Nevertheless, to reduce the numerical cost of this work and bring it within the realm of possibilities, we needed to downscale the resolution by a factor of 2 compared to similar work presented in Liska et al. (2022b) and Musoke et al. (2022) that did not take radiation into account. Despite this downscaling, the number of cells per MRI-wavelength, Q i = λ i MRI /N i , still exceeds (Q θ , Q φ ) (20, 125) and the number of cells per disk scaleheight, Z, exceeds Z 14.
This suggests, based on more controlled convergence studies (e.g. Shiokawa et al. 2012;Porth et al. 2019), that both the disk and the MRI turbulence are reasonably well resolved and close to achieving convergence. In fact, the vast majority of GRMHD simulations performed to-date (e.g. Porth et al. 2019) have similar or even slightly inferior Q and Z parameters. In addition, physical quantities like the inflow speed of the gas and the evolution of the tearing radius, are very similar to the non-radiative GRMHD models in Liska et al. (2022b); Musoke et al. (2022) suggesting the results presented in this article are robust and no significant qualitative changes are expected at higher resolutions.
Despite these arguments a full convergence study will be an interesting use-case for future generations of GPU clusters. Namely, doubling the resolution is currently computationally infeasible. Doing so, would 16−fold the cost of the presented simulations to over 4 million node hours on OLCF Summit. This would exceed the size of our entire INCITE allocation by more than an order of magnitude. On the other hand, we already know based on previous convergence studies (e.g. Shiokawa et al. 2012;Porth et al. 2019) that reducing the resolution by a factor of 2 would under-resolve both the poloidal MRI components and the scaleheight of the disk, and thus fail to reach convergence.

CONCLUSION
In this article we have presented the first radiative twotemperature GRMHD simulation of an 65 • misaligned disk accreting atṀ ∼ 0.35Ṁ edd . Similar to the idealized GRMHD models presented in Liska et al. (2019bLiska et al. ( , 2022b; Musoke et al. (2022) the radiation pressure supported accretion disk tears apart and forms a precessing disk. During the tearing process the mass accretion rate increases by a factor ∼ 10−40 and the luminosity increases by a factor ∼ 5, almost exceeding the Eddington limit at the peak of the outburst. As proposed in , the timescales and amplitudes of these luminosity swings are roughly consistent with QPEs in AGN (e.g. Miniutti et al. 2019;Giustini et al. 2020;Arcodia et al. 2021;Miniutti et al. 2022) and several heartbeat modes in BHXRBs (e.g. Belloni et al. 2000;Neilsen et al. 2011;Altamirano et al. 2011;Weng et al. 2018). Following a tearing event, the inner disk precesses for several periods with a frequency that is consistent with low-frequency QPOs (e.g. Ingram & Motta 2019). While precessing, the inner disk also exhibits a radial epicyclic oscillation of its center of mass whose frequency is consistent with high frequency QPOs (e.g. Ingram & Motta 2019). Future long-duration GRMHD simulations combined with radiative transfer calculations will need to test if disk tearing indeed produces QPOlike variability. This will require the tearing radius to remain stable over many cycles of tearing.
Warping of the disk forms two nozzle shocks directed perpendicular to the line of nodes. At the nozzle shock, the scale height of the disk decreases by a factor of ∼ 5 and the temperature increases up to T e ∼ 10 8 − 10 9 K and T i ∼ 10 9 − 10 10 K. In a companion paper (Kaaz et al. 2022), we demonstrate that dissipation in the nozzle shock leads to a much shorter accretion time scale than predicted by α-viscosity based models (e.g. Shakura & Sunyaev 1973) seeded by magnetized turbulence (e.g. Balbus & Hawley 1991). This can potentially explain the timescales associated with some changing look phenomena in AGN including QPEs. In addition, we find during tearing events that gas crossing the tear heats up to T e ∼ 10 8 K and T i ∼ 10 10 K. This is caused by an increase in the dissipation rate as the gas undergoes a rapid orbital plane change (see also Nixon & Salvesen 2013;. These results imply that the emission profile from warped accretion disks will deviate substantially from the Novikov & Thorne (1973) model of a geometrically thin accretion disk. This can, pending future radiative transfer calculations, have farreaching consequences for spectral fitting and black hole spin measurements (e.g. Zhang et al. 1997;Kulkarni et al. 2011).
While radiation pressure supported accretion disks in GRMHD simulations typically collapse (e.g. Sadowski 2016;Fragile et al. 2018;Jiang et al. 2019;Mishra et al. 2020;Liska et al. 2022a) into an infinitely thin slab due to the thermal instability (Shakura & Sunyaev 1976), this did not happen in our simulation. Here we argue that the warp and associated nozzle shock stabilize the disk against thermal collapse (Sec. 5). These shocks differ from spiral density wave induced shocks (e.g. Arzamasskiy & Rafikov 2018). In a companion paper (Kaaz et al. 2022), we demonstrate that nozzle shock dissipation drives accretion in geometrically thin warped accretion disks within at least r 50r g . This upends the standard paradigm that invokes magneto-rotational instability induced turbulence for the dissipation of orbital kinetic energy (e.g. Balbus & Hawley 1991. Future work will need to address the range of spectral states, accretion rates, and misalignment angles where nozzle shock dissipation becomes important. Past work almost exclusively invoked magnetic reconnection in coronae to explain high energy (non-)thermal emission from accretion disks (e.g. Sironi & Spitkovsky 2014;Beloborodov 2017;Ripperda et al. 2021). However, this work suggests that heated plasma formed in the nozzle and at the tearing radius can potentially be another source of high energy emission. This is especially appealing to explain nonthermal emission in the high-soft and soft-intermediate states of BHXRBs. Those states do not launch any large-scale jets and are thus unlikely to be saturated by magnetic flux (e.g. McKinney et al. 2012), which might be a key ingredient to form a corona where the prime source of dissipation comes from magnetic reconnection (e.g. Liska et al. 2022a;Ripperda et al. 2021). Test-particle methods (e.g. Bacchini et al. 2019) informed by first principles particlein-cell simulations are an interesting avenue for future research to quantify the contribution of magnetic reconnection and warp-driven shocks to the acceleration of non-thermal particles and corona-like emission.