This site uses cookies. By continuing to use this site you agree to our use of cookies. To find out more, see our Privacy and Cookies policy. Close this notification
Skip to content

The following article is Open access

Nozzle Shocks, Disk Tearing, and Streamers Drive Rapid Accretion in 3D GRMHD Simulations of Warped Thin Disks

, , , , , , and

Published 2023 September 20 © 2023. The Author(s). Published by the American Astronomical Society.
, , Citation Nicholas Kaaz et al 2023 ApJ 955 72 DOI 10.3847/1538-4357/ace051

Download Article PDF
DownloadArticle ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

0004-637X/955/1/72

Abstract

The angular momentum of gas feeding a black hole (BH) may be misaligned with respect to the BH spin, resulting in a tilted accretion disk. Rotation of the BH drags the surrounding spacetime, manifesting as Lense–Thirring torques that lead to disk precession and warping. We study these processes by simulating a thin (H/r = 0.02), highly tilted (${ \mathcal T }=65^\circ $) accretion disk around a rapidly rotating (a = 0.9375) BH at extremely high resolutions, which we performed using the general-relativistic magnetohydrodynamic code H-AMR. The disk becomes significantly warped and continuously tears into two individually precessing subdisks. We find that mass accretion rates far exceed the standard α-viscosity expectations. We identify two novel dissipation mechanisms specific to warped disks that are the main drivers of accretion, distinct from the local turbulent stresses that are usually thought to drive accretion. In particular, we identify extreme scale height oscillations that occur twice an orbit throughout our disk. When the scale height compresses, "nozzle" shocks form, dissipating orbital energy and driving accretion. Separate from this phenomenon, there is also extreme dissipation at the location of the tear. This leads to the formation of low-angular momentum "streamers" that rain down onto the inner subdisk, shocking it. The addition of low-angular momentum gas to the inner subdisk causes it to rapidly accrete, even when it is transiently aligned with the BH spin and thus unwarped. These mechanisms, if general, significantly modify the standard accretion paradigm. Additionally, they may drive structural changes on much shorter timescales than expected in α-disks, potentially explaining some of the extreme variability observed in active galactic nuclei.

Export citation and abstract BibTeX RIS

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

1. Introduction

The traditional description of an accretion disk is the axisymmetric, geometrically thin Shakura–Sunyaev model (Shakura & Sunyaev 1973; Pringle 1981 for a review; Novikov & Thorne 1973 for the relativistic treatment). While the Shakura–Sunyaev model has found widespread utility in the field, it is by no means a complete description of accretion. One of its major assumptions is that the disk angular momentum is aligned with the black hole (BH) spin. Yet, in most natural circumstances, the infalling gas that forms the disk has no prior knowledge of the BH spin orientation. Thus, many accretion disks may be at least initially misaligned. This can drastically alter the dynamics of the disk because the general-relativistic "frame-dragging" of the Kerr BH will apply Lense–Thirring (LT) torques to the disk (Lense & Thirring 1918; Misner et al. 1973). These torques induce differential precession about the BH spin vector and can lead to large-scale warps in the disk. Early analytic work found that these LT torques can align inner regions of tilted accretion disks with the BH spin (Bardeen & Petterson 1975), which has sometimes been invoked to neglect the effects of disk tilt. While recent numerical simulations have confirmed the existence of BP alignment in misaligned disks with small tilts (Nelson & Papaloizou 2000; Lodato & Pringle 2007; Perego et al. 2009; Nealon et al. 2015; Liska et al. 2019), at larger tilts the story changes dramatically. In both smoothed particle hydrodynamics (SPH; Nixon et al. 2012, 2013; Drewes & Nixon 2021; Raj et al. 2021) and general-relativistic magnetohydrodynamic (GRMHD) simulations of highly tilted disks (Liska et al. 2021, 2022; Musoke et al. 2022), the LT torques are strong enough to sometimes rupture the accretion disk, splitting it either into individually precessing annuli or into discrete subdisks. It is these highly tilted, warped and torn disks that we focus on in this paper.

There is a wealth of analytic work devoted to understanding the dynamics of warped accretion disks. This includes early work in the linearized domain of small warps (Papaloizou & Pringle 1983; Kumar & Pringle 1985, see also Pringle 1992); the fully nonlinear 1D theory that generalized the study of warped accretion disks to arbitrarily sized warps (Ogilvie 2001, 2000; Ogilvie & Latter 2013); and more recently, the more sophisticated affine model that is general to both warps and eccentricities and treats the disk as a composition of mutable fluid columns (Ogilvie 2018).

While analytic work provides a firm foundation for the understanding of warped disks, there remains only partial agreement between theory and the results of numerical simulations. In particular, disk tearing is a highly nonlinear process that results in discontinuities in the accretion flow, which are difficult to study analytically. These systems also feature anomalously high mass accretion rates that are also difficult to reconcile with theory. This was first reported in works based on SPH simulations (Nixon et al. 2012), which attributed the rapid accretion to the cancellation of misaligned angular momentum in torn regions. Rapid accretion is also found in GRMHD simulations of thin, tilted disks, which have reported effective viscosities well in excess of those expected in aligned thin disks (Liska et al. 2021).

In this work, we reveal multiple novel mechanisms that enable rapid accretion in GRMHD simulations of thin, tilted disks. In Section 2, we present the details of our simulation. In Section 3, we examine the structure of the warped accretion flow. In Section 4, we show that mass accretion occurs anomalously fast, and investigate the dissipation mechanisms that drive this rapid accretion. In Section 5, we contextualize the theoretical impact of our results, discuss their observational implications, and then summarize our findings.

2. Simulation Details

In this paper, we study a simulation of a thin, tilted accretion disk performed with the GPU-accelerated, 3D GRMHD code H-AMR (Liska et al. 2022). We work in spherical polar coordinates (r, θ and φ) and use a rapidly rotating (a = 0.9375) BH. We initialize the disk with an aspect ratio H/r = 0.02 and set the inner and outer radii to r = 6.5 rg and r = 76 rg, respectively. We initialize the velocity as circular everywhere, set the radial surface density profile to Σ ∝ r−1, and set the vertical density profile to a Gaussian profile with an FWHM equal to the local scale height of the disk. We then tilt the disk with respect to the equatorial plane by 65°. We insert into the disk a purely toroidal magnetic field described by a covariant vector potential, Aθ ∝ (ρ − 0.0005)r2, where ρ is the fluid frame gas density, which is normalized such that $\max \,\rho =1$. We normalize the magnetic field so that the average ratio of gas to magnetic pressure is β ≈ 7, such that the disk remains dominated by gas pressure. We maintain the disk thickness by using a cooling function that removes excess internal energy from the disk (e.g., Noble et al. 2009; Liska et al. 2019).

We perform the simulation on a spherical grid that is uniform in $\mathrm{log}r$ at extremely high resolutions, which are needed to resolve the turbulent motions within the disk. To achieve high resolution, we use several numerical speed-ups, including acceleration on GPUs, three levels of adaptive mesh refinement (AMR), and five levels of local adaptive time-stepping. The maximum effective resolution at ∼10 rg in the disk, where rgGM/c2 is the gravitational radius, is Nr × Nθ × Nϕ = 13440 × 4608 × 8192 cells. This resolution remains uniform within four disk scale heights to ensure that the disk structure is independent of the AMR criterion used. The three AMR levels used to achieve this resolution are added at 2, 4, and 8 rg, such that at the event horizon, the resolution is reduced to 1728 × 576 × 1024 in order to prevent the minimum time step from becoming too small, thereby speeding up the computation. To prevent the Courant condition (Courant & Hilbert 1953) in the φ direction from limiting the time step, we reduce the azimuthal resolution progressively from Nϕ = 1024 cells near the equator to 16 cells within 30° of either pole. Both inner and outer radial boundary conditions allow matter and magnetic fields to freely leave the domain. We set the inner radial boundary to be five cells within the event horizon and the outer radial boundary to be sufficiently large such that both are causally disconnected from (and thus do not affect) the flow. The polar boundary condition is transmissive, and the azimuthal boundary condition is periodic (Liska et al. 2018). We refer the reader to Liska et al. (2022) for a full description of H-AMR and Musoke et al. (2022) for an analysis of this simulation in the context of quasiperiodic oscillations in X-ray binaries (XRBs).

