# Conservative GRMHD Simulations of Moderately Thin, Tilted Accretion Disks

###### Abstract

This paper presents our latest numerical simulations of accretion disks that are misaligned with respect to the rotation axis of a Kerr black hole. In this work we use a new, fully conservative version of the Cosmos++ general relativistic magnetohydrodynamics (GRMHD) code, coupled with an ad hoc cooling function designed to control the thickness of the disk. Together these allow us to simulate the thinnest tilted accretion disks ever using a GRMHD code. In this way, we are able to probe the regime where the dimensionless stress and scale height of the disk become comparable. We present results for both prograde and retrograde cases. The simulated prograde tilted disk shows no sign of Bardeen-Petterson alignment even in the innermost parts of the disk. The simulated retrograde tilted disk, however, does show modest alignment. The implication of these results is that the parameter space associated with Bardeen-Petterson alignment for prograde disks may be rather small, only including very thin disks. Unlike our previous work, we find no evidence for standing shocks in our simulated tilted disks. We ascribe this to the combination of small black hole spin, small tilt angle, and small disk scale height in these simulations. We also add to the growing body of literature pointing out that the turbulence driven by the magnetorotational instability in global simulations of accretion disks is not isotropic. Finally, we provide a comparison between our moderately thin, untilted reference simulation and other numerical simulations of thin disks in the literature.

###### Subject headings:

accretion, accretion disks — black hole physics — relativistic processes## 1. Introduction

Twisted and tilted accretion disks have been investigated for nearly 40 years, starting from the seminal paper of Bardeen & Petterson (1975). Until relatively recently such disks were studied in the framework of different theoretical schemes, where the disk’s tilt and twist have been described by two Euler angles, and , characterizing the inclination and precession of each disk’s ring. Dynamical equations for and have then been derived under various assumptions.

Stationary configurations have also been derived. In their original work, Bardeen & Petterson (1975) proposed such a configuration in which the inclination angle decreased toward the black hole (), and the disk, accordingly, tended to align with the black hole equatorial plane – the so-called Bardeen-Petterson effect. It was shown, however, by Papaloizou & Pringle (1983) that the approach of Bardeen & Petterson (1975) was quantitatively incorrect. Papaloizou & Pringle (1983) and later Kumar & Pringle (1985) found new stationary solutions, which showed that the Bardeen-Petterson effect holds, but the characteristic radial scale of the disk’s alignment decreases when the Shakura-Sunyaev parameter gets smaller, contrary to the original Bardeen & Petterson (1975) claim. Later, a time dependent relaxation to a stationary configuration was investigated by Kumar (1990). This was later generalized by Ogilvie (1999) to take into account non-linear effects in the disk’s inclination angle . However, when the parameter becomes smaller than the ratio of the disk’s semi-thickness to its radius, the Papaloizou & Pringle (1983) theory needs to be modified to take into account both sonic effects (Papaloizou & Lin, 1995) and the effect of Einstein apsidal precession (Ivanov & Illarionov, 1997). In this regime, the stationary configurations do not necessarily show alignment of the disk with the black hole equatorial plane, as was shown by Ivanov & Illarionov (1997) and confirmed by Lubow et al. (2002) and Zhuravlev & Ivanov (2011).

The validity of the Bardeen-Petterson effect in numerical studies is also ambiguous. In the smoothed-particle hydrodynamic (SPH) simulations done by Nelson & Papaloizou (2000), the effect was observed. Those simulations, however, had a large effective viscosity and could not, therefore, test the case of small viscosity, where the analytic theory predicts that the Bardeen-Petterson effect may be invalid. In a similar way, the SPH simulations of Lodato & Pringle (2007); Lodato & Price (2010) confirmed the analytic predictions of Papaloizou & Pringle (1983)^{2}^{2}2Note that Lodato & Price (2010) also found some deviation from the simple linear theory of tilted disks in situations where the gradient of the disk’s tilt is sufficiently large. Deviations of this type are not expected in this work, however, since we deal only with rather small gradients., but were subject to the same assumptions, namely that the physics is adequately described via the normal Newtonian equations plus a torque term, that the “viscosity” is isotropic, and that the Shakura & Sunyaev (1973) parameter, , is constant. Each of these assumptions is now known to be incorrect; hence the need for more realistic numerical simulations.

More recently, Sorathia et al. (2013) performed an MHD simulation of a moderately thin tilted disk and found that the Bardeen-Petterson alignment also holds in their case. However, these authors also modeled the black hole gravitational field via two classical forces: one describing a point-like Newtonian potential and the other describing a gravitomagnetic torque. However, as was discussed by Ivanov & Illarionov (1997), in order to obtain a tilted disk configuration where the Bardeen-Petterson effect is violated, one must take into account relativistic corrections to the point-like Newtonian potential, in particular the Einstein precession of the line of apsides.

On the other hand, in the fully relativistic numerical simulations by Fragile et al. (2007, 2009); Fragile (2009), the Bardeen-Petterson effect has not been observed. However, these authors considered only the case of geometrically thick accretion disks, which are not as easily compared with analytic theory.

In this paper we present the first fully relativistic simulations of tilted thin accretion disks (), where the effective viscosity is induced by the development of the magnetorotational instability, and examine whether the Bardeen-Petterson effect holds. We also include reference simulations of untilted accretion disks; as these are still among only a small selection of such published simulations (others include Shafee et al., 2006; Reynolds & Fabian, 2008; Noble et al., 2009, 2010; Penna et al., 2010), they are worthy of mention, in and of themselves.

In our companion paper (Zhuravlev et al., 2014, hereafter referred to as Paper 2), we compare the numerical results discussed in this paper with a fully relativistic semi-analytic model of a tilted disk, based on the time-dependent formalism introduced in Zhuravlev & Ivanov (2011). The semi-analytic model allows us to clarify the physical meaning of the simulation results. Both the simulations and models show the same qualitative behavior in the disk. Namely, when the disk rotates in the same sense as the black hole (i.e. prograde), Bardeen-Petterson alignment is absent and the disk’s inclination angle actually grows slightly toward the black hole. This can be understood in terms of a standing bending wave. In the case when the disk and black hole rotations are opposite (i.e. retrograde), we observe instead a modest alignment of the disk toward the black hole equatorial plane, which may be interpreted as the Bardeen-Peterson effect. Note that this is the first fully relativistic simulation in which this type of alignment is observed.

