Tilted Disks around Black Holes: A Numerical Parameter Survey for Spin and Inclination Angle

TL;DR

This study explores how tilted accretion disks around spinning black holes behave, revealing complex precession, shock formation, and impacts on flow efficiency through extensive numerical simulations.

Contribution

It provides a comprehensive numerical survey of tilted disk properties across various black hole spins and tilt angles, highlighting new flow dynamics and shock effects.

Findings

Tilted disks precess rigidly around black hole spin axes.

Standing shocks form at high inclinations, affecting angular momentum transport.

High spins lead to more effective material redirection and reduced radiative efficiency.

Abstract

We conduct a systematic study of the properties of tilted accretion flows around spinning black holes, covering a range of tilt angles and black hole spins, using the general-relativistic magnetohydrodynamics code Athena++. The same initial magnetized torus is evolved around black holes with spins ranging from 0 to 0.9, with inclinations ranging from 0 degrees to 24 degrees. The tilted disks quickly reach a warped and twisted shape that rigidly precesses about the black hole spin axis with deformations in shape large enough to hinder the application of linear bending wave theory. Magnetized polar outflows form, oriented along the disk rotation axes. At sufficiently high inclinations a pair of standing shocks develops in the disks. These shocks dramatically affect the flow at small radii, driving angular momentum transport. At high spins they redirect material more effectively than they heat it, reducing the dissipation rate relative to the mass accretion rate and lowering the radiative efficiency of the flow.

Equations

$$Q_{\phi}\equiv\frac{2\pi(r^{3/2}+a)}{\sqrt{g_{33}}\Delta\phi}\bigg{\lparen}\frac{b_{3}b^{3}}{\rho+\Gamma/(\Gamma-1)\times p_{\mathrm{gas}}+2p_{\mathrm{mag}}}\bigg{\rparen}^{1/2}$$

$$h/r\equiv \frac{\int |{\theta }^{\prime }-\pi /2|\rho \sqrt{-g}\mathrm{d}\theta \mathrm{d}\varphi }{\int \rho \sqrt{-g}\mathrm{d}\theta \mathrm{d}\varphi },$$

$$⟨{\beta }^{-1}⟩\equiv \frac{\int {\beta }^{-1}\rho \sqrt{-g}\mathrm{d}r\mathrm{d}\theta \mathrm{d}\varphi }{\int \rho \sqrt{-g}\mathrm{d}r\mathrm{d}\theta \mathrm{d}\varphi }$$

$$\text{tensor}\ast L{}^{\hat{ı}}\equiv [ijk]\text{tensor}\ast r{}^{\hat{ȷ}}\text{tensor}\ast T{}_{\mathrm{gas}}^{0\hat{k}},$$

$$\tensor*{T{\mathstrut}}{{}_{\mathrm{gas}}^{\mu\nu}}\equiv\bigg{\lparen}\rho+\frac{\Gamma}{\Gamma-1}p_{\mathrm{gas}}\bigg{\rparen}\tensor*{u{\mathstrut}}{{}^{\mu}}\tensor*{u{\mathstrut}}{{}^{\nu}}+p_{\mathrm{gas}}\tensor*{g{\mathstrut}}{{}^{\mu\nu}}$$

$$i$$

$${\varphi }_0$$

$${\theta }^{\prime }$$

$${\varphi }^{\prime }$$

$$-\sin i\cos \theta +\cos i\cos {\varphi }_0\sin \theta \cos \varphi +\cos i\sin {\varphi }_0\sin \theta \sin \varphi ),$$

$$\theta$$

$$\varphi$$

$$\sin i\cos {\varphi }_0\cos {\theta }^{\prime }+\cos i\cos {\varphi }_0\sin {\theta }^{\prime }\cos {\varphi }^{\prime }-\sin {\varphi }_0\sin {\theta }^{\prime }\sin {\varphi }^{\prime }).$$

$$w$$

$$s$$

$$\dot{M}$$

$$\Phi$$

$$\phi$$

$${\nabla }_{\mu }(s\rho u^{\mu })=s{\nabla }_{\mu }(\rho u^{\mu })+\rho u^{\mu }{\nabla }_{\mu }s=\rho u^{\mu }{\partial }_{\mu }s.$$

$$({\nabla }_{\mu }\text{tensor}T{}_{\nu }^{\mu }{)}_{\mathrm{adi}}=-q{u}_{\nu }$$

$$Q\equiv \int -u_0{p}_{\mathrm{gas}}{u}^i({\partial }_{i}s)\sqrt{-g}\mathrm{d}r\mathrm{d}\theta \mathrm{d}\varphi ,$$

$$\text{tensor}\ast T{}_{\mathrm{Rey}}^{\mu \nu }$$

$$\text{tensor}\ast T{}_{\mathrm{Max}}^{\mu \nu }$$

$$\alpha \equiv \frac{\int |T^{x^{\prime \prime }{y}^{\prime \prime }}|/(p_{\mathrm{gas}}+p_{\mathrm{mag}})\times \rho \sqrt{-g}\mathrm{d}\theta \mathrm{d}\varphi }{\int \rho \sqrt{-g}\mathrm{d}\theta \mathrm{d}\varphi }.$$

$$\Sigma u^0$$

$$\Sigma u^1$$

$$\Sigma u^2$$

$$\Sigma u^3$$

$$\Delta$$

$$\Sigma$$

$$P$$

$$R$$

$$\Theta$$

$$\text{tensor}\ast \Gamma {}_{\mu \nu }^1\text{tensor}\ast u{}_{(0)}^{\mu }\text{tensor}\ast u{}_{(0)}^{\nu }=0,$$

$$\text{tensor}\ast g{}_{\mu \nu }\text{tensor}\ast u{}_{(0)}^{\mu }\text{tensor}\ast u{}_{(0)}^{\nu }=-1,$$

$$\displaystyle\bigg{\lparen}\frac{1}{\Delta}\big{\lparen}r^{2}+a^{2}\big{\rparen}^{2}-a^{2}\sin^{2}\!\theta\bigg{\rparen}e-\frac{2Mar}{\Delta}\ell$$