In a companion paper, Liska et al. (2023; which we will refer to as L23), we perform the same simulation that we have analyzed here, except we separately evolve the electron and ion entropies and use an M1 closure scheme for radiation rather than a predefined cooling function. The inclusion of radiation breaks the self-similarity of the flow; so in L23, we set the BH mass to 10 M and the Eddington ratio to ∼0.35. The results of the present work generally carry over in L23, and we refer the reader there for detailed comparisons of radiative versus cooled simulations of our disk.

3. Accretion Geometry

3.1. Main Features

Figure 1 shows a sequence of 3D renderings of the fluid frame gas density (ρ) in our simulation, separated in time by δ t ≈ 1000 rg/c. This figure depicts some of the main features of our simulation. The disk extends out to roughly ≈100 rg. As the central BH rotates, it drags spacetime with it, causing the surrounding matter to rotate as well. This so-called "frame-dragging" effect induces LT torques (Lense & Thirring 1918) onto orbiting fluid parcels. These torques cause particle orbits to precess with an angular frequency following the radial dependence ΩLTr−3. Were there no (magneto)hydrodynamic stresses acting on the disk, the disk would shred completely into independently precessing annuli. In reality, MHD stresses work to redistribute disk angular momentum. When MHD forces win, the disk maintains a warped structure, which we characterize as a series of concentric annuli with smoothly varying tilt and precession angles. However, when the LT torques induce a strong enough warp, the disk becomes unstable (Doǧan et al. 2018; Doğan & Nixon 2020), leading to a runaway increase in the amplitude of the warp that manifests as a "tear." This is seen most prominently in panels (b) and (c) of Figure 1, where the outer and inner regions of the disk break apart. These two "subdisks" precess almost independently of one another. The long-lived outer subdisk undergoes near rigid-body rotation, while the short-lived inner subdisk often precesses differentially and is typically quickly consumed by the BH. Afterwards, the outer subdisk refills the inner region, and the tearing process repeats.

Figure 1.

Figure 1. The LT torques induced by the rotation of the central BH cause the accretion disk to warp and, sometimes, tear into discrete subdisks. In each panel, we plot a 3D rendering of the fluid frame gas density, separated by δ t ≈ 1000 rg/c. Azimuthal oscillations in the scale height are apparent in the outer subdisk, where orbiting fluid parcels experience compressions and expansions twice an orbit (evidenced by the light-blue "spokes" in the outer subdisk). These oscillations are also apparent in the side-on view of the disk shown in the inset of panel (d).

Standard image High-resolution image

In total, we identified 10 tearing events over the course of our simulation (e.g., Table 1 of Musoke et al. 2022). A typical tearing cycle can last anywhere between ∼102 and ∼104 rg/c, and while the exact location of the tear varies, it usually occurs at r ≲ 10 − 20 rg. The snapshots in Figure 1 depict the transition between two tearing cycles. In Figure 1(a), the inner subdisk from the previous tearing cycle is about to be completely consumed. In Figure 1(b), the inner subdisk is replenished, and is about to tear again. In Figure 1(c), the inner subdisk tears off once again, and continues to precess in Figure 1(d).

Figure 1 shows that the azimuthal distribution of gas density exhibits a periodic structure: ρ increases and decreases twice an orbit, as evidenced by the twin light-blue high-density "spokes" in the outer subdisk. Whereas the volumetric density increases locally, the surface density does not: this is because the volumetric density increase is due to vertical (transverse to the disk) compression. This is also shown in the inset panel of Figure 1(d), where we show a side-on view of the disk. This might come as a surprise, as accretion disks are typically treated as axisymmetric, without any variation in the azimuthal direction. This nonaxisymmetry is fundamentally due to the warp, which we investigate in the following subsection.

3.2. Warp and Nonaxisymmetric Flow Structures

To better understand the nonaxisymmetric structures in our accretion disk, we begin by examining the local properties of a single annulus. For this, we use tilted coordinates r, $\theta ^{\prime} $, and $\varphi ^{\prime} $, which essentially "flatten" the warp of the disk and make analysis more convenient. We describe the transformation to tilted coordinates in the Appendix, with an accompanying visualization in Figure 10. The main features to note are as follows:

  • 1.  
    Because the transformation depends on the radial tilt and precession profiles, it is different at different radii.
  • 2.  
    The tilted vertical unit vector, $\hat{z}^{\prime} $, is co-aligned with the angular momentum of the disk at every radius.
  • 3.  
    The tilted polar coordinate, $\theta ^{\prime} $, is set such that the local midplane of the disk is at $\theta ^{\prime} =\pi /2$.
  • 4.  
    The tilted azimuthal coordinate, $\varphi ^{\prime} $, is set such that the local precession angle is at $\varphi ^{\prime} =0$.

In Figure 2, we depict the fluid frame gas density (ρ, panel (a)) and mass flux (ρ ur , panel (b)) as a function of $\theta ^{\prime} $ and $\varphi ^{\prime} $. This plot is shown at time t = 90, 471.8 rg/c and radius r = 13 rg. Since gas motion is predominantly azimuthal, the flow approximately follows the $\theta ^{\prime} \mbox{--}\varphi ^{\prime} $ plane. In Figure 2(a), we can immediately see that the disk scale height varies drastically as a function of $\varphi ^{\prime} $. Specifically, it undergoes periodic oscillation, compressing and expanding twice an orbit. The velocity streamlines (depicted in black) also follow this periodic motion, converging and diverging in phase with the oscillations of the disk scale height. The radial mass flux, depicted in Figure 2(b), also oscillates in phase with the scale height oscillations, reversing direction at the compression points. Additionally, these radial motions are approximately antisymmetric about the midplane. Fluid parcels below (above) the midplane move inward (outward) during the $\sim 0\lt \varphi ^{\prime} \lt \pi $ expansion and outward (inward) during the $\sim \pi \lt \varphi ^{\prime} \lt 2\pi $ expansion. This radial "sloshing" of the disk is much larger in magnitude than the average, inward radial mass flux associated with accretion.

Figure 2.

Figure 2. An annulus of the warped disk experiences vertical and radial oscillations twice an orbit. Panel (a): we depict the fluid frame gas density, ρ, at radius r = 13 rg and time t = 90471.8 rg/c. We have also drawn velocity streamlines of orbiting fluid parcels in black. The plot is depicted in tilted coordinates, $\varphi ^{\prime} $ and $\theta ^{\prime} $, where $\varphi ^{\prime} =0$ indicates the local precession angle of the disk and $\theta ^{\prime} =\pi /2$ indicates the local midplane of the disk. Panel (b): the same as the top panel, except we plot the radial mass flux ρ ur . The radial mass flux also exhibits oscillations twice an orbit, except they are antisymmetric about the local midplane of the disk.

Standard image High-resolution image

These oscillations can be qualitatively understood by considering how a warp impacts the hydrodynamic force balance of the disk (see also Lodato & Pringle 2007 for a useful description). Consider two adjacent annuli that, when unwarped, are in force equilibrium both radially and vertically. Then, impose a small tilt on one annulus with respect to the other, such that they are slightly misaligned. As fluid parcels in the two annuli orbit, at two points along the orbit they are maximally separated and at two points along the orbit they are minimally separated. In the frame of the fluid parcel, this manifests as a radial and vertical pressure gradient that oscillates twice an orbit. These pressure gradients induce corresponding oscillations of the particle orbits at the radial and vertical epicyclic frequencies, which, far from the BH, are approximately Keplerian. When the internal stresses of the disk respond to these perturbations, we are left with the oscillating patterns seen in Figure 2. Specifically, the vertical oscillations manifest as m = 2 scale height oscillations while the radial oscillations manifest as an increasing eccentricity above and below the midplane (as also seen in Deng et al. 2021). The argument of periastron above and below the midplane is out of phase by ≈180°.

