Epicyclic Motions and Standing Shocks in Numerically Simulated Tilted BlackHole Accretion Disks
Abstract
This work presents a detailed analysis of the overall flow structure and unique features of the inner region of the tilted disk simulations described in Fragile et al. (2007). The primary new feature identified in the main disk body is a latitudedependent radial epicyclic motion driven by pressure gradients attributable to the gravitomagnetic warping of the disk. The induced motion of the gas is coherent over the scale of the entire disk and is fast enough that it could be observable in features such as relativistic iron lines. The eccentricity of the associated fluid element trajectories increases with decreasing radius, leading to a crowding of orbit trajectories near their apocenters. This results in a local density enhancement akin to a compression. These compressions are sufficiently strong to produce a pair of weak shocks in the vicinity of the black hole. These shocks are roughly aligned with the lineofnodes between the blackhole symmetry plane and disk midplane, with one shock above the lineofnodes on one side of the black hole and the other below on the opposite side. These shocks enhance angular momentum transport and energy dissipation near the hole, forcing some material to plunge toward the black hole from well outside the innermost stable circular orbit. A new, extended simulation, which was evolved for more than a full disk precession period, allows us to confirm that these shocks and the previously identified “plunging streams” precess with the disk in such a way as to remain aligned relative to the lineofnodes, as expected based on our physical understanding of these phenomena. Such a precessing structure would likely present a strong quasiperiodic signal.
Subject headings:
accretion, accretion disks — black hole physics — galaxies: active — MHD — relativity — Xrays: stars1. Introduction
The purpose of this paper is to present a more detailed analysis of the numerical simulation results described in Fragile et al. (2007b) (hereafter Paper I). That paper described a global numerical simulation of an accretion disk subject to the magnetorotational instability (MRI) that was misaligned (tilted) with respect to the rotation axis of a modestly fast rotating () Kerr black hole. It was the first such numerical simulation of a tilted disk to fully incorporate the effects of the black hole spacetime and not rely on an ad hoc prescription of angular momentum transport.
A clear identification of a tilted blackhole accretion disk in Nature has yet to be made. Nevertheless, there is reason to believe they may be quite common. In any accretion disk system, the orientation of the disk is set by the net angular momentum of the gas reservoir on large scales. The orientation of the black hole, on the other hand, can either be set by its formation or its evolution. For stellarmass black holes in lowmass binary systems, only the formation is likely to matter (any evolution subsequent to the formation is unlikely to change the orientation of the black hole significantly). Since the formation of the black hole is largely independent of the angular momentum of the gas reservoir, alignment should generally not be expected in this case. For high mass binaries and active galactic nuclei (AGN), on the other hand, a large episode of misaligned accretion could reorient the spin of the black hole (Natarajan & Armitage, 1999). Therefore a more detailed understanding of the evolutionary history of the system, including mergers in the case of AGN, would be needed to know if a tilted configuration is expected.
A tilted disk is subject to differential warping due to the effect of LenseThirring precession. In disks with a large “viscous” stress (parametrized by the Shakura & Sunyaev (1973) parameter) the warping results in the BardeenPetterson configuration (Bardeen & Petterson, 1975; Kumar & Pringle, 1985), characterized by an alignment of the disk with the equatorial plane of the black hole inside some characteristic warp radius. The BardeenPetterson effect has been invoked by a number of authors to explain peculiar observations, such as misaligned jets in AGN (Kondratko et al., 2005; Caproni et al., 2006; Caproni et al., 2007) and Xray binaries (Fragile et al., 2001; Maccarone, 2002).
However, we did not see evidence for the BardeenPetterson effect in our simulation in Paper I. This was not surprising since our simulation was carried out in the low stress regime, with , where is the halfheight of the disk. In this limit warps produced in the disk propagate as waves (Papaloizou & Lin, 1995), rather than diffusively as in the Bardeen Petterson case. Instead of a smooth transition between an untilted disk at small radii and tilted disk at large radii, as with the BardeenPetterson effect, we found the tilt to be a longwavelength oscillatory function of radius.
The tilted simulation in Paper I also showed other dramatic differences from comparable simulations of untilted disks. Accretion onto the hole occurred predominantly through two opposing plunging streams that started from high latitudes with respect to both the blackhole and disk midplanes. These plunging streams also started from a larger radius than the innermost stable circular orbit (ISCO), which is often assumed to represent the inner edge of untilted disks. We interpreted this as being due to the fact that the tilted disk encounters a generalized ISCO surface at a larger cylindrical radius than an untilted disk (Fragile et al., 2007a). In this regard the tilted black hole effectively acts like an untilted black hole of lesser spin.
Because of the fast soundcrossing time in the disk, the torque of the black hole acted globally rather than differentially. Instead of strongly warping the disk, the torque caused the disk to experience global (solidbody) precession. The precession had a frequency of Hz, a value consistent with many observed lowfrequency quasiperiodic oscillations (QPOs). However, this value is strongly dependent on the size of the disk (, where is the outer radius), so this frequency may be expected to be vary over a large range. In the limit of a very large disk or in a case where the disk is strongly coupled to a gas reservoir at large radii, this precession frequency may be expected to drop to zero.
In the present work we expound on some of the features of the tilteddisk simulation which were not described in Paper I. This paper is organized as follows: In §2 we briefly review the details of the simulations presented in Paper I, in particular focusing on models 90h (untilted) and 915h (tilted), which are identical other than the initial tilt of the black hole relative to the disk ( and , respectively), and model 915m, which is a lower resolution version of 915h. In §3, we describe the largescale epicyclic motion that appears in the tilted simulation but is absent in the comparable untilted simulation. This epicyclic motion is driven by radial pressure gradients associated with the stationary bending wave created by the warping action of the gravitomagnetic torque of the black hole (e.g. Nelson & Papaloizou 1999). Next, in §4, we describe a pair of standing shocks that again have no counterparts in untilted simulations. The shocks form roughly along the lineofnodes between the disk midplane and blackhole symmetry plane, one shock on each side of the black hole. We verify that this is not a chance alignment by following the longterm evolution of model 915m and confirming that the shocks remain aligned with the lineofnodes as it precesses with the disk. Finally, in §5, we consider some of the implications of our discoveries, particularly in the context of relativistic iron lines.
2. Numerical Methods
The simulations in Paper I were carried out using the Cosmos++ astrophysical magnetohydrodynamics (MHD) code (Anninos et al., 2005). Cosmos++ includes several schemes for solving the GRMHD equations; in Paper I, the artificial viscosity formulation was used. The magnetic fields were evolved in an advectionsplit form, while using a hyperbolic divergence cleanser to maintain an approximately divergencefree magnetic field. The GRMHD equations were evolved in a “tilted” KerrSchild polar coordinate system . This coordinate system is related to the usual (untilted) KerrSchild coordinates through a simple rotation about the axis by an angle , as described in Fragile & Anninos (2005, see also Fragile & Anninos 2007).
The simulations were carried out on a spherical polar mesh with nested resolution layers. The base grid contained mesh zones and covered the full steradians. Varying levels of refinement were added on top of the base layer; each refinement level doubling the resolution relative to the previous layer. The main simulations, models 90h and 915h from Paper I, had two levels of refinement, thus achieving peak resolutions equivalent to a simulation. To demonstrate convergence, in Paper I we also presented results at higher and lower resolutions. In this paper we also discuss results of model 915m, which used a single level of refinement and had a peak resolution equivalent to a simulation, but was run for significantly longer.
In the radial direction a logarithmic coordinate of the form was used, where is the blackhole horizon radius and is the gravitational radius. The spatial resolution near the blackhole horizon was ; near the initial pressure maximum of the torus, the resolution was . In the angular direction, in addition to the nested grids, a concentrated latitude coordinate of the form was used with , which concentrates resolution toward the midplane of the disk. As a result near the midplane while it is a factor of larger for the fully refined zones near the pole. This grid is shown in Figure 1 of Paper I.
The simulations started from the analytic solution for an axisymmetric torus around a rotating black hole (Chakrabarti, 1985). The initial torus was identical to model KDP of De Villiers et al. (2003), which is the relativistic analog of model GT4 of Hawley (2000). As in model KDP, the spin of the black hole was ; the inner radius of the torus was ; the radius of the initial pressure maximum of the torus ; and the powerlaw exponent used in defining the initial specific angular momentum distribution was . An adiabatic equation of state was assumed, with . The torus was seeded with a weak dipole magnetic field in the form of poloidal loops along the isobaric contours within the torus. The field was normalized such that initially throughout the torus. For the tilted simulations (915h and 915m) the black hole was inclined by an angle relative to the disk (and the grid) through a transformation of the Kerr metric. From this starting point, simulations 90h and 915h were allowed to evolve for a time equivalent to 10 orbits at the initial pressure maximum, , corresponding to hundreds of orbits at the ISCO. Model 915m was evolved for 100 orbital times at the initial pressure maximum. Table 1 summarizes the parameters of all three simulations.
Simulation  Tilt  Equivalent 