In addition to our focus on the effects of tilt, we also consider whether the often-made assumption of an isotropic viscosity is valid. We find that although our simulations yield a reasonable value for , defined in terms of the “horizontal” components of the stress and rate-of-strain tensors, the same is not true for the other components. The anisotropy of the turbulent stresses may be important for many theoretical models of tilted disks, as they have often relied on the assumption of isotropy.

The presentation of our work is organized as follows: In Section 2 we describe the numerical simulations, including an important discussion about resolution. In Section 3, we make comparisons between our untilted simulations and earlier, similar simulations by other authors. Since simulations of moderately thin disks are still relatively new and scarce in the literature, these additional comparisons seem worthwhile. In Section 4, we get to the main result of this paper, which is a presentation of the results of our tilted simulations. Closely connected to our main results is the discussion in Section 5 of the anisotropic effective stress tensor that we get from the magnetorotational instability in our simulations. Finally, we end in Section 6 with some discussion and conclusions.

We use, hereafter, the natural system of units, setting the speed of light and gravatational constant to unity. We adopt the metric signature.

## 2. Numerical Simulations

In this work, we present seven numerical simulations, covering different resolutions, tilted and untilted cases, prograde and retrograde disks, and testing different procedures for introducing the tilt. The main difference between these simulations and our earlier work (Fragile et al., 2007, 2009) is that here we use the new, fully-conservative, high-resolution shock capturing (HRSC) method of the Cosmos++ astrophysical fluid dynamics code (Fragile et al., 2012a), plus an ad-hoc cooling function to control the scale height of the disk (Fragile et al., 2012b). In all of the simulations presented here, we set the target relative scale height to , making these the thinnest tilted disks ever simulated using GRMHD. Our motivation for choosing this value, as well as for the black hole spin, are twofold: first, these values should be reasonably small to facilitate comparison with our semi-analytic model, as discussed in Paper 2; second, we want values that give us a reasonable chance to capture the spatial scales associated with Bardeen-Petterson disks within our computational domain. This also motivates how long we run our simulations, since we wish, at a minimum, to cover the relevant timescales in the inner parts of the disk.

Following our previous GRMHD simulations, we use Kerr-Schild coordinates as the principal coordinate system in our numerical work. To introduce the tilt, we make a rotation of the coordinate system by an amount about the y-axis, as in equation (11) of Fragile et al. (2007). This results in a change of the spherical angles such that initial disk mid-plane always coincides with the equatorial plane of the tilted coordinates, .

### 2.1. Initialization

There are many possible starting configurations that one could consider. Motivated by our own earlier numerical work, we have chosen to initialize the simulations with an axisymmetric torus orbiting around a rotating black hole (Abramowicz et al., 1978), following the procedure described in Chakrabarti (1985). Through the action of turbulence generated by the magnetorotational instability (MRI Balbus & Hawley, 1991), this torus rapidly evolves into an accretion disk. The hope is that the resulting disk will settle to the Bardeen-Petterson solution. Another option would be to try to start the simulation from one of the stationary, tilted-disk solutions in the literature (e.g. Papaloizou & Pringle, 1983; Kumar & Pringle, 1985) and see if the simulation retains this solution. We plan to consider this in future work.

The initial torus configuration is specified by its inner radius (), the radius of its pressure maximum (), the black hole spin (), and a parameter () that is used to define the angular momentum distribution

(1) |

where

(2) |

Here is the covariant 4-velocity of the gas, are the Kerr-Schild metric coefficients in terms of the tilted coordinates, and is a coefficient that is fixed by the requirement that at . Notice that we are not accounting for the effects of tilt in initializing our torus, rather assuming that the equatorial coordinate plane is equivalent to the black hole equatorial plane, which of course it is not whenever . However, since the torus just serves as a convenient starting point and will be unstable because of the MRI anyway, this discrepancy is not of significant concern. Furthermore, as we explain below, most of the simulations start from an untilted () configuration anyway.

In this work we fix , , and . We mention here that we will often use the orbital period of a test particle at , i.e. , as a convenient time unit when discussing our results. The specific internal energy of the gas is fixed by the spacetime and knowledge of (see e.g. Fragile et al., 2007, and references therein for details). This fixes the gas density, , up to an arbitrary constant, where we assume a polytrope. Finally, the gas pressure is fixed by assuming an ideal gas equation of state (EOS),

(3) |

In order to seed the MRI we add a weak initial magnetic field, composed of a single set of poloidal loops that follow the isodensity contours of the torus, such that initially everywhere in the torus, where and are the gas and magnetic pressures, respectively. This field configuration does not contain enough magnetic flux to lead to a magnetically-arrested state over the course of the simulation (McKinney et al., 2012); tilted accretion disks in such a state are considered in McKinney et al. (2013).

In the background region where the torus solution does not apply, we set up a rarefied non-magnetic plasma accreting into the black hole (Komissarov, 2006). The density and internal energy density of this gas are given by and , where is the maximum initial density in the torus. Because there is a free scale in the problem, the absolute numerical value of this density is irrelevant. These profiles also serve as numerical floors on the values of and during the evolution of the simulation. The initial radial velocity profile of this background gas is given by

(4) |

The simulations are carried out on a uniformly spaced grid of spatial coordinates {, , }. All curvature, both real and coordinate, is handled via the metric. This procedure relies on the following transformations between the grid coordinates and and the corresponding Kerr-Schild spatial ones:

(5) |

and (Noble et al., 2010)

(6) |

where is any positive, odd integer (taken to be 9 in our case) and is the amplitude of the nonlinear term. Note that we do not use a cut out region around the pole. Because we use horizon-penetrating Kerr-Schild coordinates, we are able to place the inner radial boundary of the grid inside the black hole event horizon, ( for ), thus isolating it from the physical domain of interest. The inner radial boundary is at , giving 7 zones inside the event horizon for our high-resolution simulations (and 3 for the medium ones). The outer boundary is placed at . We use outflow boundary conditions at both the inner and outer radial boundaries, reflecting boundaries on the pole, and periodic boundaries in the azimuthal direction.