These oscillations were also seen in the thick, tilted disk simulations of Fragile & Blaes (2008), and were recently computed analytically. Fairbairn & Ogilvie (2021b) developed a theory of oscillating fluid tori that can describe annuli of warped disks, and then in Fairbairn & Ogilvie (2021a) they applied this theory to nonlinear warps in inviscid, Keplerian disks. They identified a bouncing regime above a critical warp amplitude, leading to large-scale height variations (Figure 2 of Fairbairn & Ogilvie 2021a) that are remarkably similar to those seen here. We refer the reader to these works for an analytical analysis of this behavior.

Before continuing, we must describe some of our diagnostics. This includes several averages, generally using tilted coordinates, which we define as

Equation (1)

where x is the coordinate over which we are performing the average, and W is the weight (if we use one). gxx is one of the components of the covariant metric tensor. We sometimes also use this notation for multiple directions, i.e., $\langle \cdots {\rangle }_{\rho }^{\theta ^{\prime} ,\varphi ^{\prime} }$ would indicate a density-weighted average over the tilted coordinates $\theta ^{\prime} $ and $\varphi ^{\prime} $.

We analyze the flow structures depicted in Figure 2 by measuring the pressure scale height,

Equation (2)

where pg is the fluid frame gas pressure. This expression for ${H}_{{p}_{{\rm{g}}}}$ returns the exact scale height, H, for an isothermal thin accretion disk in vertical hydrostatic equilibrium (i.e., for the vertical density and pressure profile $\propto \exp (-{z}^{2}/2{H}^{2})$).

We explore the time-dependent structure of the scale height oscillations in Figure 3, where we show spacetime ($\varphi ^{\prime} -t$) diagrams of ${H}_{{p}_{{\rm{g}}}}/r$ at fixed radii r = 10 rg (panel (b)) and 40 rg (panel (d)). These diagrams are shown from times t ≈ 72–90 × 104 rg/c, during which the inner disk is shrinking and the tearing radius is decreasing. We have depicted our target scale height of H/r = 0.02 with white colors, such that compressed regions are purple and expanded regions are orange. We note, however, that since the "target thickness" of the disk assumes vertical hydrostatic equilibrium, our cooling prescription is a function only of the temperature profile and does not consider vertical oscillations. Additionally, it occurs on a Keplerian timescale, effectively averaging cooling over the annulus. While this is an ad hoc treatment of the disk thermodynamics, we find much of the same behavior in L23 where we self-consistently evolve radiation, reassuring us that our scale height evaluation is robust.

Figure 3.

Figure 3. The phase offset of the m = 2 scale height oscillations is equal to the precession angle (set to $\varphi ^{\prime} =0$) in the outer disk but not in the inner disk. We depict this with spacetime diagrams of ${H}_{{p}_{{\rm{g}}}}/r$ at fixed r = 10 rg (panel (b)) and r = 40 rg (panel (d)). In the time period visualized, the tearing radius is initially 10 rg < rtear < 40 rg, but drifts inward, crossing r = 10 rg at t ≈ 79 × 104 rg/c. This is reflected in the warp amplitude (ψ) plots accompanying each spacetime diagram. In panel (a), we have drawn a red line when the disk tears. We have also drawn a blue dotted ψc = 0.089 line as an estimate for when extreme scale height oscillations begin occurring, calculated using Equation (154) of Fairbairn & Ogilvie (2021a). When t < 79 × 104 rg/c, panel (b) shows the inner disk, which exhibits radially dependent phase offsets (depicted in dark blue) due to strong differential precession. At all times at r = 40 rg, the depicted annulus is part of the outer disk, and the phase of the scale height oscillations is locked with the precession angle. At r = 10 rg, the scale height is also generally larger than the target scale height (=0.02, depicted in white) due to the enhanced dissipation and because vertical oscillations cause a departure from hydrostatic equilibrium.

Standard image High-resolution image

At t = 79 × 104 rg/c, the tearing radius crosses r = 10 rg. Thus, the annulus depicted in the left panel belongs to the inner subdisk before this time and the outer subdisk afterward. Since our tilted coordinate system sets $\varphi ^{\prime} =0$ to the precession angle at every radius, it is physically meaningful to analyze the phase offsets of our scale height oscillations. We derive the phase offsets by fitting ${H}_{{p}_{{\rm{g}}}}(r,\varphi ^{\prime} )$ (Equation (2)) to a sinusoidal dependence that we define at a given radius as,

Equation (3)

While we only use our fit for ${\varphi }_{0}^{\prime} $ in Figure 3, we will return to this expression in Section 4.2. Figure 3(b) highlights the phase offset of the depicted annulus in blue. In the inner subdisk, the offset between the precession angle and the scale height oscillations varies in time, while in the outer subdisk the phase offset is constant. This is also clear in the annulus depicted in Figure 3(d), which at all times belongs to the outer subdisk and has a roughly constant $\varphi ^{\prime} $ dependence. It is interesting that these offsets are constant, because this implies that they evolve with the warp, i.e., $\varphi {{\prime} }_{0}\propto {e}^{i{\omega }_{{ \mathcal P }}}$, where ${\omega }_{{ \mathcal P }}(r)$ is the precession rate, which is nearly uniform in the outer subdisk. Yet, $\varphi {{\prime} }_{0}$ is a function of the warp, and points approximately at the "local" line of nodes between misaligned adjacent annuli (e.g., it is approximately pointed to by the unit vector $\hat{l}\times \partial \hat{l}/\partial r$, where $\hat{l}(r)$ is the angular momentum unit vector of an annulus). This indicates that the entire outer subdisk precesses in a rigid, yet warped, geometry (as expected analytically; e.g., Fairbairn & Ogilvie 2021a).

To support our statements about the tear, in Figures 3(a) and (c) we depict the warp amplitude, $\psi \equiv | \tfrac{\partial \hat{l}}{\partial \mathrm{ln}r}| $, for each spacetime diagram. In the left ψ plot, we have drawn a red line at t ≈ 79 × 104 rg/c to indicate when the disk tears. We have also drawn blue dotted lines at ψc = 0.089, which is an estimate of the "critical" warp amplitude above which the extreme scale height oscillations activate. We obtained this value using Equation (154) of Fairbairn & Ogilvie (2021a) for our target scale height of H/r = 0.02. Figure 3(c) shows that the ψψc criterion is marginally satisfied in the outer subdisk at r = 40 rg. This is a relatively mild warp, suggesting that it may be easy to activate large-scale height oscillations even in disks with initial tilt angles that are much smaller than the 65° angle considered here.

3.3. Tearing Region

When the disk tears, the inner and outer subdisks begin to precess independently. As the two subdisks evolve, they expand until they interact. This interaction is strongest at the line of nodes, where gas parcels orbiting in misaligned planes collide. These gas parcels shock, leading to significant dissipation (discussed later in Section 4.3). This dissipation partially cancels the angular momentum of colliding gas parcels, leading to the formation of low-angular momentum streams of gas that then fall radially inwards. We refer to these as "streamers," and they are an important feature of the inner accretion flow. To better examine the flow in this region, we will define the following diagnostics.

First, we split the mass flow into inward and outward components,

Equation (4)

Equation (5)

where Θ(x) is the Heaviside step function. We also define the average temperature of the flow as,

Equation (6)

where we have used a "$\,\hat{}\,$" symbol as a reminder that the units of $\hat{T}$ are nonphysical. This expression also assumes that the flow is gas-pressure dominated. We expect this is true for the coronal regions, which are generally optically thin, but not for the disk. In Figure 4, we depict ${\dot{M}}_{\mathrm{in}}$, ${\dot{M}}_{\mathrm{out}}$, and $\hat{T}$ in the R'–z' plane (where $R^{\prime} $ is the cylindrical radius in tilted coordinates) at time t = 78854.3 rg/c. We have chosen this time because this is when the inner subdisk is transiently aligned with the BH spin. The inner subdisk then has zero warp and thus no scale height oscillations, allowing us to isolate the effects of the streamers from the tear on the inner subdisk.