$$\displaystyle\frac{2Mar}{\Delta}e+\bigg{\lparen}\frac{1}{\sin^{2}\!\theta}-\frac{a^{2}}{\Delta}\bigg{\rparen}\ell$$

$$\frac{\mathrm{d}\phi}{\mathrm{d}\theta}\equiv\frac{u^{3}}{u^{2}}=\pm\bigg{\lparen}\frac{\ell}{\sin^{2}\!\theta}-ae+\frac{ae(r^{2}+a^{2})-\ell a^{2}}{r^{2}-2Mr+a^{2}}\bigg{\rparen}\bigg{\lparen}q-\bigg{\lparen}a^{2}\big{\lparen}1-e^{2}\big{\rparen}+\frac{\ell^{2}}{\sin^{2}\!\theta}\bigg{\rparen}\cos^{2}\!\theta\bigg{\rparen}^{-1/2}.$$

$$\Delta \varphi =4{\int }_{\pi /2+i}^{\pi /2}\frac{\mathrm{d}\varphi }{\mathrm{d}\theta }\mathrm{d}\theta -2\pi .$$

Full text

accretion disks, black hole physics, hydrodynamics, MHD, relativity

††software: Athena++ (White et al., 2016)

1 Introduction

The general-relativistic (GR) effects of nonzero spin can have important consequences for accretion flows onto Kerr black holes. These effects are most dramatic when the spin of the black hole is not aligned with the rotational axis of the accretion disk, breaking the axisymmetry of the system. Such tilted disks are expected to occur in nature, for instance in supermassive black holes where the material falling from large radii in general has no reason to be aligned with the black hole spin. Lense–Thirring precession of a tilted flow is also a popular model for type C low-frequency quasi-periodic oscillations in black hole X-ray binaries (Ingram et al., 2009).

First simulated in GR by Fragile & Anninos (2005) and Fragile et al. (2007) with the Cosmos and Cosmos++ codes, tilted accretion disks display dynamical behavior not seen in untilted systems. For example, the monotonic relation between spin and the inner edge of the disk breaks down (Fragile, 2009), complicating black hole spin measurements. Even more noteworthy is the shape of tilted disks, as radial variation in the strength of Lense–Thirring precession results in disks warping and twisting as they precess around the black hole spin axis.

Simulations have shown that tilted disks can develop a pair of standing shocks that can alter the transport of angular momentum and dissipation of energy in the disk (Fragile & Blaes, 2008; Generozov et al., 2014). Such structures may directly affect the variability observed in accreting systems (Henisey et al., 2012). Tilt also comes into play when relativistic jets are considered. For example, Polko & McKinney (2017) have studied the effect of relativistic jets on orienting tilted disks, and Liska et al. (2018) have examined how jets in such systems are oriented and how they precess.

The evolution of the shape of a tilted disk is often categorized as either as a diffusive process (Bardeen & Petterson, 1975; Papaloizou & Pringle, 1983) or in terms of bending waves (Papaloizou & Lin, 1995), depending on whether the scale height $h/r$ is small or large compared to the effective $\alpha$-viscosity parameter. In the diffusive regime, the inner disk may align with the black hole spin, with the outer disk remaining misaligned. This will generally not happen in the bending wave regime, which is where most numerical work is done. Even when simulations can resolve geometrically thin accretion flows ($h/r\lesssim 0.1$), it can be prohibitively expensive to resolve flows so thin as to have $h/r<\alpha$.

Krolik & Hawley (2015) contest that the distinction between diffusion and bending waves is based on isotropic viscosity models, and that in magnetohydrodynamics (MHD) the only important distinction is the degree of nonlinearity in the bending waves. As shown by Sorathia et al. (2013) and Krolik et al. (2014) using finite-volume Zeus simulations, MHD effects cannot be neglected in studying these disks. Even though magnetic forces on tilted disks are generally small compared to hydrodynamic ones, the anisotropy of MHD stresses and the turbulent nature of flows subject to the magnetorotational instability (MRI, Balbus & Hawley, 1991) have important consequences for the transport of angular momentum, including the component of angular momentum responsible for alignment.

Even where numerical work has been done in the traditional bending-wave regime, there is a lack of consensus on what the final shape of the disk is. In contrast to the aforementioned authors, Nelson & Papaloizou (2000), and later Nealon et al. (2015) and Nealon et al. (2016) using the smoothed-particle hydrodynamics Phantom code, find that disks can break, with fragments at different radii inclined relative to one another. This may be a result of the disk being unable to radially transport differentially applied Lense–Thirring torques rapidly enough. While direct comparison between such different numerical methods using the same physical conditions is difficult, the existence of these discrepancies between simulations indicates that the codes may not be in agreement regarding the rate of angular momentum transport relative to GR precession.

Neither Zeus nor Phantom employ GR directly, but rather add post-Newtonian corrections to capture the leading-order nodal (Lense–Thirring) and possibly apsidal (Einstein) precession terms, the magnitudes of which are sometimes artificially increased in order to observe the effects of differential precession in fewer orbital periods. Phantom is limited to simulating hydrodynamics with an artificial, isotropic viscosity, while the Zeus results employ MHD to self-consistently induce turbulence and drive angular momentum transport. As it is the interplay between GR dynamics and angular momentum transport that shapes the disk, these concerns need to be resolved. This motivates the use of a code that employs both GR and MHD, naturally including the post-Newtonian terms at all orders and directly modeling the MRI turbulence we know to be essential in accretion processes.

Other studies on tilted disk dynamics have gone so far as to employ full numerical relativity and thus capture the gravitational effect of the disk itself (Mewes et al., 2016a, b), though these are still done in viscous hydrodynamics rather than MHD. While the Cosmos++ code uses both MHD and stationary GR, the aforementioned results were only able to marginally resolve the MRI. Moreover, those studies, like most of those done to date, examine one or two tilted configurations and compare them to an untilted control disk. As the shape and state of disks are apparently sensitive enough to GR effects and angular momentum transport that oversimplified approximations can lead to notably different outcomes, it is natural to ask how outcomes vary with physical parameters, assuming numerics are handled adequately.