Duration 


Angle  Peak  
Resolution  
90h  0  0.9  25  10  
915m  0.9  25  100  
915h  0.9  25  10 
3. Epicyclic Motion
As noted above, LenseThirring precession causes differential warping in tilted blackhole accretion disks. In the thickdisk regime, appropriate for the simulations described in Paper I, warp disturbances are expected to propagate in a wavelike, rather than diffusive, manner. For nearly Keplerian disks, such as those resulting from MRI turbulent simulations, the resulting bending wave is expected to produce large horizontal motion within the disk (Nelson & Papaloizou, 1999; Torkelsson et al., 2000). The radial and azimuthal components of this motion are odd functions of , leading to large vertical shear (, ).
Figures 1 and 2 compare the radial and azimuthal motions of the gas for our tilted and untilted simulation (915h and 90h). Clearly the fluid velocity in the tilted disk is no longer dominated by turbulent motion as it is for untilted disks (De Villiers & Hawley, 2003), but has an ordered sense about it. For instance, in Figure 1a we see that in the righthand side of the image the gas in the upper layers of the disk is moving radially outward, whereas the gas in the lower layers is moving radially inward. The sense of motion is reversed on the lefthand side of the image. Similarly, on the righthand side of Figure 2a the gas in the upper half of the disk is moving slower than the bulk angular velocity, whereas the gas in the lower half is moving faster than this average. Again, the sense of the flow is reversed on the lefthand side of the image. This pattern of motion is reminiscent of epicyclic motion, with the upper and lower halves of the disk executing epicycles that are out of phase with one another. Notice that no such organized motion is apparent in our untilted simulation (Figures 1b and 2b).
We find that the pattern of the epicyclic motion is tied to the orientation of the disk relative to the black hole. As our simulated tilted disk 915m precesses, the velocity pattern of the epicyclic motion changes as represented in Figures 3 and 4, where we show the velocity patterns at and , approximately 1/2 precession period apart. Note that the sense of the epicyclic motion is reversed between the two different evolution times.
Figure 5 shows the epicyclic motion caused by the bending wave from a different perspective. In the figure, a lattice of streamlines starts in the plane. The lattice is centered above and below the symmetry plane of the black hole and the viewer is looking almost down the blackhole spin axis (the blackhole spin axis is tilted away from the viewer’s lineofsight to give a better perspective). As expected from the previous figures of simulation 915h at , the streamlines that begin above the symmetry plane of the disk are initially moving radially outward, whereas those that begin below are moving radially inward. After approximately onequarter of an orbit, the upper streamlines have reached their apocenter and the radial motion changes direction. This is also where material encounters one of the standing shocks, described in the next section, which explains the sudden change in direction of the streamlines. Much of this material then begins plunging toward the hole, ultimately accreting within a couple orbits. Thus, as noted in Fragile & Anninos (2005), highlatitude material is preferentially being drained from the disk.
Figure 5 also reveals that the eccentricity of the particle streamlines depends on the height of the streamline above the disk midplane (the top streamlines are more eccentric that those one row below). This is consistent with the expectation that (Ivanov & Illarionov, 1997), where is the vertical distance measured perpendicular to the midplane of the disk ( is the component of a cylindrical coordinate system (, , ) aligned with each concentric ring of the disk, what is sometimes called the twisting coordinate system).
4. Standing Shocks
The full expression for the eccentricity of the orbit of each fluid element is (Ivanov & Illarionov, 1997)
(1) 
where and are the tilt and twist of the disk, respectively, defined for each concentric ring. In Figure 6 we plot , , and (for computational convenience we have replaced the cylindrical twisting coordinate with the sphericalpolar radius ). It is important to note the increase in with decreasing . Such behavior, with orbits getting more eccentric closer to the black hole, results in a concentration of fluid element trajectories near their respective apocenters (see Fig. 5a of Ivanov & Illarionov, 1997). Notice in Figure 5, for instance, that each of the top four streamlines makes its closest approach to the next streamline out when it is at its own apocenter.
The crowding of fluid element trajectories near their apocenters produces a local density enhancement. This will be most pronounced away from the disk midplane because of the dependence of on . Because the fluid elements are traveling supersonically, this local density enhancement produces a weak shock in the flow. In Figure 7 we identify the standing shocks by the sudden change in the magnitude and direction of the velocity vectors along a roughly linear feature associated with the leading edge of the each plunging stream.
Another way to identify the standing shocks is by plotting the magnitude of the vorticity, indicated by , which increases at a shock. Figure 8 shows an isosurface of overlaid on a density isosurface. Very similar results are obtained if the gradient of the entropy, another indicator of the presence of a shock, is plotted instead. Notice that one of the shock surfaces lies mostly above the chosen isodensity surface, while the other lies below it. This is consistent with the linear dependence of on . We don’t expect the shocks to extend into the disk midplane where epicyclic motion ceases. The location of the shocks in azimuth is also consistent with the apocenters of the fluid elements trajectories, as inferred from Figures 1 – 5.
Based on our physical understanding of the plunging streams identified in Paper I and the standing shocks identified in this paper, we expect both features to maintain a constant orientation vis a vis the lineofnodes between the disk midplane and blackhole symmetry plane. To confirm this we have evolved our “medium” resolution simulation (model 915m of Paper I) for , which is more than a full precession period. Figure 9 shows the plunging streams and standing shocks at 0.1, 0.35, 0.6, and 1.1 precession periods. The orientations of both features closely track the precession of the disk and maintain roughly constant alignments relative to the lineofnodes.
4.1. Shock Effects
The presence of nonaxisymmetric standing shocks in the accretion disk can have important consequences, principally including enhanced angular momentum transport and dissipation. To illustrate the former, in Figure 10 we plot angular momentum residual weighted by the circular orbit value, i.e. we plot as a function of radius, where is the angular momentum of circular orbits with inclinations of and for simulations 915h and 90h, respectively, calculated from the following expression (equation 26 of Paper I)
(2) 
where
(3) 
(4) 
(5) 
and
(6) 
and is the densityweighted shell average of the specific angular momentum. Shell averaged quantities are computed as:
(7) 
where is the surface area of the shell. The data have also been timeaveraged over the interval, , where time averages are defined as
(8) 
The sharp downturn of the specific angular momentum inside in simulation 915h is indicative of extra angular momentum being removed from the flow. This suggests that the standing shock plays a significant role in the transport of angular momentum in the tilted disk.
The shock also enhances the energy dissipation in the disk. This is shown in Figure 11, where we plot the densityweighted shell averages of the fluid entropy in the tilted and untilted disks. Actually, since we are only interested in relative changes in entropy, for simplicity we plot . The significant enhancement in entropy generation inside for the tilted disk indicates the extra dissipation provided by the shock. This upturn nicely coincides with the downturn in Figure 10, further supporting the association of both effects with the shock.
The extra energy dissipation and angular momentum transport significantly affect the inner regions of the disk. In Figure 12, we see that the radial plunging region, defined by the sharp upturn in , begins at a considerably larger radius for the tilted simulation (915h) than for the untilted (90h). This can have important implications for disk observations (Krolik & Hawley, 2002) since most of the radiated energy from a disk comes from just outside the plunging radius. Of special concern is the common use of the inner edge of the disk as a direct indicator of the spin of the black hole. Clearly tilt must be taken into account if one wants to relate the plunging radius of the disk to the spin of the black hole for tilted disks. Figure 12 also compares other characteristic velocities associated with the disks. All of the velocities are density weighted shell averages . The local sound speed is recovered from . The Alfvén speed is
(9) 
We approximate the turbulent velocity as , where
(10) 
is the dimensionless stress parameter. The most notable difference is in the radial inflow velocity , although there is also some indication that the turbulent stress rises faster in the interior of the tilted disk.
4.2. PostShock Magnetic Field
We found previously that magnetic fields remain subthermal everywhere in the disk throughout the simulation (see Figure 11 of Paper I). They, therefore, must not play a significant role in the dynamics of the epicyclic motion, plunging streams, or standing shocks. This statement is further supported by the fact that very similar features were seen in earlier hydrodynamic simulations of tilted blackhole accretion disks (Fragile & Anninos, 2005). Despite the minor role of the magnetic fields in the dynamics, the presence of a standing shock can enhance the strength of the magnetic field. From Fragile et al. (2005) we note that the postshock value of in the Newtonian limit can depend sensitively on the strength of the shock
(11) 
where is the preshock value and we assume the field is oriented perpendicular to the shock normal. For , as in our work, this gives
(12) 
Therefore, for , the postshock will be lower than the preshock value. In Figure 13 we show that there is, indeed, a thin layer of magnetically dominated plasma just behind the shock. To get as Figure 13 starting from (a reasonable guess for the high latitude material), we need , which is consistent with what we see in the simulation. Again, these do not appear to be tremendously strong shocks.
5. Discussion
In this paper we have compared two simulations that are identical to one another in every respect except one: the initial tilt between the blackhole and disk angular momenta. Despite the apparent similarity, we observed remarkable differences in the evolution of these two simulations.
The primary new feature that we describe in the main disk body is a strong epicyclic driving attributable to the gravitomagnetic torque of the misaligned (tilted) black hole. The induced motion of the gas is coherent over the scale of the entire disk. An interesting point about this epicyclic motion that has not been made before is that it could be detectable, for instance in the profile of relativistically broadened iron lines. We have previously pointed out the importance of the iron lines in directly probing tilted accretion disks (Fragile et al., 2005), but at the time we were not in a position to recognize the importance of the epicyclic motion. From Figure 2 we get that the velocities associated with the epicyclic motion represent a significant fraction () of the orbital velocity of the gas. The corresponding shift in a reflection feature such as the iron line should be of a similar magnitude. The interesting thing is to note that an observer viewing the disk from the vantage point of Figure 2a and seeing the “top” of the disk would see a smaller than expected red shift (the gas going away is not moving as fast as expected) and a larger than expected blue shift (the approaching gas is moving faster than expected). An observer viewing the same disk from the same vantage point but seeing the “bottom” of the disk would see exactly the opposite shift. Therefore, depending on the viewing angle, the entire line profile of a misaligned disk could be shifted toward the blue or the red relative to an aligned disk. On the redshifted side this effect might be confused with gravitational redshifting, making it appear that the line is coming from deeper in the potential well than is actually the case.
If the disk actually precesses then there is a clear and simple way to disentangle this effect because the epicyclic motion is phased with the orientation of the disk, as discussed in §3 above. This means that the line shift we are describing would reverse itself twice per precession period, appearing for half a precession period as an overall blueshift and for the other half as an overall redshift. If the integration time of the detector is longer than the precession period, the net effect would generally be to smear or broaden the line. This may provide an alternative explanation for an exceptionally broad line profile such as GX 3394 (Miller et al., 2004), without requiring a high blackhole spin or small disk radius. On the other hand, if the integration time is shorter than the precession period, then one would expect the iron line to vary in phase with changes in the Xray flux (i.e. in conjunction with a low frequency QPO that would correspond to the precession frequency of the disk). Such behavior has been identified in at least one source, GRS 1915+105 (Miller & Homan, 2005).
The second new feature of our tilted disk simulation that we described is a pair of standing shocks, roughly aligned with the lineofnodes between the disk and black hole symmetry planes. This asymmetric shock provides additional angular momentum transport and energy dissipation in the tilted disk, relative to the untilted one. This enhanced dissipation may help compensate for the loss of radiative efficiency (relative to an untilted disk) due to the plunging region starting further out than the equatorial ISCO radius of the spinning black hole. In collisionless accretion flows, such as those thought to be relevant for Sgr A* (e.g. Quataert, 2003), it is conceivable that such shocks could also be sites of particle acceleration. The precession of these shocks might then result in periodic variations of nonthermal radiation.
In considering the magnetic fields in the inner part of the tilted disk, as noted in Paper I, the magnetic fields remain largely subthermal. They do not play an important role in the physics of the epicyclic motion, plunging streams, or standing shocks, although we do find some enhancement of associated with the shocks.
We finish with a figure (14) which provides a visual summary of our findings. We note the overall consistency of our interpretation of the results: 1) The observed epicyclic motion is in agreement with expectations for warped disks in the wavelike propagation limit; and 2) The shocks are located near the apocenters of the epicyclic motion, as one expects when the eccentricity increases with decreasing radius.
Footnotes
 slugcomment: To appear in ApJ
 Radius of initial pressure maximum.
 In units of , the geodesic orbital period at the initial pressure maximum .