### 2.2. Resolution

Before continuing, we need to address the significant drawback of us trying to simulate the thinnest tilted disks ever using GRMHD – numerical resolution becomes a major issue. Thin accretion disks are challenging in general. Resolving the vertical scale height of the disk while simultaneously maintaining comparable resolution in all three spatial dimensions (a requirement of most grid-based numerical techniques), requires many more zones in the radial and azimuthal directions, roughly . For a thin disk, with , this can be prohibitive. Further, in the case of a titled disk, there are no symmetries to the problem that can be exploited; the full three-dimensional domain must be treated.

To keep the computational expense of this project reasonable, we exploit the fact that we are considering a very small tilt and a moderately thin disk by concentrating most of the zones very near the mid-plane. Even so, our simulations are admittedly very under-resolved with respect to the MRI. In the vertical direction we have approximately 3 zones per MRI wavelength (), and in the azimuthal direction we have about 9 () for our high-resolution simulations, at relevant times () over the radial range (), where

(7) |

, , and the fastest growing MRI wavelength is given by

(8) |

where is the Alfvén speed in the direction of interest; are components of magnetic field defined in the co-moving frame, where is the dual of the Faraday tensor; is the co-moving gas thermal energy per unit of mass; and is the magnetic pressure. The values improve somewhat inside of .

Our modest values mean that the MRI turbulence is not fully developed in these simulations (see Hawley et al., 2011; McKinney et al., 2012; Hawley et al., 2013, for discussions of MRI resolution in numerical simulations). However, we make the following arguments for why our results are still meaningful: First, although this likely means that our effective values for may be smaller than they would be in better resolved simulations, we, nevertheless, have measurable effective viscosities acting in our simulations, and we can meaningfully discuss the response of the disks to them. Second, it has been shown in previous studies (e.g. Fragile & Blaes, 2008; Sorathia et al., 2013) that the main responses within tilted disks are hydrodynamic, not MHD. Therefore, despite the poor resolution of the MRI, the leading order effects in tilted disks are still captured. Additionally, in Section 3 we show that our untilted simulations reproduce many of the key features seen in other simulations of thin accretion disks. Finally, in this and previous studies, we have performed our simulations at different resolutions and confirmed that we find substantially the same behavior at all resolutions.

Nevertheless, there may be one important consequence of our poor resolution of the MRI: the lack of well-developed turbulence may allow the disk to precess more nearly as a solid body than would be the case otherwise. In comparing their hydrodynamic and MHD simulations of Lense-Thirring precession, Sorathia et al. (2013) found that the MHD simulations showed less solid-body precession. They argued that the turbulence in the MHD simulation acted to break up the coupling between neighboring rings that allows for solid-body precession. This may also play a role in how effectively the disk can align with the black hole.

### 2.3. Evolution

Our numerical scheme is explained in more detail in Anninos et al. (2005); Fragile et al. (2012a). Here, we provide only a brief summary of the most relevant details. First, we solve the following set of coupled GRMHD equations:

(9) |

(10) |

(11) |

(12) |

where

(13) |

is the usual MHD stress-energy tensor, is the generalized fluid density, is the total energy density, is the covariant momentum density, is the boosted magnetic field three-vector, , is the generalized boost factor, and is the metric determinant. While most of these variables are defined at the respective cell centers, the magnetic field components, , are staggered, residing at the respective cell faces. This facilitates application of our constrained transport procedure for maintaining a divergence-free field, as described in Fragile et al. (2012a).

The relative scale height of the disk () is controlled through a cooling function of the form (Noble et al., 2009)

(14) |

where is the ratio of the actual to the target temperature

(15) |

is the relativistic orbital frequency, and is the target scale height (set to 0.08 throughout this work).

The basic parameters of each simulation are given in Table 1. The naming convention is in keeping with our previous work and such that the first number indicates the dimensionless black hole spin, in units of tenths; the second number indicates the tilt of the disk, in units of degrees; if followed by the letter “r”, then the simulation is retrograde; the next letter indicates the resolution, either “m” for medium or “h” for high; and finally, the letter “t” indicates a simulation that started with a tilt. Most cases started with an untilted torus. Simulations 10m, 10rm, and 10h are run to their final termination time, , in this untilted configuration. These simulations serve as reference simulations, allowing us to better isolate the effects of tilt. They are also used as background models for the semi-analytic approach, as described in Paper 2. For simulations 110m, 110rm, and 110h, we restart simulations 10m, 10rm, and 10h, respectively, from a time , but with a tilt of . We then allow these simulations to run out to a cumulative evolution time of ( for 110m, for 110rm, and for 110h). This means that the MRI is already established and the disk is accreting before the tilt is introduced. For the final simulation, 110mt, we use the procedure from Fragile et al. (2007), whereby the tilt is introduced immediately after setting up the initial torus. Thus, even the early evolution of this disk includes the effects of tilt. The simulations 10rm and 110rm are our retrograde cases (the black hole spins counter to the disk angular momentum). These are included since retrograde cases are always expected to show a tendency toward (counter-) alignment.

Simulation | Resolution | Tilt () | |||
---|---|---|---|---|---|

10m | 30 | 1033 | 34 | ||

10rm | 33.2 | 1202 | 22 | ||

110m | 30 | 1033 | 34 | ||

110rm | 33.2 | 1202 | 22 | ||

110mt | 30 | 1033 | 22 | ||

10h | 30 | 1033 | 12.5 | ||

110h | 30 | 1033 | 12.5 |

^{a}

^{a}The following parameters remain fixed for all simulations: , , and .

## 3. Comparison with Previous Untilted Simulations

Since these are the first 3D global, GRMHD simulations we have done of thin disks, before considering our tilted simulations, we provide a brief comparison of our untilted ones with earlier GRMHD simulations of comparable thickness, notably the work of Noble et al. (2010) and Penna et al. (2010). Global GRMHD simulations of thin disks are also still relatively new, so it is worthwhile to provide these extra data.