Here we more thoroughly explore a slice of parameter space, varying black hole spin and disk tilt angle, focusing on how these parameters affect the dynamics of tilted accretion disks. We use the finite-volume GR ideal MHD capabilities of Athena++ (White et al., 2016) to evolve ten similar initial conditions with different combinations of spin and initial tilt. Our goal is to systematically vary these parameters keeping all else fixed, so that we may quantify how the accretion flows depend on these important variables. At the same time, we want to be sure that we are including all relevant physics (GR and MHD) in order to accurately simulate the behavior of tilted disks.

The simulations are described in §2, with important definitions regarding tilted coordinate systems given in §3. The results are analyzed in terms of the shape of the disk midplane (§4), the resulting poloidal structure (§5), and the standing shocks we find (§6). Concluding remarks are made in §7. As these are non-radiative, non-self-gravitating, ideal MHD simulations, the results scale to any black hole mass $M$, with all lengths implicitly in units of $GM/c^{2}$ and times in units of $GM/c^{3}$.

2 Numerical Setup

We choose to use spherical Kerr–Schild coordinates (Font et al., 1998) aligned with the black hole axis. Our coordinates cover the entire sphere, including the appropriate transmissive polar boundary, and extend from inside the horizon to $r=70$. In order to keep the timesteps reasonable, we use static mesh refinement to keep the part of the grid near the poles at low resolution. The regions within $22.5^{\circ}$ of the poles are at an effective resolution of $56\times 32\times 64$ in $r$, $\theta$, $\phi$; regions between $22.5^{\circ}$ and $33.75^{\circ}$ from the poles are at $112\times 64\times 128$; and the grid within $56.25^{\circ}$ of the midplane is at $224\times 128\times 256$. For comparison, the fiducial resolution of Fragile et al. (2007) was $128^{3}$ (with the outer radius of the simulation at $r=120$). Improvement of the azimuthal resolution in particular means that at late times we achieve a relativistic MRI quality factor

[TABLE]

in the inner parts of the disk ranging from approximately $70$ (no spin) to $90$ (highest spin). Hawley et al. (2011) suggest convergence sets in around values of $15$ to $20$. We note that the results of §6 indicate that while the MRI may be necessary to have any accretion, once the flow develops standing shocks it is the adequate resolution of these hydrodynamic features that may become more important.

In choosing to incline the disk relative to the grid, rather than have the two be initially aligned, we are limited in how extreme a tilt we can achieve before the disk lies in the lower resolution polar region. However, this choice means the disk is free to precess about the spin axis without approaching the grid axis. For example, if we had the disk and grid be initially aligned, and the spin axis misaligned from these by $24^{\circ}$, then after half a precession the disk would still be $24^{\circ}$ off the spin axis but $180^{\circ}$ out of phase with the grid, with the result being the disk and grid misaligned by $48^{\circ}$ and much of the disk worryingly close to the grid polar axis. Even introducing a global precession into the grid itself will not solve this, given that the disk can precess differently at different radii.

Athena++ is run with a second-order van Leer integrator at a Courant–Friedrichs–Lewy number of $0.3$, where $1/3$ is the maximum for this integrator in $3$ spatial dimensions. We use the second-order modified van Leer spatial reconstruction from Mignone (2014). The HLLE Riemann solver is used to calculate fluxes.

The initial conditions are those of the hydrostatic equilibrium solution of Fishbone & Moncrief (1976) with inner edge at $r=15$, pressure maximum at $r=25$, adiabatic index of $\Gamma=4/3$, and peak rest mass density of $\rho=1$. These are the same inner edge and pressure maximum values as in Fragile et al. (2007), though we have a more compressible adiabatic index compared to their $\Gamma=5/3$. These are thick tori with no cooling added, so at late times the scale height

[TABLE]

ranges from $0.1$ to $0.2$ inside $r=20$. This solution is calculated for an untilted disk and is then given the desired inclination (so it is no longer exactly in equilibrium).

We next add a magnetic field from a vector potential that has only an azimuthal component in disk-aligned coordinates proportional to $\max(\rho-0.2,0)$. The strength of the magnetic field is normalized with respect to plasma $\beta^{-1}\equiv p_{\mathrm{mag}}/p_{\mathrm{gas}}$ such that

[TABLE]

is $10^{-3}$.

We perform ten simulations out to a time of $t=4000$, with snapshots made every $10$ time units. Two of these are around a Schwarzschild black hole, one with an initial disk inclination of $0^{\circ}$, and another with an inclination of $24^{\circ}$ used to verify that the results do not depend on the orientation of the disk with respect to the grid. We also run simulations with $0^{\circ}$, $8^{\circ}$, $16^{\circ}$, and $24^{\circ}$ inclinations around black holes with dimensionless spins $a=0.5$ and $a=0.9$. Steady-state inflow is established out to approximately $r=10$ by the end of the simulations.

3 Measuring Disk Properties

In order to analyze the properties of a tilted accretion disk, we must have a definition of the disk’s orientation at a moment in time, either globally or as a function of radial coordinate. We choose the gas angular momentum defined to be

[TABLE]

where hats indicate Cartesian coordinates related to our underlying spherical coordinates in the usual way, the brackets denote the antisymmetric Levi–Civita symbol, and

[TABLE]

are the components of the hydrodynamic stress-energy tensor. The results are much the same if instead of $T_{\mathrm{gas}}^{0\mu}$ we choose $\rho u^{0}u^{\mu}$ or even $\rho u^{\mu}/u^{0}$ as is done by some authors.

From the angular momentum we can construct the Euler angles $i$ (inclination, sometimes called $\beta$ in the literature) and $\phi_{0}$ ($90^{\circ}$ less than the longitude of ascending node, sometimes called $\alpha$). They are given by

[TABLE]

These angles relate the coordinates $\theta,\phi$ aligned with the grid and black hole to the coordinates $\theta^{\prime},\phi^{\prime}$ aligned with the disk according to

[TABLE]