References
 Anninos, P., Fragile, P. C., & Salmonson, J. D. 2005, ApJ, 635, 723
 Bardeen, J. M., & Petterson, J. A. 1975, ApJ, 195, L65
 Caproni, A., Abraham, Z., Livio, M., & Mosquera Cuesta, H. J. 2007, MNRAS, 379, 135
 Caproni, A., Abraham, Z., & Mosquera Cuesta, H. J. 2006, ApJ, 638, 120
 Chakrabarti, S. K. 1985, ApJ, 288, 1
 De Villiers, J., & Hawley, J. F. 2003, ApJ, 592, 1060
 De Villiers, J., Hawley, J. F., & Krolik, J. H. 2003, ApJ, 599, 1238
 Fragile, P. C., & Anninos, P. 2005, ApJ, 623, 347
 Fragile, P. C., Anninos, P., Blaes, O. M., & Salmonson, J. D. 2007a, in proceedings of the 11th Marcel Grossmann Meeting on General Relativity (astroph/0701272)
 Fragile, P. C., Anninos, P., Gustafson, K., & Murray, S. D. 2005, ApJ, 619, 327
 Fragile, P. C., Blaes, O. M., Anninos, P., & Salmonson, J. D. 2007b, ApJ, 668, 417
 Fragile, P. C., Mathews, G. J., & Wilson, J. R. 2001, ApJ, 553, 955
 Fragile, P. C., Miller, W. A., & Vandernoot, E. 2005, ApJ, 635, 157
 Hawley, J. F. 2000, ApJ, 528, 462
 Ivanov, P. B., & Illarionov, A. F. 1997, MNRAS, 285, 394
 Kondratko, P. T., Greenhill, L. J., & Moran, J. M. 2005, ApJ, 618, 618
 Krolik, J. H., & Hawley, J. F. 2002, ApJ, 573, 754
 Kumar, S., & Pringle, J. E. 1985, MNRAS, 213, 435
 Maccarone, T. J. 2002, MNRAS, 336, 1371
 Miller, J. M., & Homan, J. 2005, ApJ, 618, L107
 Miller, J. M., et al. 2004, ApJ, 601, 450
 Natarajan, P., & Armitage, P. J. 1999, MNRAS, 309, 961
 Nelson, R. P., & Papaloizou, J. C. B. 1999, MNRAS, 309, 929
 Papaloizou, J. C. B., & Lin, D. N. C. 1995, ApJ, 438, 841
 Quataert, E. 2003, Astronomische Nachrichten Supplement, 324, 435
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
 Torkelsson, U., Ogilvie, G. I., Brandenburg, A., Pringle, J. E., Nordlund, Å., & Stein, R. F. 2000, MNRAS, 318, 47