First, to demonstrate that our cooling function is performing as intended, and that we indeed have a moderately thin disk, in Figure 1 we plot a time-averaged, radial profile of the disk scale height, . To estimate it, we use the following expression from Penna et al. (2010)

(16) |

where the scale height is weighted by the square of the density and the integrals are over all cells within a given radial shell. This is a different weighting than we have used in earlier works (e.g. Fragile et al., 2007). This new expression is preferred because it gives a better agreement with the target that we use in our cooling function. Note, however, that the value we get from this estimate of does not necessarily agree quantitatively with other estimates. For example, one might consider (see Section 4.2 for definitions), which comes from the relationship for the vertically-integrated sound speed, . Figure 11 would then suggest .

Figure 1 also includes the time-averaged profile of the effective horizontal viscosity, (see Section 5 for a complete description of how we define ). Here, a normal density weighting is used:

(17) |

This profile of shows similar properties to other such profiles extracted from GRMHD simulations (see, e.g., Penna et al., 2013), in that it has a strong radial dependence at small radii and then asymptotes to a value between at large radii. One difference is that our profile does not turn over and begin decreasing inside the innermost stable circular orbit (ISCO), as it does in Penna et al. (2013), but we must emphasize that we are using a different definition of . Most other work, including our own previous papers, have defined only in terms of the MHD stress tensor in the co-moving frame, normalized to the total pressure, i.e. . This simplification relies on the fact that for purely Keplerian disks, . Here, if we use the same definition, we get profiles that look very similar to profiles from other GRMHD simulations. However, it is not obvious that should be proportional to in the same way for tilted disks. Furthermore, one of our goals in Section 5 is to test whether , a common assumption in most previous analytic work on tilted accretion disks. Therefore, we use the more general expression above.

In Figure 2, we present additional time-averaged radial profiles for simulations 10h and 110h. The top panel presents the time-averaged, radial mass accretion rate,

(18) |

This shows that both simulations have achieved a reasonable inflow equilibrium inside of , with the tilted simulation accreting at a slightly higher rate, consistent with our earlier work (Fragile et al., 2007; Fragile & Blaes, 2008; Generozov et al., 2014). There is, however, a slight rise in inside the ISCO, which we comment on below. The second and third panels show the radial profiles of the shell-averaged gas density, , and density-weighted, shell-averaged gas and magnetic pressures, and , respectively. The fourth panel shows the density-weighted, shell-averaged radial velocity profile, . This profile is almost exactly the same as the comparable profile in Figure 6 of Penna et al. (2010), and follows fairly closely the predictions of Novikov & Thorne (1973, NT73) outside of the ISCO. Finally, the bottom panel shows the density-weighted, shell-averaged specific angular momentum profile, . Again, this plot is very similar to the comparable plot in Figure 14 of Noble et al. (2010). It also tracks fairly closely the model of NT73. The main discrepancy is that continues to drop inside the ISCO, indicating that stresses are still acting to extract angular momentum from the flow, in contrast to the assumption of NT73. This same behavior was noted in Noble et al. (2010) and Penna et al. (2010). The results of Penna et al. (2010) suggest that this discrepancy would become negligible if the disk thickness were to approach zero. Their results also suggest that the discrepancy would be smaller if we had initialized our problem with multiple, alternating islands of poloidal magnetic field loops instead of a single set. Additional radial profiles of , , the surface density , and at different times, for both the high and medium resolution untilted simulations (10h and 10m), are provided in Paper 2. There the profiles serve as inputs to our semi-analytic models, thus allowing the models to capture the effects of both temporal and spatial variability.

Time histories of the mass and angular momentum fluxes at the event horizon are shown in Figure 3. The top panel shows the mass accretion rate, , while the middle panel shows the flux of specific angular momentum at the horizon

(19) |

This value can be compared to the prediction of NT73, which is also shown in Figure 3. We find that our value is 13% less than the prediction. Since NT73 assumes that internal stresses in the disk to vanish at the ISCO, the angular momentum flux through the event horizon should be equal to its value at the ISCO. However, as Figure 2 shows, the specific angular momentum in our simulations continues to fall inside the ISCO, suggesting that stresses are still acting on the gas. Thus, not as much angular momentum reaches the black hole. The bottom panel of Figure 3 includes the “nominal” efficiency, , which represents the total loss of specific energy into the black hole, where

(20) |

Here the time histories of both simulations are consistent with the results of other comparable simulations (for example, see Figure 7 of Penna et al., 2010) and the predictions of NT73. In particular, the nominal efficiency for a thin disk around a slowly rotating black hole is expected to be . Our simulations give values slightly higher than this, consistent with additional energy being liberated from the accreting gas after it passes the ISCO. The discrepancies in and are comparable to those reported in Noble et al. (2010) and Penna et al. (2010) for similar disk thicknesses and field topologies. Again, the discrepancies would likely be smaller had we used smaller, alternating poloidal magnetic field loops in our initialization.

Noting that the mass accretion rate increases slightly inside the ISCO, in contradiction with the assumption of NT73, one could instead compare the angular momentum flux and nominal efficiency using the ISCO value of . For model 10h, the mass accretion rate at the ISCO is 96% of that at the event horizon (0.708 vs. 0.734). Using the ISCO value, instead of the event horizon one, would bring the data up slightly, but it would still lie below the predictions of NT73. For the nominal efficiency, changing in this way would actually drop the simulation results below the NT73 prediction (to about 0.035). However, it is obviously inconsistent to combine different fluxes from different radii in this way; we merely mean to illustrate the sensitivity of these results to the exact fluxes. As for the cause of the increase in inside the ISCO, this is not a systematic effect, but rather reflects time variations in the mass accretion rate that are not completely smoothed out by time averaging that we use.

Finally, in Figure 4, we provide a volume-visualization of the end state of our high-resolution, untilted simulation 10h. This will be useful for making a basic qualitative comparison with our tilted disk results, presented in the next section.

## 4. Results of Tilted Disk Simulations