Note that for a circular geodesic, $i$ is not exactly the maximum deviation of $\theta$ from the midplane over the orbit. However, these two senses of inclination are in close agreement, never deviating by more than approximately a degree over all radii in the range of spins and tilts considered.

4 Disk Shape

In each of our simulations, the disk orientation as a function of radius stops evolving by $t=1500$, except for a slow, rigid-body precession. This precession rate is $2.4^{\circ}$ per $1000$ time units in the $a=0.5$ cases, and $3.8\text{--}4.1^{\circ}$ per $1000$ time units in the $a=0.9$ cases. For comparison, Fragile et al. (2007), whose initial conditions match ours in terms of the inner edge and pressure maximum radii, find a rigid-body precession period of approximately $0.3\ (M/M_{\odot})\ \mathrm{s}$ for $a=0.9$ with a $15^{\circ}$ inclination. This corresponds to $6^{\circ}$ per $1000$ time units. As with Fragile et al., our precession rate is consistent with the total angular momentum of the rigid body and the Lense–Thirring torque applied to it. As a measure of conservation of angular momentum when the disk and grid are not aligned, the $24^{\circ}$ inclined disk around the Schwarzschild black hole precesses less than $0.1^{\circ}$ per $1000$ time units.

The radially dependent orientation defined by the Euler angles $i$ and $\phi_{0}$, averaged over $3000\leq t\leq 4000$, is shown in Figures 1 and 2. We mark the equatorial innermost stable circular orbit (ISCO) for reference. As the disk is rotating rigidly at a rate much slower than the precession of geodesics at small radii, we can time-average $\phi_{0}$ without worrying that we are averaging over many full precession periods at these radii. In fact, over the period of the time averaging, the solid body rotates less than $10^{\circ}$ at all radii.

The small deviation in $i$ of the initially untilted disk with $a=0.9$ (Figure 1, right panel, black line) measures the imprecision of this definition due to the stochastic nature of the accretion flow, as well as magnetic field angular momentum not included in (4). The relative flatness of the initial $24^{\circ}$ inclination when there is no spin (left panel, colored line) demonstrates that the code does well at conserving angular momentum even when the flow is not aligned with the grid or coordinate system.

The most striking feature of Figure 1 is that the inclinations increase with decreasing radius. That is, the disks are warped even further out of the plane near the black hole. Only at the highest spin and lowest inclination is this not seen. This additional warping is observed in other simulations of tilted thick disks; see for example Figure 12 of Fragile et al. (2007) or Figure 7 of Mewes et al. (2016b).

One way this effect can arise is if angular momentum exchange between adjacent radii occurs primarily at the ascending and descending nodes. (See §6 for a description of the standing shocks found in our simulation, whose locations near the nodes can cause this to happen.) The inclination angle defined in (6a) is notably different from the inclination of the orbit’s direction while crossing the line of nodes, $\pm\tan^{-1}(u^{2}/u^{3})$. In particular, for a constant inclination as defined by angular momentum (or by antinodal latitude), this nodal inclination decreases with decreasing $r$. Conversely, if nodal inclination is kept constant with radius, the angular momentum inclination of circular geodesics must increase toward the black hole. Thus if the stress responsible for transporting angular momentum tends to align the trajectories of adjacent fluid elements, and if this mechanism operates primarily near the line of nodes, one expects to see inclination increasing toward smaller radii in steady state.

With Figure 2 we again have a measure of the numerical non-conservation of angular momentum, given by the small deviation of the $24^{\circ}$ inclined disk around the Schwarzschild black hole (left panel, colored line). For the inclined flows around spinning black holes, we see the disks being twisted by differential precession. This effect is strong enough in the $a=0.9$ case that the very inner parts of the disk can be twisted beyond $180^{\circ}$ relative to the outer parts.

While we observe a large amount of differential precession in these simulations, we note that the effect would be much stronger if not for internal stresses redistributing angular momentum. To see this, we consider a model in which particles orbit on circular geodesics, accumulating a precession angle at a rate $\mathrm{d}\phi_{0}/\mathrm{d}\tau$ with respect to their proper time (see the appendix, §A). This can be converted to a differential precession $\mathrm{d}\phi_{0}/\mathrm{d}r$ by dividing by a prescribed radial velocity $u^{1}\equiv\mathrm{d}r/\mathrm{d}\tau$.

We show the results of this calculation in Figure 3. The solid lines are the observed precession angles as shown in Figure 2. The dashed lines are computed by integrating $\mathrm{d}\phi_{0}/\mathrm{d}r$ inward from $r=20$, where the curves are matched to the observed values. The values of $u^{1}$ used are taken from the respective simulations; they are density-weighted shell averages within one scale height of the disk midplane and averaged in time over $3000\leq t\leq 4000$. The fact that the dashed lines are higher than the solid lines indicate the disk is transporting angular momentum induced by the Lense–Thirring effect from small radii to large radii.

The spread in the observed $\phi_{0}$ profiles matches that of this geodesic calculation, and this is due almost entirely to differences in $u^{1}$ between simulations at different inclinations. That is, the differences in differential precession between simulations with the same value of $a$ are dominated by differences in radial velocity, not by minute differences in precession of geodesics at different inclinations. While nodal precession of circular geodesics does have a small dependence on inclination (beyond the linear order used by Lense and Thirring), this effect is hardly noticeable at the scale of the plot. Instead of using the radial velocities from each simulation we can use the values from the untilted simulation at the appropriate spin, thus removing the effect of different $u^{1}$ values, and the result is shown as the black dotted line in Figure 3.

As there are internal aligning and precessing torques in these thick disks, one might consider comparing our results to the theory of bending waves. Convenient formulas for the solution to the linearized bending wave problem are given in Foucart & Lai (2014). This framework is developed in a Newtonian context with a nonspherical potential. The solution depends on the orbital, epicyclic, and vertical frequencies of orbits, as well as the integrated density and pressure of the disk and an effective viscosity for the disk.