Figure 4.

Figure 4. During a tear, the inner and outer subdisks collide, causing streamers of low-angular gas to rain down onto the inner subdisk. Panel (a): we plot contours of inward flowing mass in tilted coordinates at time t = 78854.3rg/c. The mass flux has been integrated in $\varphi ^{\prime} $ and is depicted in the plane of tilted cylindrical ($R^{\prime} $) and vertical ($z^{\prime} $) coordinates. At the tearing radius (rtear ≈ 10 rg), the mass flux is distributed roughly uniformly over a shell. Within the tear, gas plunges radially onto the inner disk, shocking it and increasing the disk mass. Panel (b): same as panel (a), except we depict the outward flowing mass flux. Panel (c): same as panel (a), except we plot the temperature of the gas. The streamer-populated region within the tear is roughly 2 orders of magnitude hotter than the inner subdisk. Panel (d): here, we plot the radial profile of the dimensionless circularization radius of the gas ≡rc/r. Values above (below) unity indicate that the orbital velocity is above (below) its value for a circular orbit. At rtear, rc/r is small, causing streamlines to plunge. Just above the tear, rc/r is above unity. Panel (e): here, we plot the radial profile of the gas pressure in the disk. Since the disk is depleted of gas at the tear, the gas pressure must have a positive radial gradient at rrtear. This causes an inward force that is compensated by the super-Keplerian motion in this region.

Standard image High-resolution image

In Figure 4(a), we can see that ${\dot{M}}_{\mathrm{in}}$ is uniformly distributed about the tearing radius rtear ≈ 10 rg. This is because most of the angular momentum at this radius is dissipated, causing the gas distribution to spread more evenly over a spherical shell, instead of being confined to an annulus. This low-angular momentum gas then forms streamers that rain down onto the inner subdisk, seen as radially extended filamentary structures that sandwich the inner subdisk. Figure 4(b) shows that ${\dot{M}}_{\mathrm{out}}$ follows the structure of ${\dot{M}}_{\mathrm{in}}$ at the tearing radius. This is a consequence of angular momentum conservation; since low-angular momentum streams fall inward, there must also be an outward transport of angular momentum. Since transport by magnetic fields is subdominant (discussed later in Section 4.1), angular momentum must be carried outward by mass. Inside the tear, where there are free-falling streamers, there is essentially zero outwardly flowing gas, as expected. In Figure 4(c), we see that the streamer-populated regions within the tear are roughly 2 orders of magnitude hotter than the outer disk. This is essentially because they evolve on a dynamical timescale, making it difficult for them to cool. They also collide with the inner disk, causing them to shock, heating this region further. The hot, streamer-filled atmosphere is reminiscent of the usual corona that surrounds many accretion disks, suggesting that streamers may lead to hard X-ray emission.

Moving on to the remaining panels of Figure 4, we first define the specific angular momentum of a particle on a circular orbit at r = rc in a Kerr metric as (Shapiro & Teukolsky 1983),

Equation (7)

At every radius, we calculate the BH spin aligned specific angular momentum of the gas, ${{\ell }}_{c}=-\langle {u}_{\varphi }/{u}_{t}{\rangle }_{\rho }^{\theta ^{\prime} ,\varphi ^{\prime} }$, and numerically invert Equation (7) to find the corresponding circularization radius of the gas. We express this dimensionlessly as rc/r and plot it as a function of radius in Figure 4(d). At most radii, rc/r ≈ 1, and the orbital motion is approximately circular. However, at rtear ∼ 10 rg (shown in red), rc drops substantially. This is caused by the cancellation of misaligned angular momenta. The gas at this radius can fall inward until it reaches its local value of rc and form streamers in the process.

Figure 4 also shows that just outside rtear, we have rcirc/r > 1; this means the gas there is super-Keplerian. We can understand why by turning to Figure 4(e), where we depict the gas pressure in the disk, $\langle {p}_{{\rm{g}}}{\rangle }_{\rho }^{\theta ^{\prime} ,\varphi ^{\prime} }$, as a function of radius. The depletion of gas at rtear causes a dip in pressure, which results in a positive pressure gradient at radii ≳rtear. This results in a pressure force that is pointed inward, which is what compensates the super-Keplerian centrifugal force of the gas.

4. Accretion Mechanisms

4.1. Why Do Highly Tilted Accretion Disks Accrete So Rapidly?

In the previous section, we focused on the structure of the warped accretion flow in our simulations. Now, we will study how this structure determines the accretion mechanisms in our disk. In aligned disks, angular momentum transport is usually parameterized by the α parameter, which sets the strength of an effective viscosity and is, in reality, thought to represent magnetized turbulence driven by the magnetorotational instability (MRI; Balbus & Hawley 1991). It is theoretically expected that α < 1. This is because the turbulence should be subsonic and confined by the scale height of the disk; i.e., viscosity is vleddy Veddyα Hcs where leddy < H and Veddy < cs are the characteristic length and velocity scales of the eddies, respectively (Pringle 1981).

Highly tilted and warped accretion disks can, however, accrete at much higher rates (see the reported effective α parameters in Liska et al. 2021, or the dynamically driven accretion reported in Nixon et al. 2012). Figure 5 shows the radial profiles of the effective α parameter, αeff, and the α parameter derived from the Maxwell's stress, αM. Both quantities are averaged over the duration of a tearing event, t ∼ 72 − 90 × 104 rg/c. It is instructive to first look at αeff, which is the effective α parameter that results when assuming mass accretion is fully driven by local viscous processes,

Equation (8)

We can interpret αeff as the α parameter associated with local turbulent stresses if accretion is, in fact, driven by local turbulent stresses. However, as we can see from Figure 5, αeff is ≈1–100 through much of the disk, which is much larger than the α ≤ 1 theoretical limit imposed on turbulent stresses. Additionally, we can compare αeff to αM, which is defined as

Equation (9)

where bμ is the four-vector of the fluid frame magnetic field. Here, αM represents angular momentum transport driven by magnetized turbulence, thought to be seeded by the MRI, and Figure 5 shows that it is bounded below unity as expected. However, since αM is 2–3 orders of magnitude below αeff, we can definitively say that magnetic stresses do not drive accretion in our simulation. 6

Figure 5.

Figure 5. The effective α parameter far exceeds both unity and the magnetic α parameter, suggesting that accretion is neither driven by local turbulence (where α ≲ 1) nor magnetic fields. We show this by plotting radial profiles of the effective α parameter, αeff (Equation (8)), and the α parameter associated with the Maxwell's stress, αM (Equation (9)). These profiles are averaged from times ∼72–90 × 104 rg/c. For both quantities, we also plot their instantaneous profiles at each time, which are depicted in lighter colors. Our time averages are not true steady-state measures of α because disk tearing makes the accretion flow inherently transient. The imprint of disk tearing can be seen in the instantaneous αeff curves, where we have labeled the radius of the tear, which moves inward as the inner disk is accreted.

Standard image High-resolution image

The main takeaways from Figure 5 are as follows: (i) mass accretion occurs much faster in warped accretion disks than in equatorial accretion disks, (ii) the transport is likely nonlocal in nature, since αeff ≫ 1, and (iii) accretion must occur mainly via nonmagnetic stresses, since αMαeff. In the following subsections, we identify the two main accretion mechanisms in our simulation.

4.2. Nozzle Shocks

In Section 3.2 we identified extreme scale height variations that occur twice an orbit. We will now show that these compressions lead to the formation of shocks that dissipate orbital energy. These "nozzle shocks" are conceptually similar to those that occur in tidal disruption events (TDEs; Kochanek 1994; Rees 1988). They were also seen in the thicker tilted disk simulations of Fragile & Blaes (2008), who dubbed them "standing shocks." More recently, Fairbairn & Ogilvie (2021a) also suggested nozzle shocks may form in warped disks.