Figure 5 shows a volume visualization of our high-resolution tilted simulation 110h at the same evolution time as Figure 4. The most remarkable thing to note is how similar the disks appear in both cases, giving the first indication that, even for these much thinner accretion disks (as compared to the ones in Fragile et al. (2007)), there is no sign of Bardeen-Petterson alignment in the prograde case. This is despite the fact that Figure 1 shows that for most of the time-resolved part of the disk (; see Figure 11), , such that one might, conventionally, expect some Bardeen-Petterson-like behavior.

A more quantitative illustration of this statement is given in Figure 6, which plots the tilt, , of model 110h as a function of radius for all times during the simulation (see Paper 2 for details on how we extract from the simulations). Rather than alignment of the disk toward the black hole symmetry plane, this figure shows a tendency for to evolve away from zero, at least at small radii. Furthermore, no part of the disk shows more than a 5% alignment at any time. The same results hold true for model 110m, which was run to a much later stop time. The clear implication is that Lense-Thirring precession does not work to align the disk in this case.

However, Lense-Thirring precession does still cause a twisting of the disk, although again, not in the simple manner that might be expected. Figure 7 shows a spacetime plot for the disk twist, , similar to the previous figure for tilt. This plot reveals two important behaviors in the disk: First, at small radii (), after an initial period of strong differential precession, the twist saturates at a modest value and later relaxes toward a smaller one. Meanwhile, the twist gradually builds up at larger radii, as the rest of the disk “catches up” with the twist of the inner disk. By the end of simulation 110h, the twist front has moved out to about . A roughly similar pattern of strong initial differential precession, followed by more gradual global precession was seen in Sorathia et al. (2013). Generally, somewhat stronger differential precession has been seen in SPH simulations (e.g. Nelson & Papaloizou, 2000), although direct comparisons are difficult as the parameters of SPH and GRMHD simulations are often quite different.

Unlike the prograde cases, the retrograde case, 110rm, does exhibit some tendency to (counter-) align, as shown in Figure 8, although the effect is still not strong. The alignment inside of is only about 10%, and this remains true for most of the duration of the simulation. Although this is not a large value, it is consistent with the predictions of our semi-analytic model, as we show in Paper 2. Using those results as a guide, we predict that thinner disks should exhibit even greater alignment, with full alignment at the inner edge of the disk () likely to occur when , for and . Thus, complete alignment of a disk plane with a black hole symmetry plane may be more likely to occur for retrograde systems than for prograde ones.

All of these results are consistent with the linear analytic theory of tilted disks whenever the effective viscosity is small, (see e.g. Ivanov & Illarionov, 1997; Nelson & Papaloizou, 2000; Lubow et al., 2002; Zhuravlev & Ivanov, 2011). In particular, in this limit, analytic theory predicts alignment is only possible whenever the directions of nodal and apsidal precessions are opposite to each other, while whenever they are directed in the same sense, either growth of the disk inclination angle toward the black hole or radial oscillations of this angle are predicted. Furthermore, in Paper 2 we show that all these results are qualitatively consistent with the predictions of our semi-analytic model, although there are some quantitative discrepancies. We also differentiate which physical effects in the disk are most responsible for our results.

### 4.1. Standing Shocks

One important defining characteristic of the tilted accretion disks in our previous numerical studies was the presence of standing shocks at small radii (Fragile & Blaes, 2008). These shocks were found to be responsible for a number of unique phenomena. They lead to a larger disk truncation radius in simulated tilted disks (Fragile, 2009; Dexter & Fragile, 2011); they may help amplify the inherent variability of accretion disks (Henisey et al., 2012); and they dissipate a significant fraction of the accreted rest mass energy (Generozov et al., 2014). However, we find no evidence for similar shocks in the simulations we present in this paper. One measure of this is to look at the density-weighted shell average of entropy, , as in Figure 9. Since shocks generate entropy, the association of shocks with tilted simulations should manifest itself through an excess of entropy (see, e.g., Figure 14 of Dexter & Fragile, 2011). However, such an excess is not seen in Figure 9; instead the entropy of the untilted and tilted simulations track each other quite closely at all radii. We also point out that the general rise in entropy toward smaller radii, starting from around , is actually associated with the gradual decline of disk density (see Figure 2) and dissipation of magnetic field energy toward the black hole. Other signatures of standing shocks, such as regions of post-shock density enhancement, are also missing in simulation 110h.

The physical explanation put forward in Fragile & Blaes (2008) for the standing shocks is that they are a response of the disk to the crowding of particle trajectories near their apocenters whenever the eccentricity of those trajectories increases significantly toward smaller radii (Ivanov & Illarionov, 1997). The orbital eccentricity, itself, is a manifestation of the epicyclic driving of the gas due to unbalanced radial pressure gradients found at high latitudes in tilted disks. To quantify these statements, we can use the criterion proposed in Ivanov & Illarionov (1997) that shocks will form whenever the eccentricity, , becomes comparable to . This is based upon the assumption that a typical value for the velocity perturbations associated with the eccentricity is and that the sound speed is . For tilted, warped disks, we can estimate the eccentricity at one scale height in the disk as

(21) |

In order to evaluate this expression, we first fit the shell-averaged and with power laws (for this part only, we use the and obtained from Equations (32) and (41), respectively, of Fragile et al., 2007). The result for simulation 110h is shown in Figure 10. Unlike our previous results (see Figure 16 of Dexter & Fragile, 2011), the eccentricity in this simulation is very small and the inequality is satisfied by a wide margin. Thus, the lack of shocks in this simulation is perfectly consistent with our understanding of what generates the standing shocks in our other tilted simulations. In this particular case, it seems the combination of small black hole spin, small tilt, and small scale height prevent the standing shocks from forming.

### 4.2. Timescales