The result of this comparison, however, only highlights the inapplicability of linear bending wave theory to most realistic tilted disks around black holes. The theory makes a number of assumptions that are violated here: that the epicyclic and vertical frequencies deviate only slightly from their Keplerian values (see the appendix, §B, for the definitions and plots of these in Kerr spacetime), that there is no stress at the ISCO, and that $h/r\gg\alpha$. In order to illustrate this incompatibility, we show the observed and theoretical disk inclinations in Figure 4.

We can quantify the out-of-plane bending of the disk, defining the dimensionless warp and twist to be

[TABLE]

Here $\vec{n}$ is the Cartesian vector of angular momentum components $L^{\hat{\imath}}$, scaled to have unit norm. These are similar to the definitions given in Shiokawa (2013), but we use the radial coordinate rather than scale height to make them dimensionless. The nontrivial warps and twists from our simulations are plotted in Figures 5 and 6.

5 Poloidal Structure

Here we examine the structure of our tilted accretion flows in the poloidal plane. For the plots in this section we slice the simulation along the plane orthogonal to the line of nodes, orienting all images such that the black hole spin axis is vertical. Thus the disks will always be tilted down to the right. These slices are averaged in time for $3000\leq t\leq 4000$. Streamlines are constructed from the $r$- and $\theta$-components of vector quantities. In all cases we show the inner $40\times 40$ gravitational radius patch of the simulation, with the black hole masked in the center.

Slices of rest-frame density $\rho$ are shown in Figure 7. From these it is clear that most of the disks display large bending out of their original planes. For the highest inclinations, the disk even appears to split into two streams: a thinner one close to vertical in the image, and a broader one. The fact that these streams remain clearly visible even after the long time averaging indicates they are not transient phenomena.

The density plots indicate that an evacuated, low-density region forms along the direction of the disk’s angular momentum rather than parallel to the black hole spin. Untilted simulations generically see outflows along the common axis, so this raises the question of what direction the material is flowing in the tilted case. For this reason we construct Figure 8 showing the velocity streamlines of $u^{1}$ and $u^{2}$, with colors indicating the radial velocity.

In the untilted cases, we clearly see the polar outflow strengthening with increasing spin. Comparing tilted disks to untilted disks at a given nonzero spin, we see that this outflow is overall suppressed in the former. Still, its direction is generally aligned with the low-density regions.

By comparing the lower right panels in Figures 7 and 8, we can analyze the velocity structure in relation to the disk location. In particular, there is a pattern of positive $u^{1}$ above the disk on the right-hand side and below the disk on the left, with negative $u^{1}$ in the opposite positions. The outward radial velocity is not an unbound outflow, but rather the part of an eccentric orbit going from periapsis to apoapsis. Above the disk, the matter is orbiting with apoapsis behind the figure and periapsis in front, with the reverse holding below the disk. This pattern persists through time averaging, showing that it is fixed in place relative to the rigidly precessing disk. These eccentric orbits, with those above the disk $180^{\circ}$ out of phase relative to those below and with apsides near to the disk’s line of nodes, are observed in Fragile & Blaes (2008, cf. Figure 1), where they are found to be signatures of standing shocks. The coherence of the velocity pattern decreases with decreasing spin and tilt, so that we expect shocks to be strongly influencing the accretion flows only in our simulations with higher tilts and spins.

Also worth viewing in poloidal slices is the magnetization of the material. For this we choose $\sigma\equiv 2p_{\mathrm{mag}}/\rho$, shown in Figure 9. Here streamlines are constructed from $B^{1}$ and $B^{2}$. In a time-averaged sense, there is a current sheet in a relatively narrow region roughly aligned with the midplane of the disk. A radial current is expected whenever a disk with a poloidal magnetic field rotates faster at its midplane, shearing the field azimuthally with a reversal at the midplane. Tilting the disk does not change this general feature except by warping the current sheet.

Unsurprisingly, the low-density polar regions seen in Figure 7 correspond to the higher magnetizations in Figure 9. It is more noteworthy that the average magnetic field orientation in this region appears to have less inclination to the black hole spin axis (vertical in the figures) than the disk’s rotational axis does. This is despite the initial field being aligned with the disk, and in contrast with the average velocity structure in the same regions, which is close to perpendicular to the plane of the disk. Thus the velocity and magnetic structures of a jet may carry information about the relative inclination of the disk to the black hole.

Note these are only weak polar outflows, as expected given that we did not use initial conditions designed to yield a magnetically arrested disk (MAD). Define

[TABLE]

where the time average in the last equation is taken over $3000\leq t\leq 4000$. We find $\varphi\lesssim 8$ in all our simulations, whereas the prototypical strong-jet-producing MAD simulations of Tchekhovskoy et al. (2011) have $\varphi\approx 47$.

6 Standing Shocks

In order to investigate standing shocks, we first construct averages of density, gas pressure, and velocity over the time interval $3000\leq t\leq 4000$, rotating each snapshot in $\phi$ in order to compensate for the global precession rate. The quantities other than density are density weighted in this average. At each point in this average state, assumed to be stationary, we compute the nonadiabatic source term for generation of entropy per unit coordinate volume per unit coordinate time, $\rho u^{i}\partial_{i}s$, where $s\equiv\log(p_{\mathrm{gas}}/\rho^{\Gamma})/(\Gamma-1)$. This follows from writing a “conservation law” for entropy per unit volume analogous to that for mass per unit volume:

[TABLE]

The second equality follows from mass conservation $\nabla_{\mu}(\rho u^{\mu})=0$. The right-hand side must be the nonadiabatic entropy source per unit time and volume in the rest frame of the fluid, whose direct evaluation is normally hindered by not knowing $\partial_{0}s$ given only $s$ at a single time slice. However in steady state this term must vanish, and we are left with the source term $\rho u^{i}\partial_{i}s$. Dividing by $u^{0}$ converts the rate from that measured in the fluid frame to that measured with respect to Kerr–Schild time $t$, but we must also multiply by $u^{0}$ to obtain a rate per unit coordinate volume rather than per unit fluid-frame volume. Multiplying by temperature $T\equiv p_{\mathrm{gas}}/\rho$ gives the irreversible heating rate per unit coordinate time per unit coordinate volume $q\equiv p_{\mathrm{gas}}u^{i}\partial_{i}s$. This term should have large positive values in regions of shocks.