Figure 6 shows the specific entropy, κg = pg/ργ , in the $\theta ^{\prime} -\varphi ^{\prime} $ plane at radius r = 21 rg and time t = 115, 923.3 rg/c. In a steady, laminar flow, κg should be conserved along streamlines. However, if the gas shocks, then κg will increase, making it a useful quantity for tracking dissipation. We can see that throughout most of the depicted annulus, κg is roughly constant, suggesting that the compressions are mostly adiabatic. We say "mostly," however, because at the points of maximum compression, there is in fact dissipation occurring. To ascertain this, we start by expressing the heating rate per unit volume of a fluid parcel by its the change in entropy 7 (Ressler et al. 2015),

Equation (10)

In practice, we use only positive values of Q, since negative values result in a decrease in entropy along a streamline and are unphysical in a steady flow. 8 While this is a somewhat crude fix, a more precise treatment would require a dedicated heating scheme (i.e., Ressler et al. 2015). To analyze dissipation in our nozzles, we are concerned only with the dissipation of the $\hat{\varphi }^{\prime} $ velocity components, so we will also define ${Q}_{\varphi ^{\prime} }\equiv {\rho }^{\gamma }{(\gamma -1)}^{-1}{u}^{\varphi ^{\prime} }{\partial }_{\varphi ^{\prime} }{\kappa }_{{\rm{g}}}$. We then integrate the azimuthal heating rate vertically,

Equation (11)

Then, we perform a cumulative integral of this quantity in azimuth, ${\displaystyle \int }_{0}^{\varphi ^{\prime} }{F}_{\varphi ^{\prime} }d\varphi ^{\prime} $. We plot this in the bottom panel of Figure 2, where we have normalized our integral to the estimated orbital energy per unit area, ${\rm{\Sigma }}\langle {u}_{t}{\rangle }_{\rho }^{\theta ^{\prime} ,\varphi ^{\prime} }$, where ${\rm{\Sigma }}\equiv {\displaystyle \int }_{0}^{\pi }\sqrt{{g}_{\theta ^{\prime} \theta ^{\prime} }}\rho d\theta ^{\prime} $ is the surface density of the disk. We can then see a very sharp discontinuity almost exactly at $\varphi ^{\prime} \,=\pi $ of scale ≈0.018, indicating that about 1.8% of orbital energy is lost in the nozzle shock. Extrapolating, this would suggest that ∼3.6% of orbital energy is lost every orbit, since there are two nozzle shocks. This is a significant dissipation rate, which we will now show is enough to power the rapid accretion that we reported in Section 4.1.

Figure 6.

Figure 6. The vertical compressions of the warped disk lead to "nozzle shocks" twice an orbit, where significant dissipation occurs. Panel (a): here, we depict the fluid frame entropy (κg), at radius r = 21 rg and time t = 115923.3 rg/c in the $\theta ^{\prime} -\varphi ^{\prime} $ plane. We see that the entropy spikes at $\varphi ^{\prime} \approx 0$ (or 2π) and π, where the disk is most compressed. Panel (b): we plot the cumulative fraction of dissipated energy along the annulus (Equation (11) and following text), normalized to the orbital energy of the annulus. Across the $\varphi ^{\prime} =\pi $ nozzle shock, ≈1.8% of the orbital energy is dissipated.

Standard image High-resolution image

We would like to isolate the dissipation associated with the nozzle shocks, so we will use a criterion to select compressed regions in our integration defined in Equation (11). To do this, we use our sinusoidal fit to ${H}_{{p}_{{\rm{g}}}}(r,\varphi ^{\prime} )$, defined in Equation (3). Then, we only consider the energy dissipation where ${H}_{{p}_{{\rm{g}}}}\lt A(f)$, where $A(f)\equiv {\tilde{H}}_{{p}_{{\rm{g}}}}^{(\mathrm{mean})}-f{\tilde{H}}_{{p}_{{\rm{g}}}}^{(\mathrm{amp})}$ and 0 ≤ f ≤ 1. For the specific choice of $f=\tfrac{1}{\sqrt{2}}$, A is one standard deviation below the mean of Equation (3). We then rewrite Equation (11) with our A(f) criterion,

Equation (12)

Next, we would like to relate this to a predicted mass accretion rate that we can compare with the true mass accretion rate. To do this, we treat FNZ as an axisymmetric, Newtonian dissipation rate associated with a shear viscosity. We can then relate FNZ to a predicted accretion rate (for details, see Pringle 1981) by writing,

Equation (13)

In Figure 7(a), we compare ${\dot{M}}_{\mathrm{NZ}}$ to the simulated $\dot{M}$ for various values of f. Here, we have time-averaged ${\dot{M}}_{\mathrm{NZ}}$ over the period ∼72 − 90 × 104 rg/c (same as Figure 5). Higher values of f indicate a stricter criterion for selecting compressed regions along an annulus. We can see that for f = 0–0.75, ${\dot{M}}_{\mathrm{NZ}}$ matches $\dot{M}$ within an order-unity factor everywhere except the innermost regions. At f = 1, ${\dot{M}}_{\mathrm{NZ}}$ starts departing farther from $\dot{M}$ because our criterion for selecting the dissipation region becomes too strict. To aid our intuition for the "strictness" of our criterion, Figure 7(b) shows the percentage of the disk used in our integration at every radius for each f. We see that when f = 0.75, for which ${\dot{M}}_{\mathrm{NZ}}$ largely accounts for $\dot{M}$ at most radii, we only select roughly ∼20% of an annulus at any given radius, which is an already rather small fraction of the disk. Note that ${\dot{M}}_{\mathrm{NZ}}$ underestimates $\dot{M}$ in the inner regions. This is because the streamers discussed in Section 3.3 can additionally drive accretion, which we explore further in Section 4.3.

Figure 7.

Figure 7. Dissipation in nozzle shocks can mostly account for the measured mass accretion rate everywhere except the innermost regions of the disk. Panel (a): here, we compare $\dot{M}$ to the estimated mass accretion rate due to nozzle shocks, ${\dot{M}}_{\mathrm{NZ}}$ (Equation (13)). We do this for multiple values of f (Equation (12)), where higher f indicates a stricter cutoff for the degree of compression. Each quantity is averaged from times ∼72–90 × 104 rg/c. Panel (b): here, we show what percentage of a given annulus passes our compression criterion when calculating ${\dot{M}}_{\mathrm{NZ}}$, for each value of f. Panel (c): here, we show our fits for the dimensionless mean (${\tilde{H}}_{{p}_{{\rm{g}}}}^{(\mathrm{mean})}/r$) and amplitude (${\tilde{H}}_{{p}_{{\rm{g}}}}^{(\mathrm{amp})}/r$) of the scale height oscillations (Equation (3)). The quantities are time-averaged over the same period as the other panels, with the corresponding instantaneous curves depicted in lighter colors.

Standard image High-resolution image

Figure 7(c) shows our time-averaged estimates for the dimensionless scale height mean, ${\tilde{H}}_{{p}_{{\rm{g}}}}^{(\mathrm{mean})}$, and amplitude, ${\tilde{H}}_{{p}_{{\rm{g}}}}^{(\mathrm{amp})}$. For both, we also plot instantaneous curves in lighter colors to give the reader a better sense of the time-variability of the scale height. At all radii, the mean exceeds the amplitude. In general, the mean tends to be larger than our target scale height H/r = 0.02. The reason for this is twofold: first, the enhanced dissipation rate makes cooling to target thickness more difficult, and second, the cooling prescription does not account for the vertical oscillations, which also inflate the scale height. Both the mean and the amplitude, however, become small at r ≲ 5 rg because the inner region aligns with the BH spin during the phase we have time-averaged, removing the warp.

4.3. Streamers