Another way to think of tilted disks is in terms of timescales. For tilted black hole accretion disks, there are a number of relevant ones to consider: 1) the dynamical timescale, ; 2) the radial sound-crossing time, ; 3) the accretion timescale, ; 4) the viscous timescale, ; and 5) the Lense-Thirring precession time, , where is the Lense-Thirring precession frequency. Detailed explanations of these and other timescales are provided in Paper 2. Figure 11 shows a plot of these timescales as a function of radius for simulation 110h. Other than , each of the timescales is generated from density-weighted averages of the relevant fluid variables, , , , , and . From this figure, we see that the duration of our high-resolution simulation 110h, , is enough to cover: the viscous timescales out to ; the Lense-Thirring and accretion timescales out to ; and the two remaining timescales at all radii.

According to the original Bardeen-Petterson picture, a tilted disk should align with the black hole wherever , which is most likely to occur close to the black hole, where is shortest. In Figure 11, we see that this is the case all the way out beyond , though at such large radii both timescales are significantly longer than the duration of our high-resolution simulation, so caution should be used when considering those data. The fact that over much of the disk, would seem to suggest the disk should align with the black hole. However, as mentioned in Section 1, the Bardeen-Petterson argument is known to be quantitatively incorrect. In Paper 2 we introduce a relaxation radius and timescale that provide a better explanation for what is going on here.

The behavior of in our simulations can be understood in terms of the sound-crossing time, . Based on our previous work (Fragile & Anninos, 2005; Fragile et al., 2007), we have argued that if the sound speed in tilted disks is high enough, or equivalently the sound-crossing time short enough, then the Lense-Thirring torque rapidly redistributes itself throughout the body of the disk.^{3}^{3}3This same effect has been discussed in the context of non-relativistic twisted disks in binary systems (see e.g. Papaloizou & Lin, 1995; Larwood et al., 1996). The result is that the disk as a whole experiences rigid-body precession. As we commented in Section 4, the spacetime diagram of in Figure 7 shows that, after an initial period of differential precession, the radial profile of does not change significantly at small radii, while the precession continues to grow gradually at larger radii. Our previous results suggest that once the rest of the disk starts precessing, the whole disk will precess together (e.g. Fragile & Blaes, 2008).

This motivates us to define a total angular momentum vector for the disk (defined over the region ), and see if we find evidence of this vector precessing. The result is shown in Figure 12, where we plot the instantaneous precession angle of the disk angular momentum vector as a function of time. For comparison, we can also estimate the precession timescale, , for the disk (see Paper 2 for details). For small values of the average twist, , when the disk is close to , this timescale is

(22) |

where we have assumed solid-body precession and ignored relativistic corrections and corrections proportional to . Since the surface density is a function of time, this precession timescale is also formally time dependent. However, we have checked that for all times available in our high resolution simulations, the value is close to . From this we get that the average precession angle should depend on time as

(23) |

This prediction is plotted against the simulation data in Figure 12. Although the agreement between the predicted and measured slopes is not great, it is within the implied (order unity) uncertainties in equation (23). These uncertainties come from the fact that the full calculation of involves evaluating a fraction where the numerator is the difference of two large, nearly equal numbers, and the denominator is small (again, see Paper 2 for details). This type of operation is prone to the accumulation of numerical errors. Furthermore, these errors are sensitive to our choice of a small black hole spin and small tilt. To confirm that our method is, nevertheless, sound, we performed a test on a purely hydrodynamic simulation of a tilted disk with a much larger black hole spin () and found that the measured value for precession matched our analytic prediction to better than 1%.

## 5. Turbulent Stress Tensor

Much of the analytic, and even numerical SPH, work on tilted disks relies in some way on the validity of the Boussinesq approximation, where an effective stress tensor describing the action of turbulent motions on the mean flow is postulated to be proportional to the sum of an isotropic tensor, giving an additional contribution to the pressure term, and the rate-of-strain tensor. The latter part describes dissipation of the energy of the mean flow and is proportional to an effective viscosity coefficient, which is, in turn, proportional to the Shakura-Sunyaev parameter . The Boussinesq approximation is based on the assumption that turbulent motions are isotropic in a statistical sense. However, recent work has called this assumption into question (Sorathia et al., 2013; Nauman & Blackman, 2014).

We test the isotropy of the turbulence here for the first time using GRMHD simulations. To do so, we calculate both the stress () and rate-of-strain () tensors in the co-moving frame of the fluid (see Appendix A for how we do this). If the stress tensor were indeed isotropic, then, as we mentioned above, it should be proportional to a multiple of the rate-of-strain tensor plus an isotropic tensor. We would, therefore, expect the ratio between the off-diagonal components of the stress tensor and the same off-diagonal components of the rate-of-strain tensor () to be a constant, which we could use to define an “alpha” for our disk (, where ). This is probably too strong a statement, as the turbulence is more likely to be nearly isotropic, rather than exactly so. In this case, it might be sufficient for us to find that all the have a consistent sign and similar magnitude. Instead, we find that, while behaves as expected for thin accretion disks (with a consistent, positive sign and magnitude in the range 0.01-0.1 inside the disk), the other components, and , vary wildly in magnitude, and are not even consistent in sign.

For example, Figures 13 and 14 show how and vary, respectively, with spatial position over the radial shell at the end of simulation 110h. A plot of is not included, though looks very similar to that of . One sees in Figure 13 that is nearly everywhere consistent in sign, with magnitudes between 0.01 and 1 over most of the shell. This plot is from a single time slice, so, as expected, there is some spatial variability, even in the azimuthal direction. Furthermore, none of the components are shell-averaged in this figure, as they were in Figure 1. The two panels of Figure 14, on the other hand, look very different from Figure 13. The very fact that we have to show two panels, one to represent positive values of and one to represent negative, indicates that, unlike , the signs of the other components of are almost equally likely to be negative as positive. We also see that the magnitudes vary over a larger range, .

These results are broadly consistent with those presented in Sorathia et al. (2013), the only other paper of which we are aware that reported on the isotropy of the stress tensor in global MHD simulations. Note that their only considers the component of the Maxwell stress, and does not include the Reynolds stress. It generally has a magnitude in the range , while they state elsewhere in their paper that the Reynolds stress for this component can often approach . Since our includes both the Maxwell and Reynolds contributions, it seems our results are consistent in both a qualitative and quantitative sense with theirs. About the only difference is that the figures in Sorathia et al. (2013) show much finer spatial structure, as would be expected given their higher resolution.