Plots of the spatial distribution of this source term are shown in Figure 10. In order to make these plots, we orient each shell of material according to its own angular momentum. The longitude of ascending node at $\phi^{\prime}=\pi/2$ is the upward vertical axis in the figures, with the angular momentum vector coming out of the page. As a result the plane of the figure is warped and twisted in $r,\theta,\phi$ coordinates. The source term is volume-averaged over a scale height in $\theta^{\prime}$ at each point $r,\phi^{\prime}$. The solid black circles indicate the event horizon, the gray shaded regions are inside the prograde photon sphere radius $r_{\mathrm{ph}}$, and the dashed black lines denote the ISCO, all calculated in the black hole’s midplane.

From the figure we can see a trend of stronger, more well-defined shocks at higher inclinations. Much of this happens outside the photon sphere, which can be used as a rough estimate for the boundary inside of which extra dissipation will have only negligible observable consequences (Generozov et al., 2014). In fact, the shock structure extends well beyond the ISCO.

More quantitatively, we can integrate the heat source term over the domain to gauge the extra heating caused by the shock. To do this, we must first incorporate our heating term into a proper tensor equation for energy evolution. In the fluid rest frame, we know $q$ must be the source for contravariant energy density $T^{00}$. As this source term should be the time component of a four-vector, it is more appropriate to say $qu^{0}$ is the source for $T^{00}$, as $u^{0}=1$ in this frame. More generally, the adiabatic evolution of conserved (covariant) energy-momentum density must then obey

[TABLE]

in all frames. Physically, the factor of $-u_{0}$ multiplying the right-hand side when considering energy ($\nu=0$) accounts for gravitationally redshifting and Doppler shifting the released energy to a stationary observer at infinity. Thus we define the integrated source term

[TABLE]

which agrees with Ressler et al. (2015, \tagform@45), where a fluid-frame heating rate is integrated to find a global quantity measured at infinity. We integrate from $r_{\mathrm{ph}}$ to $r=10$, with the upper bound chosen in order to focus on the shock in the region of the flow in steady state. $Q$ has units of energy per unit coordinate time. The results are plotted in the left panel of Figure 11. At low spins, the heating rate is not strongly affected by disk inclination, implying that whatever nonaxisymmetric shocks are created in these cases are too weak or too limited in extent to significantly alter the flow. For $a=0.5$ the heating rate is $19\%$ larger at $24^{\circ}$ compared to $0^{\circ}$. However, the $a=0.9$ cases show strong tilt dependence, with the $24^{\circ}$ heating rate $89\%$ larger than found at $0^{\circ}$. Thus there is significant extra heating due to this mechanism, but it only applies when spins and inclinations are sufficiently high. This agrees with expectations from Teixeira et al. (2014), who found that no shocks form for $a=0.1$.

The right panel normalizes $Q$ by the time average of $\dot{M}$ over the same time interval over which $Q$ was measured, so the plotted quantity has units of energy per unit mass. This provides a measure of the heating efficiency of the accretion flow (a proxy for the radiative efficiency). The result of increasing initial tilt is that both shock heating and accretion rate increase, and in fact that latter has a steeper dependence on tilt, resulting in lower heating efficiencies. In going from $0^{\circ}$ to $24^{\circ}$, $Q/\langle\dot{M}\rangle$ decreases by $7\%$ for $a=0.5$ and by $26\%$ for $a=0.9$. At higher inclinations the shocks induce more angular momentum transport per unit dissipation. We note however that this trend only applies when considering heating at all latitudes. When dissipation is only measured near the midplane, it increases with inclination faster than $\dot{M}$. That is, the increased heating is more concentrated in the disk midplane than the increased mass flux.

Figure 12 shows the three-dimensional structure of the shock in the $a=0.9$, $24^{\circ}$ inclined case. Both panels are oriented looking down the black hole spin axis, with the disk ascending node vertical in the top half of the figure. The data is averaged over $3000\leq t\leq 4000$. In both panels values are shown along the $\rho=0.05$ isosurface within $r=20$. On the left we show the heating term $p_{\mathrm{gas}}u^{i}\partial_{i}s$ on a logarithmic scale, and on the right we show radial velocity $u^{1}$.

The radial velocity pattern indicates fluid elements are on eccentric orbits, in agreement with Figure 8 but more clearly seen here. This pattern is antisymmetric with respect to reflection through the midplane of the disk. At the disk antinodes, material that is at high latitudes with respect to the disk but near the black hole midplane is moving outward. When these fluid elements reach apoapsis near the disk’s line of nodes, they encounter the shock and are abruptly redirected inward. At the same time we can consider material at the disk antinodes and above the disk but at very high black hole latitudes. This material is moving inward, but when it reaches periapsis near the disk line of nodes it changes radial direction smoothly with no shock. This picture of material shocking where eccentric orbits’ apoapsides cluster (due to radially varying eccentricity) matches that of Fragile & Blaes (2008).

Only one standing shock is visible in Figure 12, with the other one hidden on the underside of the disk, $180^{\circ}$ away in azimuth. The shock locations correspond to where the material abruptly changes from moving outward to moving inward, indicating the shocks are significant causes of angular momentum transport in this simulation.

The relative importance of the shocks can be examined by partitioning the stress-energy tensor into Reynolds and Maxwell components:

[TABLE]

We define $\alpha_{\mathrm{Rey}}$ and $\alpha_{\mathrm{Max}}$ using the appropriate stress-energy terms in the definition

[TABLE]

Here the tensor is evaluated in a Cartesian frame comoving with the average fluid as described in the appendix, §C. The coordinates are oriented such that $z^{\prime\prime}$ points along the angular momentum direction at that radius and $y^{\prime\prime}$ points along the ascending node, and so this definition (corresponding to $\alpha_{\mathrm{acc}}$ in the appendix) captures the radial transport of angular momentum that is aligned with the disk axis and drives accretion. These values are plotted in Figure 13. In the nontrivially tilted cases Reynolds stress is generally larger than Maxwell stress. This difference is especially striking in the innermost parts of the $a=0.9$ tilted simulations, where the shock is strongest.