In Section 3.3, we showed than when our disk tears in two, the subdisks can interact and lead to dissipation. This dissipation results in the formation of low-angular momentum streamers that rain down onto the inner subdisk. In this subsection we will study how this process can enhance accretion in the inner region, as we saw in Figure 4. We follow our approach in Section 3.3 where we focused on a tearing cycle in which the inner subdisk transiently aligns with the BH spin axis. During this aligned phase, the inner disk is unwarped and thus has no nozzle shocks. This allows us to isolate the effects of streamers on the accretion process.

We would like to first qualitatively depict the impact of streamers on the inner subdisk. We do this in Figure 8(a), where we have plotted the entropy density ρ κg in the $\theta ^{\prime} \mbox{--}\varphi ^{\prime} $ plane at radius r = 6 rg. This is aided by Figure 8(c), which shows a 3D rendering of ρ κg, where we have excised gas at radii r > 20 rg to focus on the inner regions. Both snapshots are at time t = 81, 062 rg/c when the tear is at rtear = 7.8 rg. We find entropy density, rather than specific entropy, a useful quantity to plot as it encodes both dissipation in shocks and adiabatic compression. In both 2D and 3D renderings, we see that the streamers collide with the subdisk, resulting in an increase of ρ κg. After they collide, some of the material from the streamers "spills" over the inner subdisk.

Figure 8.

Figure 8. Low-angular momentum streamers produced at the tear crash into the inner subdisk, shocking it and leading to significant dissipation. Panel (a): we depict the fluid frame entropy density (ρ κg) at radius r = 6 rg and time t = 81, 062rg/c. Streamers that originated at rtear = 7.8 rg rain down onto a transiently aligned, unwarped inner subdisk, shocking and dissipating their kinetic energy. Panel (b): we depict the cumulative dissipation, analogous to Figure 6(b), except we show the dissipation in each direction ($r,\theta ^{\prime} ,\varphi ^{\prime} $). We see that we get dissipation of a similar strength to a nozzle shock, except the dissipation is more spread out on the inner disk. The dissipation is still centered at two points along the annulus, but is no longer located at $\varphi ^{\prime} \approx 0$ and π. Panel (c): we depict a 3D rendering of ρ κg at the same time. Here, we have excised gas at radii >20 rg to focus on the inner regions. We can see that streamers from the outer subdisk rain down onto the inner subdisk from either side, and then "spill" over the top.

Standard image High-resolution image

In Figure 8(b), we depict the cumulative dissipation along the azimuthal direction. This is analogous to our calculation in Figure 6 (see also Equation (11) and following text), except that we show dissipation rates for each direction (r, $\theta ^{\prime} $, and $\varphi ^{\prime} $). We do this because the streamer trajectories are significantly altered from the mean flow of the disk. We see that as a streamer crashes into the disk, it shocks, leading to significant dissipation. If we compare this to the nozzle shock dissipation in Figure 2(b), the dissipation is still concentrated at two points along the annulus, because at each radius, streamers will make landfall both above and below the inner subdisk. However, the dissipation is less "peaked" than in a nozzle shock, more spread in azimuth, and not concentrated near $\varphi ^{\prime} =0$ and π.

To make this analysis more quantitative, we will invoke the conservation of mass and BH-aligned angular momentum, and track the fluxes of each entering and leaving the inner subdisk. We start by writing down the mass of the inner subdisk,

Equation (14)

We also track the mass and (BH spin aligned) angular momentum accreted by the inner subdisk from the outer subdisk,

Equation (15)

Equation (16)

We also will use the mass accreted by the BH,

Equation (17)

Finally, we define use the specific angular momentum of the gas accreted onto the inner subdisk as,

Equation (18)

In Figure 9(a), we plot a time series accr normalized to the specific angular momentum of a circular orbit at r = rtear, ≡c (Equation (7)). The time series is plotted over the course of the same tearing cycle depicted in Figures 4 and 8, where the inner subdisk transiently aligns with the BH spin. This is depicted in Figure 9(c), where we plot a spacetime diagram of the tilt angle, ${ \mathcal T }$. Here, the discontinuity between ${ \mathcal T }\approx 0^\circ $ and 65° indicates the tear. At all times during this phase, accr/c ≈ 0.5 − 0.7. This means that all of the gas the inner subdisk accretes will try to circularize at smaller and smaller radii. If the mass of the accreted low-angular momentum gas becomes comparable to the mass of the inner subdisk, then this will cause the inner subdisk to shrink.

Figure 9.

Figure 9. Dissipation at the tearing radius causes the formation of low-angular momentum streams of gas ("streamers") that plunge onto the inner subdisk, causing it to rapidly reduce in size. Panel (a): we depict a time series of the specific angular momentum of gas accreted at the tear (accr, Equation (18)) normalized to the specific angular momentum of a circular orbit at the tearing radius (c, Equation (7)). At all times during the depicted tearing cycle, accr < c, indicating that the gas will try to circularize at smaller and smaller radii. Panel (b): we show a time series of the mass of the inner subdisk (Mdisk), the mass accreted by the BH (∫dMBH, Equation (17)), and the mass accreted by the inner subdisk from the tear (∫dMtear, Equation (15)). Each is normalized to the initial mass of the inner subdisk during the depicted phase (Mdisk,0). At t ≈ 83000 rg/c, there is a transition, when low-angular momentum gas from the tear fully replenishes the initial mass of the inner subdisk. Panel (c): here, we plot a corresponding spacetime diagram of the tilt angle. The tear is delineated by the discontinuity between red and blue regions. By t ≈ 77500 rg/c, the inner subdisk fully aligns with BH spin. The inner subdisk is consumed by the BH at the end of the depicted phase and is quickly replaced by outer disk material that has tilt angle ≈60° − 65°. Panel (d): we show a spacetime diagram of the shell-integrated dissipation rate (∫QdAθ φ ). There is generally dissipation at the tear and at radii r ≲ 5 − 6 rg, where streamers crash into the inner subdisk. Dissipation in these regions is correlated, since dissipation in the tear is associated with an increased formation rate of streamers.

Standard image High-resolution image

In Figure 9(b), we plot three masses normalized to the mass of the inner subdisk at the beginning of the depicted phase, Mdisk,0. We show the mass of the inner subdisk, Mdisk, the total mass the inner subdisk accretes during the phase, $\displaystyle \int {{dM}}_{\mathrm{tear}}=\displaystyle \int {\dot{M}}_{\mathrm{tear}}{dt}$, and the total mass the BH accretes from the inner subdisk, $\displaystyle \int {{dM}}_{\mathrm{BH}}=\displaystyle \int {\dot{M}}_{\mathrm{BH}}{dt}$. Before about t ≈ 83,000 rg/c, the BH accretion rate essentially traces the tear accretion rate, suggesting the inner subdisk acts as a "conveyor belt" of material. Since the accreted gas from the tear has low specific angular momentum, the inner subdisk begins to shrink during this phase. At t ≈ 83,000 rg/c, each curve intersects; the inner subdisk has been depleted of its initial mass by the BH but has also been replenished by the tear. Since all of this replenished material will try to circularize at smaller radii at the time of accretion, the subdisk must shrink. After this transition, the inner subdisk rapidly plunges into the BH, marking the end of the tearing cycle.

Dissipation must be happening for the rapid accretion of the inner subdisk to occur. We know that some of this dissipation occurs at the tear, since streamlines that pass through it lose a large fraction of their angular momentum. However, as we showed in Figure 8, dissipation is also occurring where these streamers merge with the inner subdisk. To get a better sense of the positional dependence of the dissipation rate, we show a spacetime diagram of ∫QdAθ φ in Figure 9(d). We can see that along the tear, there is generally a peak in the dissipation rate, but this peak can vary in strength by about 1–2 orders of magnitude. Correlated with this peak is a spread in the dissipation rate at radii ≲5 rg; these are due to the streamers. In the final stages of the tearing cycle, the streamer dissipation rate is particularly enhanced. Together, this suggests that dissipation at the tear and where the streamers collide with the inner subdisks is comparably important. We also argue that this effect is the cause for the discrepancy between $\dot{M}$ and ${\dot{M}}_{\mathrm{NZ}}$ seen at radii ≲5–10 rg in Figure 7(a), indicating that streamer- and tear-induced dissipation are essential contributors to accretion in torn disks.