There are several possible explanations for the apparent anisotropy of the turbulence: First, even assuming that there is a local-in-time correlation between the stress and rate-of-shear tensors, the relationship between them may not be a direct proportionality of their respective components. In an anisotropic medium, a general linear relation could be of the form , where the tensor is not necessarily diagonal. Secondly, there could be non-local-in-time correlations between and . For example, Hirose et al. (2009) found in shearing-box simulations that there could be a time lag between build-ups of and pressure in the disk. This situation can either be described by expressing the relationship between and in terms of a convolution with an appropriate kernel^{4}^{4}4Such convolutions have been discussed by, e.g. Ivanov & Papaloizou (2004), in the context of tidal forcing of a convective star. or using a dynamical equation for obtained from a suitable closure procedure (see, e.g. Kato, 1994; Ogilvie, 2003, for a discussion).

## 6. Conclusions

In this paper we have presented results of the thinnest tilted disk simulations yet done with a GRMHD code. Performing simulations of thinner disks is important, as the behavior of tilted disks is expected to be qualitatively different for thin disks from thick ones. The present simulations are roughly where we expect to begin to see this change of behavior.

Since these are some of the thinnest simulations we have thus far performed, we took this opportunity to compare our untilted reference simulations with earlier such simulations in the literature. Despite our very modest resolution of the MRI, we generally find that our simulations reproduce the expected trends. Furthermore, our numerical results approach the limit of the Novikov & Thorne (1973) disk model over the region where we have achieved a reasonable inflow equilibrium. Evidence of this include the relatively constant value of the specific angular momentum inside the ISCO, as well as an observed accretion efficiency close to the expected value.

The main focus of this work, however, is on tilted disk behavior. Here we were particularly interested in whether or not our simulations would exhibit Bardeen-Petterson-like alignment, as the disk scale height was chosen to be comparable to the effective . There are obviously many ways one could set up such a numerical experiment. Motivated by our own previous work, we chose to initialize the simulations with an isolated, tilted torus, threaded by weak poloidal magnetic fields. This triggered the formation of an MRI turbulent disk that we then cooled, as needed, to maintain the chosen scale height. As with any numerical experiment, one hopes that the results are not strongly sensitive to the initial conditions, although this can not generally be proven without further numerical experiments.

For our prograde simulations, we did not observe Bardeen-Petterson alignment, instead finding that the tilt of the disk grows larger at small radii. This type of behavior is more consistent with the bending wave behavior expected in thicker disks. We tentatively conclude that we are not yet in the parameter space of the Bardeen-Petterson solution. However, we can not be certain that the Bardeen-Petterson solution applies in MRI turbulent disks, as it relies on the assumption of an isotropic viscosity, which is not supported by our simulations. In contrast, we did observe some modest alignment of the disk in our retrograde simulation. In this case, however, even bending wave theory predicts some alignment. Overall, our results appear to be consistent with those analytic and semi-analytic models that fully account for all relativistic effects, a point we explore further in Paper 2.

In this paper, we saw some evidence for solid-body precession of the disk in each simulation, consistent with our earlier work. Such behavior may be important for explaining some types of quasi-periodic oscillations (QPOs) seen in black hole X-ray binaries (Ingram et al., 2009; Motta et al., 2014). However, this should be explored further, as there are significant unanswered questions about this behavior. First, it is not clear how long a tilted disk can continue to precess. In Fragile & Blaes (2008), we confirmed that an isolated thick disk can precess more or less as a solid body for at least one precession period, but to explain QPOs, disks would need to be able to sustain this behavior for very many precession periods. The second question has to do with what influence a large outer disk, and even a binary companion, would have on this precession (see Tremaine & Davis, 2014, for a recent discussion). All of these are issues that should be explored further in future work.

Unlike our earlier work, we saw no evidence for standing shocks in the tilted simulations presented in this paper. We attribute this to the small black hole spin, small tilt, and small disk scale height explored in this paper. Together, these prevent the formation of significant, unbalanced radial pressure gradients in the disk. Such radial pressure gradients are crucial for driving epicyclic motion, leading to non-circular orbits near the black hole. If the eccentricity of these orbits increases toward smaller radii, then there is a crowding of orbits, leading to the formation of standing shocks. However, we find that the eccentricity is very small, and nearly constant, in the thin, slightly tilted disks we explore in this paper, thus preventing the formation of such shocks.

Finally, we confirmed a point that has recently been highlighted by other authors (e.g. Sorathia et al., 2013; Nauman & Blackman, 2014), namely that the turbulent stresses in these simulated disks is definitely not isotropic. At this stage, it is difficult to guess what implications this result has for the theory of tilted and warped accretion disks. We explore this issue to some degree in Paper 2, although a detailed analysis must be postponed until a deeper theoretical understanding of this phenomenon is achieved.

Our results are, of course, tentative until they can be confirmed by further numerical work. Our main concerns are the relatively poor resolution and short evolution times of our simulations. However, remedying these will require significantly greater amounts of computing time than were used here, and so will have to wait.

## Appendix A Calculation of Turbulent Stress Tensor

The components of the stress and rate-of-strain tensors (in the co-moving frame) have the form and , where are the basis vectors describing the local rest frame of the fluid. Note that, to be consistent with previous estimates of the effective viscosities (e.g. Beckwith et al., 2008; Penna et al., 2010, 2013), we calculate the stress-energy and rate-of-strain tensors in Boyer-Lindquist coordinates. Since our calculations are done in Kerr-Schild coordinates, we need to transform the stress energy (13) as follows

(A1) |

where

(A2) |

with . For reference, we reproduce the basis vectors from Beckwith et al. (2008)

(A3) | |||||

(A4) | |||||

(A5) |

where

(A6) | |||||

(A7) | |||||

(A8) |

To calculate the rate-of-strain tensor components, we follow the same procedure, but replace the MHD stress-energy tensor (13) with the covariant rate-of-strain tensor:

(A9) |

where is the projection tensor.

## References