7 Discussion

We have run ten simulations showing how the same initial magnetized fluid evolves as a function of black hole spin and disk inclination. In all cases with nontrivial inclinations, the disks evolve to a highly warped and twisted steady state that precesses rigidly. While it is clear that torques are being transmitted through the disks, redistributing angular momentum, the shapes are not readily described by a straightforward application of linear bending wave theory. Given that we did not go to the most extreme spins, nor to very high inclinations, we expect many naturally occurring systems to challenge such a theory with their nonlinear behavior.

For these inefficiently cooled, thick accretion flows, there is no observable Bardeen–Petterson effect, nor do the disks break. Rather, the inner parts of the disk smoothly become more inclined than the outer parts (Figure 7). This can be explained by angular momentum transport not being evenly distributed in azimuth, as happens when standing shocks form at the line of nodes, combined with the peculiar shape of inclined circular geodesics in Kerr spacetime. Analytic models that assume azimuthally symmetric coupling between adjacent rings of material can miss this effect.

All of our simulations have a low-density polar region. When there is no spin the relatively little amount of net vertical magnetic flux means this region consists of material falling inward, though an outflow turns on with even moderate spins. At $a=0.5$, this moderate outflow can be entirely disrupted by any inclination in the disk; tilted disks can suppress otherwise moderate outflows. With $a=0.9$ the outflow is strong enough to be formed even in the presence of a tilt. These outflows align with the disk rather than the spin of the black hole, in agreement with Liska et al. (2018). However, the large-scale, average magnetic field in some cases evolves from being initially aligned with the disk to partially aligning with the black hole. Further investigation with simulations that achieve high resolution in the polar regions can help quantify to what extent the velocity and magnetic structure of the outflow depend on the tilt of the disk.

Our results confirm the presence of standing shocks, as found in Fragile & Blaes (2008) and later works. These shocks considerably increase the heating of material and the hydrodynamic angular momentum transport in the inner $5\text{--}10$ gravitational radii of the disk, but only at sufficiently high spin and only when the inclination is larger than about $8^{\circ}$. The effects are more pronounced at higher inclination angles, and they can be the dominant driving mechanism for both dissipation and accretion.

Importantly, these shocks can be better at redirecting fluid elements than stopping them. This can lead to more angular momentum transport and thus accretion per unit dissipation, mimicking the lower heating efficiencies of systems with lower black hole spins. As noted above, however, this trend reverses if the higher latitudes of these accretion flows are neglected. It is possible that a given heating efficiency can correspond to a range of black hole spins. While tilted disk simulations with in-situ radiation will be required to verify these trends, but it is clear that one should be cautious when interpreting possibly tilted flows using untilted models.

In the limit where accretion is dominated by these shocks ($a\gtrsim 0.9$ and $i\gtrsim 16^{\circ}$), it is possible that simple analytic models can be developed to capture the correct behavior. That is, there may be regimes where the complexity of MRI turbulence can be neglected, not because it can be replaced by an isotropic viscosity, but because in reality transport becomes dominated by a mechanism more amenable to modeling. Other processes like the MRI are, however, still needed to drive material close to the black hole in the first place.

Given the presence of standing shocks, we expect to see time variability in accretion, based on the results of Henisey et al. (2012). Here we have focused on large-scale properties that stay constant in time, leaving the examination of the variability present in these simulations for a future work. We also suggest that a parameter survey similar to this be done with more efficiently cooled disks. This can help address whether the lack of disk breaking we see here, where GR and MHD processes are being modeled accurately, extends to thinner disks.

In summary, the magnitude of tilt does not strongly affect dimensionless warping and twisting, nor does it qualitatively change the behavior of inclination and precession angles. It affects the velocity of polar outflows, while having less of an effect on their magnetization. Most importantly, large tilts can increase the efficacy of standing shocks in heating the material and transporting angular momentum.

We thank S. M. Ressler for pointing out the subtleties of entropy and heating in a GR context. This research was supported in part by the National Science Foundation under grants NSF PHY-1748958, NSF AST 13-33612, and NSF AST 1715054; by Chandra theory grant TM7-18006X from the Smithsonian Institution; and by a Simons Investigator award from the Simons Foundation (EQ). This work used the Extreme Science and Engineering Discovery Environment (XSEDE) Comet at the San Diego Supercomputer Center through allocation AST170012.

Appendix A Nodal Precession of Geodesics in Kerr Spacetime

Here we outline the numerical evaluation of the nodal precession rate in Kerr spacetime. The procedure here is exact, as opposed to the well-known Lense–Thirring formula, which is formulated in the linear regime and breaks down for large $a$ and small $r$.

We will work in Boyer–Lindquist coordinates until the end of this section. In this system, the geodesics of massive particles obey the well-known equations (cf. Wilkins, 1972, \tagform@1–3)

[TABLE]

where we define

[TABLE]

and $e$, $l$, and $q$ are constants of the motion.

Consider a particle initially at position $x_{(0)}^{\mu}=(0,r,\pi/2+i,0)$, for fixed radius $r$ and inclination $i$, and with velocity $u_{(0)}^{\mu}=(u_{(0)}^{0},0,0,u_{(0)}^{3})$. We want the orbit to stay at the same radius for all time, and so the radial geodesic equation becomes

[TABLE]

which is a quadratic equation for $u_{(0)}^{3}/u_{(0)}^{0}$ (taking the positive root for prograde motion). This can be used with the velocity normalization

[TABLE]

to find $u_{(0)}^{0}$ and $u_{(0)}^{3}$ independently.

Equations (A1a) and (A1d) can be rearranged, together with (A2a) and (A2c), to yield

[TABLE]

This linear system can easily be solved to find the specific energy and angular momentum constants $e$ and $\ell$. The Carter constant $q$ then follows from (A2e) and knowing $\Theta=0$ from (A1c) (or from (A2c) and (A2d) knowing $R=0$ from (A1b)).