5. Discussion

5.1. Expanding the Standard Model of Accretion Disks

The results of Section 4 have profound implications for the understanding of accretion disks. For decades, the standard picture has been that thin, magnetized disks are subject to the MRI that systematically drives angular momentum outward and mass inward, thus enabling accretion (Balbus & Hawley 1991). Although we have demonstrated that magnetic stresses do contribute to accretion (αM curve in Figure 5), they are highly subdominant to the warp-induced dissipation from nozzle shocks (Section 4.2, disk tearing and streamers; Section 4.3). This is not to say the MRI can be neglected in our simulations; were there no magnetized turbulence at all, we expect that the disk would immediately shred because there would be nothing to initially withstand the differential LT torques induced by the BH (as discussed in Section 3.1). Furthermore, the accretion mechanisms studied here are strongest in the inner regions of the disk, as evidenced by the steep radial dependence in Figure 5. Most astrophysical disks are, in reality, much larger in radial extent than the one we have simulated. Active galactic nuclei (AGNs) can harbor Novikov & Thorne (1973) disks up to radii ∼103 − 106 rg (depending on the mass accretion rate) before becoming Toomre-unstable (Toomre 1964; Sirko & Goodman 2003; Thompson et al. 2005). XRB disks are fundamentally limited in size by their Hill radius, which is for instance ∼105 − 106 rg for estimated orbital parameters of Cygnus X-1 (Miller-Jones et al. 2021), although the accretion flow will likely circularize at much smaller radii. So, in misaligned disks hosted by BHs across the mass spectrum, it may be that the MRI drives accretion through most of the disk while the warp drives accretion in the inner regions.

A difficulty of Novikov & Thorne (1973) disks is that their radiation-dominated inner regions are both thermally (Pringle et al. 1973) and viscously (Lightman & Eardley 1974) unstable. This is essentially because any increase (decrease) in temperature results in an increase (decrease) in viscous heating, thus triggering a runaway process. This instability is specific to α-disks, since it is assumed that the viscous stress is ∼α Ptot, which is ∝T4 for a radiation pressure dominated gas. Yet, there is little observational evidence of the thermal-viscous instability of thin disks, which appear to be stable up to significant fractions of the Eddington luminosity (Done et al. 2004). We suggest that nozzle shock-driven accretion is unlikely to be subject to this instability. We expect this because thinner disks are more compressible and thus more susceptible to dissipation in nozzle shocks. If dissipation increases, then the disk will puff up, increasing the scale height and thus decreasing the nozzle shock dissipation rate and maintaining stability. We are unable to probe the stabilizing effect of nozzle shocks in this work since we employ a predefined cooling rate. However, in L23 where we re-run our simulation using an M1 closure scheme for radiation, we find that our disk is in fact thermally stable.

It is important to recognize that we have analyzed a single simulation in this work. Paired with the results of Fragile & Blaes (2008), Nixon et al. (2012), and Liska et al. (2021), it does appear that rapid accretion may be a generic feature of warped disks, but the parameter dependence must be explored before we can ascertain this. One of the most important parameters is the tilt angle. However, if accretion happens at random angles (this is the case if the BH spin and the gas supply have no prior knowledge of one another), then the average tilt angle is 60°. This is only 5° shy of the tilt angle used here, so the accretion mechanisms we have studied may be quite general.

Another critical parameter is the scale height of the disk. Thinner disks will be more strongly influenced by a warp, both for geometric reasons and because their efficient cooling makes them more compressible. So, if we increase the aspect ratio of our disk, the disk will be harder to tear, and warps may be less efficient at driving accretion. It is possible then that the mechanisms discussed in this work are corrections to the accretion process rather than than being the primary drivers of accretion. The scale height is generally set by the Eddington ratio (the ratio of the mass accretion rate of the BH to its Eddington luminosity for a given radiative efficiency), which provides us a more "astrophysical" scale for examining the aspect ratio of the disk. Disks at Eddingtion ratios of ∼0.1%–10% generally cool efficiently and thus lead to thinner disks. At higher Eddington ratios, radiation pressure begins dominating, puffing the disk up. At lower Eddington ratios, cooling becomes inefficient, also puffing the disk up. So, the degree to which warps can affect accretion can also be taken as a function of Eddington ratio.

A particularly interesting laboratory for our accretion mechanisms may be TDE debris disks. This is for a few reasons. First, TDEs have randomly oriented inclinations (and thus an average tilt angle of ∼60°). Second, debris disks are formed from material supplied between the stream self-intersection point and the periastron of the tidal disruption, making them much smaller in radial extent than AGNs or XRBs (Dai et al. 2018; Andalman et al. 2022). Third, although the accretion rate is initially highly super-Eddington, it falls off as ∝t−5/3 after its peak, leading to a sharp drop in the scale height (Shen & Matzner 2014; Tchekhovskoy et al. 2014; Piran et al. 2015). These considerations suggest that warps may be a primary driver of accretion in TDE debris disks, meriting further study.

5.2. Observational Implications

The accretion structure and mechanisms that we have described in the preceding sections will necessarily alter the emission of accreting BHs. One of the most interesting effects of tearing and warp-driven accretion is the resulting variability. While the fractional (say, ∼20%–40%) broadband variability of accretion disks (Gaskell 2004; Uttley et al. 2005; Uttley & McHardy 2005) can occur between viscous and dynamical timescales (by, for instance, a stochastic dynamo action; e.g., Hogg & Reynolds 2016), larger structural changes in an accretion disk are typically expected to occur on a viscous timescale. This is, in turn, limited by the value of α, which is thought to be ≲0.1–1 in standard Novikov & Thorne (1973) disks. However, there is a growing sample of observations of AGNs that exhibit extreme luminosity variations on timescales of months to years, while the viscous timescale can be hundreds to thousands of years. This extreme variability is also not rare; it is exhibited in an estimated ∼30%–50% of quasars (Rumbaugh et al. 2018). These so-called "changing-look" AGNs (CL AGNs; Matt et al. 2003; see also "quasiperiodic eruptions"; Miniutti et al. 2019) are difficult to reconcile with theory, leading some to proclaim a "viscosity crisis" in AGN disks (Lawrence 2018). To explain CL AGNs, it is thought that there must be some instability that leads to the catastrophic and rapid destruction of the inner accretion flow. Currently proposed theories generally invoke a radiation pressure instability (Lightman & Eardley 1974) acting on the inner disk (Janiuk et al. 2002; Sniegowska et al. 2020; Śniegowska et al. 2022a, 2022b), but others have suggested it is the result of a sudden magnetic flux inversion in the disk (Scepi et al. 2021). We argue that the accretion mechanisms presented here may be a natural way of producing CL AGNs. First, as demonstrated in Figure 5, accretion in our simulation happens on timescales that are at least 10–100 times shorter than the usual viscous timescale. The repeated depletion of the inner subdisk will also cause a precipitous drop in the luminosity. This supports the tantalizing hypothesis that some CL AGNs may be the observational result of the tearing process (see also Nixon et al. 2012; Raj et al. 2021). We plan to perform a dedicated comparison of GRMHD disk tearing to CL AGNs in an upcoming work.

Quasiperiodic oscillations (QPOs) form another commonly observed, yet poorly understood, class of accretion disk variability. QPOs are variable signals usually observed in the power spectra of XRBs (van der Klis et al. 1985; Mucciarelli et al. 2006; Gierliński et al. 2008), but have also been observed in TDEs and AGNs as well (Pasham et al. 2019; Smith et al. 2021). The underlying cause of the various kinds of QPOs remains elusive, but possible explanations include the LT precession of tilted disks (Stella & Vietri 1998; Stella et al. 1999; Fragile et al. 2016) or trapped modes excited by warped or eccentric disks (Okazaki et al. 1987; Kato 2004; Ferreira & Ogilvie 2008, 2009; Dewberry et al. 2020a, 2020b). A recent work in our collaboration, Musoke et al. (2022), has performed a separate analysis on this same simulation and has found evidence of both low-frequency (LF) and high-frequency (HF) QPOs. The HFQPOs were associated with radial epicyclic oscillations of the inner subdisk (which were not analyzed in this work), and the LFQPOs were associated with geometric effects due to the precession of the inner subdisk. The variability associated with QPOs is related to, but separate from, any longer-term variability due to the recurrent depletion of the inner subdisk due to streamers from the tear.