- Abramowicz et al. (1978) Abramowicz, M., Jaroszynski, M., & Sikora, M. 1978, A&A, 63, 221
- Anninos et al. (2005) Anninos, P., Fragile, P. C., & Salmonson, J. D. 2005, ApJ, 635, 723
- Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
- Bardeen & Petterson (1975) Bardeen, J. M., & Petterson, J. A. 1975, ApJ, 195, L65+
- Beckwith et al. (2008) Beckwith, K., Hawley, J. F., & Krolik, J. H. 2008, MNRAS, 390, 21
- Chakrabarti (1985) Chakrabarti, S. K. 1985, ApJ, 288, 1
- Dexter & Fragile (2011) Dexter, J., & Fragile, P. C. 2011, ApJ, 730, 36
- Fragile (2009) Fragile, P. C. 2009, ApJ, 706, L246
- Fragile & Anninos (2005) Fragile, P. C., & Anninos, P. 2005, ApJ, 623, 347
- Fragile & Blaes (2008) Fragile, P. C., & Blaes, O. M. 2008, ApJ, 687, 757
- Fragile et al. (2007) Fragile, P. C., Blaes, O. M., Anninos, P., & Salmonson, J. D. 2007, ApJ, 668, 417
- Fragile et al. (2012a) Fragile, P. C., Gillespie, A., Monahan, T., Rodriguez, M., & Anninos, P. 2012a, ApJS, 201, 9
- Fragile et al. (2009) Fragile, P. C., Lindner, C. C., Anninos, P., & Salmonson, J. D. 2009, ApJ, 691, 482
- Fragile et al. (2012b) Fragile, P. C., Wilson, J., & Rodriguez, M. 2012b, MNRAS, 424, 524
- Generozov et al. (2014) Generozov, A., Blaes, O., Fragile, P. C., & Henisey, K. B. 2014, ApJ, 780, 81
- Hawley et al. (2011) Hawley, J. F., Guan, X., & Krolik, J. H. 2011, ApJ, 738, 84
- Hawley et al. (2013) Hawley, J. F., Richers, S. A., Guan, X., & Krolik, J. H. 2013, ApJ, 772, 102
- Henisey et al. (2012) Henisey, K. B., Blaes, O. M., & Fragile, P. C. 2012, ApJ, 761, 18
- Hirose et al. (2009) Hirose, S., Krolik, J. H., & Blaes, O. 2009, ApJ, 691, 16
- Ingram et al. (2009) Ingram, A., Done, C., & Fragile, P. C. 2009, MNRAS, 397, L101
- Ivanov & Illarionov (1997) Ivanov, P. B., & Illarionov, A. F. 1997, MNRAS, 285, 394
- Ivanov & Papaloizou (2004) Ivanov, P. B., & Papaloizou, J. C. B. 2004, MNRAS, 353, 1161
- Kato (1994) Kato, S. 1994, PASJ, 46, 589
- Komissarov (2006) Komissarov, S. S. 2006, MNRAS, 368, 993
- Kumar (1990) Kumar, S. 1990, MNRAS, 245, 670
- Kumar & Pringle (1985) Kumar, S., & Pringle, J. E. 1985, MNRAS, 213, 435
- Larwood et al. (1996) Larwood, J. D., Nelson, R. P., Papaloizou, J. C. B., & Terquem, C. 1996, MNRAS, 282, 597
- Lodato & Price (2010) Lodato, G., & Price, D. J. 2010, MNRAS, 405, 1212
- Lodato & Pringle (2007) Lodato, G., & Pringle, J. E. 2007, MNRAS, 381, 1287
- Lubow et al. (2002) Lubow, S. H., Ogilvie, G. I., & Pringle, J. E. 2002, MNRAS, 337, 706
- McKinney et al. (2012) McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083
- McKinney et al. (2013) —. 2013, Science, 339, 49
- Motta et al. (2014) Motta, S. E., Casella, P., Muñoz-Darias, T., et al. 2014, ArXiv e-prints
- Nauman & Blackman (2014) Nauman, F., & Blackman, E. G. 2014, MNRAS, 441, 1855
- Nelson & Papaloizou (2000) Nelson, R. P., & Papaloizou, J. C. B. 2000, MNRAS, 315, 570
- Noble et al. (2009) Noble, S. C., Krolik, J. H., & Hawley, J. F. 2009, ApJ, 692, 411
- Noble et al. (2010) Noble, S. C., Krolik, J. H., & Hawley, J. H. 2010, ApJ, 711, 959
- Novikov & Thorne (1973) Novikov, I., & Thorne, K. 1973, in Black Holes, ed. C. DeWitt & B. DeWitt (New York: Gordon and Breach), 343–450
- Ogilvie (1999) Ogilvie, G. I. 1999, MNRAS, 304, 557
- Ogilvie (2003) —. 2003, MNRAS, 340, 969
- Papaloizou & Lin (1995) Papaloizou, J. C. B., & Lin, D. N. C. 1995, ApJ, 438, 841
- Papaloizou & Pringle (1983) Papaloizou, J. C. B., & Pringle, J. E. 1983, MNRAS, 202, 1181
- Penna et al. (2010) Penna, R. F., McKinney, J. C., Narayan, R., et al. 2010, MNRAS, 408, 752
- Penna et al. (2013) Penna, R. F., Sa̧dowski, A., Kulkarni, A. K., & Narayan, R. 2013, MNRAS, 428, 2255
- Reynolds & Fabian (2008) Reynolds, C. S., & Fabian, A. C. 2008, ApJ, 675, 1048
- Shafee et al. (2006) Shafee, R., McClintock, J. E., Narayan, R., et al. 2006, ApJ, 636, L113
- Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
- Sorathia et al. (2013) Sorathia, K. A., Krolik, J. H., & Hawley, J. F. 2013, ApJ, 777, 21
- Tremaine & Davis (2014) Tremaine, S., & Davis, S. W. 2014, MNRAS, 441, 1408
- Zhuravlev & Ivanov (2011) Zhuravlev, V. V., & Ivanov, P. B. 2011, MNRAS, 415, 2122
- Zhuravlev et al. (2014) Zhuravlev, V. V., Ivanov, P. B., Fragile, P. C., & Teixeira, D. M. 2014