With the constants of motion in hand, we can compute the relative rate of change between the two angular coordinates:

[TABLE]

This formula is similar to \tagform@25 from Wilkins (1972), but we include the spin dependence rather than assume an extremal black hole. The overall negative sign should be taken for prograde motion with $i>0$, as the particle will initially start to move toward smaller $\theta$ as it moves forward in $\phi$. By symmetry we know $\phi$ advances by the same amount in the first quarter of the vertical oscillation as in all other quarters, and so the total azimuthal precession in one vertical period is

[TABLE]

The corresponding change in proper time for the particle is similarly calculated:

[TABLE]

Dividing $\Delta\phi$ by $\Delta\tau$ gives us $\mathrm{d}\phi/\mathrm{d}\tau$, the average nodal precession of particles at a given radius and inclination around a black hole with a given spin. This is easily converted to Kerr–Schild coordinates via

[TABLE]

In particular, we can neglect the difference between Kerr–Schild and Boyer–Lindquist coordinates to the extent our approximation of vanishing radial velocity is appropriate.

For comparison, the Lense–Thirring result is a gravitomagnetic nodal advance of $4\pi ar^{-3/2}$ per orbit. Alternatively, a more accurate estimate of $2\pi(1-\Omega_{\theta}/\Omega)$ can be used (see §B for the definitions of $\Omega$ and $\Omega_{\theta}$), taking into account the retrograde contribution from the black hole’s quadrupole moment. Around a black hole of spin $a=0.9$ at a radius of $r=5$, these formulas predict $\Delta\phi\approx 58.0^{\circ}$ and $\Delta\phi\approx 43.0^{\circ}$, while the exact result ranges from approximately $49.0^{\circ}$ to $49.9^{\circ}$ as the inclination ranges from $8^{\circ}$ to $24^{\circ}$.

Appendix B Dynamical Frequencies for Circular Geodesics in Kerr Spacetime

For reference we give the dynamical frequencies for prograde circular equatorial orbits in Kerr spacetime. All formulas here are the same in both Kerr–Schild and Boyer–Lindquist coordinates, since the two $3{+}1$ metrics induce the same $2{+}1$ metric on any given hypersurface of constant $r$.

The orbital, epicyclic, and vertical frequencies with respect to coordinate time are given by Okazaki et al. (1987) and Kato (1990):

[TABLE]

We can convert these to frequencies with respect to an orbiting particle’s proper time via the Lorentz factor $\gamma\equiv\mathrm{d}t/\mathrm{d}\tau$. As the particle’s velocity is $u^{\mu}=(\gamma,0,0,\gamma\Omega)$, the normalization $g_{\mu\nu}u^{\mu}u^{\nu}=-1$ yields

[TABLE]

With $\gamma$ we can convert $\Omega$ to the orbital frequency per unit proper time $\omega\equiv\mathrm{d}\phi/\mathrm{d}\tau$, and likewise for the other frequencies:

[TABLE]

Plots of these frequencies for the spins considered in our simulations are shown in Figure 14. These enter into the linear warp theory as deviations of epicyclic and vertical frequencies from the orbital frequency. Since the epicyclic frequency necessarily diverges at the ISCO, its deviations are particularly large for the inner parts of the disk, to the point where linear theory may no longer be applicable even for very small inclinations.

Appendix C Definition of Effective Viscosity

In order to define an $\alpha$-viscosity in a way appropriate for tilted disks and GR, we adopt and modify a procedure described in Penna et al. (2013).

We begin with the fluid velocity in Kerr–Schild coordinates $u_{\mathrm{KS}}^{\mu}$. These are easily transformed to Boyer–Lindquist coordinates via

[TABLE]

Knowing the Euler angles $i$ and $\phi_{0}$ at each radius (see §3), we can convert to the disk-aligned frame. Specifically, we need the components

[TABLE]

For completeness, we give all the components of the Jacobian corresponding to (7):

[TABLE]

We next define the mean velocity components

[TABLE]

This averaging only yields meaningful results with tilted components and along rings that vary only in the coordinate $\phi^{\prime}$. We then find the corresponding time components $\bar{u}_{\mathrm{BL}}^{0^{\prime}}$ via

[TABLE]

Note that while the radial and azimuthal components of the mean velocity vary with just $r$ and $\theta^{\prime}$, the time component will in general vary with $\phi^{\prime}$ as well. These four components are then transformed back into the untilted coordinates:

[TABLE]

We can now construct the transformation from untilted Boyer–Lindquist coordinates to an orthonormal frame following an appropriately defined mean flow, as was first done in Krolik et al. (2005). (There is still a degree of arbitrariness in rotating the radial and poloidal directions of this frame, as has been noted by previous authors.) We present the formulas in notation more closely following Penna et al. (2013), correcting one of the normalizations. The transformation components are

[TABLE]

with the definitions

[TABLE]

The velocity at any point can be transformed into this fluid frame according to

[TABLE]

The magnetic field components $b^{i^{\prime\prime}}$ are found in the same way.

We perform one last transformation to Cartesian coordinates:

[TABLE]

(and similarly for the magnetic field). Here the $z^{\prime\prime}$-direction, like the $z^{\prime}$-direction, points in the direction of the disk angular momentum. Thus $x^{\prime\prime}y^{\prime\prime}$-stress corresponds to radial transport of $z^{\prime\prime}$-angular momentum, driving accretion. Similarly the $y^{\prime\prime}$-direction points to the ascending node, so $x^{\prime\prime}z^{\prime\prime}$-stress governs rotation about this axis, aligning the disk with the black hole. Finally, $y^{\prime\prime}z^{\prime\prime}$-stress corresponds to precessing the disk angular momentum about the black hole axis. These considerations justify the definitions

[TABLE]

for the accretion, precession, and alignment effective viscosities. Though the calculation of these relativistic effective viscosities are rather involved, they still broadly agree with the Newtonian equivalents. For example, $\alpha_{\mathrm{pre}}$ is small in magnitude and suffers large cancellations when averaging, as was found with $\alpha_{*}$ derived from $r\theta$-stresses in Sorathia et al. (2013).