Another consideration is the emission produced by streamers. Streamers naturally produce hot, low-density features that surround the inner subdisk (Figure 4). This is reminiscent of the usual accretion disk corona and may result in enhanced hard X-ray emission due to the up-scattering of thermal photons emitted by the inner subdisk. This may help explain the rapid evolution of the X-ray corona observed in some CL AGNs (e.g., Ricci et al. 2020) or contribute to the hard emission observed in XRB state transitions (Esin et al. 1997; Remillard & McClintock 2006), which merits a dedicated study of coronal emission during tearing events.

5.3. Summary

We have performed an analysis of a 3D GRMHD simulation of a highly tilted accretion disk around a rapidly rotating BH, performed at extremely high resolution. We have focused mainly on how the warping and subsequent tearing of the accretion disk impact its geometry and introduce new dissipation mechanisms that drive rapid accretion. Our main findings in this work are as follows:

  • 1.  
    Warped accretion disks drive structural oscillations both vertically and radially. The vertical oscillations manifest as extreme expansions and compressions of the scale height twice an orbit (Section 3.2). The radial oscillations manifest as eccentric streamlines above and below the midplane of the disk. The argument of periapsis for fluid parcels above and below the disk midplane is out of phase by ≈180°. This oscillating solution precesses rigidly with the disk at larger distances (≳10 − 20 rg), but becomes twisted in the rapidly evolving inner subdisk.
  • 2.  
    The oscillating scale height of the disk results in nozzle shocks twice an orbit, dissipating orbital energy and driving rapid accretion. Extreme compressions can shock the gas, leading to enhanced dissipation (Section 4.2). We refer to these as "nozzle shocks," as they are similar to those usually studied in TDEs. These nozzle shocks can lead to rapid accretion, with αeff ∼ 10 − 100, well in excess of that predicted in standard thin disk models, where it is thought that α ≲ 0.1 − 1.
  • 3.  
    Disk tearing results in the formation of low-angular momentum streamers, which rain down on the inner subdisk and drive further accretion. When the disk tears, it can lead to the rapid accretion of the inner subdisk at the tearing radius (Section 4.3). We have attributed this to the cancellation of misaligned angular momentum. This causes low-angular momentum "streamers" to form at the tear and rain down on the inner subdisk. The addition of low-angular momentum gas to the inner subdisk causes it shrink to conserve total angular momentum. This can result in high accretion rates, even in the absence of nozzle shocks or a warped inner disk.

There are several future directions that we would like to consider before concluding. First, the evolution of warped disks is governed by torques acting on both the parallel and perpendicular components of angular momenta, which are mainly determined by the local dissipation mechanisms in the disk. While we have demonstrated that novel dissipation mechanisms (nozzle shocks, tearing, and streamers) play an important role in driving the evolution of our disk, we have done so for a single simulation, and it is unknown how these mechanisms depend on parameters such as the initial tilt and thickness of the disk. Insight could be provided by performing more simulations across the relevant parameter space and by the expansion of existing analytic models to include these dissipation mechanisms. A major frontier to explore is also radiation physics. In this work, we use a predefined cooling function, which simplifies the thermodynamics of the disk. This merits the inclusion of dedicated radiation schemes, which we do in L23. Finally, it is ultimately most important to be able to connect our results to observations, which can be accomplished by producing synthetic observations from simulation results such as those presented here.

Acknowledgments

We thank Y. Lithwick and G. Lodato for useful discussions. An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE), OLCF Director's Discretionary Allocation, and ASCR Leadership Computing Challenge (ALCC) programs under award PHY129. This research used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under contract DE-AC05-00OR22725. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under contract No. DE-AC02-05CH11231 using NERSC award ALCC-ERCAP0022634. We acknowledge PRACE for awarding us access to JUWELS Booster at GCS@JSC, Germany. M.L. is supported by the John Harvard Distinguished Science Fellowship and the ITC fellowship. N.K. is supported by an NSF Graduate Research Fellowship. J.J. and A.T. acknowledge support by the NSF AST-2009884 and NASA 80NSSC21K1746 grants. G.M. is supported by a Netherlands Research School for Astronomy (NOVA), Virtual Institute of Accretion (VIA) postdoctoral fellowship. A.T. was supported by the National Science Foundation grants AST-2206471, AST-2009884, AST-2107839, AST-1815304, OAC-2031997, and AST-1911080.

Appendix: Analyzing Misaligned Accretion Disks

The complicated geometry of warped accretion disks can make them difficult to visualize and analyze. To alleviate this, we perform conversions from "true" (r, θ, φ) to "tilted" ($r,\theta ^{\prime} ,\varphi ^{\prime} $) coordinates. There are two parts to this conversion. First, we have to convert the grid itself. To do this, we calculate radial profiles of the local angular momentum vector and use its orientation to determine the corresponding tilt and precession angles, following Fragile & Anninos (2005). Then, we rotate a θφ grid at each radius such that the precession angle is at $\varphi ^{\prime} =0$ and the tilt angle is =0° (i.e., the disk midplane is at $\theta ^{\prime} =\pi /2$). We then interpolate all relevant quantities to the rotated grid. This is visualized for the fluid frame gas density ρ in Figure 10. In the top row, we can see that an xy slice in true coordinates becomes a face-down image of the disk in tilted coordinates. In the bottom row, we see that a θφ slice in true coordinates shows an annulus rocking up and down due to its tilt angle, while in tilted coordinates the annulus is flattened.

Figure 10.

Figure 10. We depict the conversion of global ("true") coordinates to our radially dependent tilted coordinates. We do so by showing snapshots of ρ in both coordinate systems. The tilted coordinates are indicated by primes and are a function of the local tilt and precession angles (Equation (A1)). At all radii in the tilted coordinate system, $\hat{z}^{\prime} $ is parallel to the local angular momentum vector, $\varphi ^{\prime} =0$ is set to the local precession angle, and $\theta ^{\prime} =\pi /2$ is set to the local midplane of the disk.

Standard image High-resolution image

The second part of this conversion is the coordinate transformation of four-vectors from true to tilted coordinates, which for instance reads ${u}^{\mu ^{\prime} }=\tfrac{\partial {\chi }^{\mu ^{\prime} }}{\partial {\chi }^{\mu }}{u}^{\mu }$ for the four-velocity. The specific matrix elements of the transformation are

Equation (A1)

where ${ \mathcal T }$ and ${ \mathcal P }$ are the tilt and precession angles, respectively. This is a purely Newtonian rotation that leaves the time component of our four-vectors unchanged. In tilted coordinates, the azimuthal velocity ${u}^{\varphi ^{\prime} }$ is aligned with the rotation of its annulus in an average sense. Correspondingly, the vertical velocity ${u}^{\theta ^{\prime} }$ averaged over an annulus is approximately zero.

Footnotes

  • 6  

    We note the absence of an α parameter associated with Reynolds' stresses in Figure 5; this was an intentional choice. The oscillating flow structures in our warped disk break the ergodicity of the disk in the $\hat{\varphi }^{\prime} $-direction and would merit a spectral Reynolds' decomposition of the flow into its Fourier modes, which is beyond the scope of this work.

  • 7  

    Entropy is also lost via our cooling function. However, we do not need to consider it in this expression, as we are only concerned with the entropy dissipated in shocks—not the entropy that's removed by the cooling function.

  • 8  

    While our flow is generally transient, this is a fine assumption in the $\hat{\varphi ^{\prime} }$-direction since the orbital timescale is generally much shorter than any other timescale.

Please wait… references are loading.
10.3847/1538-4357/ace051