THE TIME VARIABILITY OF GEOMETRICALLY THIN BLACK HOLE ACCRETION DISKS. I. THE SEARCH FOR MODES IN SIMULATED DISKS

There is a very extensive literature on oscillations and waves in accretion disks. Here, we focus on just those aspects of the field that turn out to be relevant for the interpretation of our simulation, which have been elucidated most clearly by LP93.

LP93 studied three-dimensional wave propagation in accretion disks ignoring self-gravity. In the case where one ignores vertical motions, they show that radial waves obey the well known dispersion relation:

$$\begin{equation} \omega ^2=\kappa ^2+c_s^2k^2, \end{equation} \tag{ 1 }$$

where c_s is the sound speed (assumed, in this case, to be purely a function of r) and κ is the radial epicyclic frequency given by κ² = 4Ω² + r ∂Ω²/∂r (also see Binney & Tremaine 1987). Here, Ω(r) is the angular frequency of the background Keplerian flow. LP93 proceeded to study the propagation of axisymmetric waves in the case where the atmosphere has a locally isothermal vertical structure. They find two types of waves. There are low-frequency gravity waves for which

$$\begin{equation} 0<\omega <\Omega. \end{equation} \tag{ 2 }$$

There are also high-frequency acoustic waves that have

$$\begin{equation} \omega ^2>(n\gamma +1)\Omega ^2, \end{equation} \tag{ 3 }$$

where n = 0, 1, 2, ... and γ is the adiabatic index. In the special case of purely vertical perturbations, the inequality becomes an equality:

$$\begin{equation} \omega ^2=(n\gamma +1)\Omega ^2. \end{equation} \tag{ 4 }$$

The n = 0 mode corresponds to a bulk vertical displacement of the disk and subsequent vertical oscillation at the vertical epicyclic frequency, which, in the analysis of LP93 and in all analyses performed in this paper, coincides with the orbital frequency. In general, the nth mode has n vertical nodes (i.e., locations where the vertical velocity perturbation vanishes), and is either even or odd depending on whether n is even or odd, respectively.

As discussed in Section 4, our simulations demonstrate the effectiveness with which MHD turbulence excites these local acoustic waves.

As discussed in Section 1, several groups have studied the global oscillation of black hole accretion disks using linear perturbation theory, identifying three classes of the normal mode (g-modes, p-modes and c-modes). Trapped g-modes have received particular attention as a possible source of the HFQPO, although the other families of modes may well be relevant. Here we review some of the basic results of these analyses, following the approach of Nowak & Wagoner (1991, 1992; hereafter NW91 and NW92, respectively). The NW91 and NW92 analyses are not fully relativistic; instead, they employ a pseudo-Newtonian potential. Thus, these analyses can be readily compared with our pseudo-Newtonian simulations. We also note that full general relativistic MHD simulations have typically yielded results for slowly rotating black holes that are very similar to those obtained with pseudo-Newtonian potentials (e.g., Gammie et al. 2003; De Villiers & Hawley 2003).

NW91 and NW92 used a Lagrangian formalism (Friedman & Schutz 1978) and a WKBJ approximation to derive the linearized equations describing perturbations of an inviscid hydrodynamic thin accretion accretion disk about a pure Keplerian background state. They also initially examined the special case of purely radial oscillations and found the standard dispersion relation of disk theory (Equation (1)). Given that the gravity-modified p-modes described by this dispersion relation become evanescent when ω² < κ², global p-modes can be trapped between the ISCO (where κ = 0) and the radius at which ω = κ. In practice, however, the "leaky" nature of the ISCO would seem to make the trapping of these modes ineffective.

The more general case, including perturbations that have vertical as well as radial motions, yields more promising results. NW92 examined the linearized equations describing the behavior of the scalar potential

$$\begin{equation} \delta u\equiv \delta P/\rho, \end{equation} \tag{ 5 }$$

where δP is the Eulerian variation in the pressure. They showed that the linearized equations are approximately separable into radial and vertical equations, with the separation constant being a slowly varying function of r, ϒ(r). The general dispersion relation for these modes becomes

$$\begin{equation} [\omega ^2-\gamma \Upsilon (r)\Omega ^2](\omega ^2-\kappa ^2)=\omega ^2c_s^2k^2\;. \end{equation} \tag{ 6 }$$

Assuming that γϒΩ² > κ², nonevanescent solutions exist for ω² > γϒΩ² (predominantly radial p-modes) or ω² < κ² (predominantly vertical g-modes). Through this analysis, NW92 identified a class of global g-modes that are trapped between two evanescent regions, r < r₋ and r > r₊, where κ(r_±) = ω. In other words, these modes are trapped under the peak of the epicyclic frequency. They principally focused on the m = 0 (axisymmetric) modes, and showed that the mode frequency is only slightly smaller than the maximum radial epicyclic frequency κ_max. Radial harmonics of these modes are very closely spaced. Thus, the inner radius at which the mode becomes evanescent is still a finite distance (and, in plausible settings, several vertical scale heights) from the ISCO. This raises the interesting possibility of having appreciable power in such modes without significant leakage across the ISCO.

A major issue, however, is the effect of the turbulent MHD background state on these modes. The diskoseismic mode frequencies are comparable with the frequencies characterizing the expected MHD turbulent fluctuations (which is very different to the situation in the Sun, for example, where the observed helioseismic modes have frequencies that are four orders of magnitude higher than the turbulent turnover frequency). Thus, an MHD turbulent disk is likely to be a hostile environment for any diskoseismic modes. Furthermore, magnetic forces can lead to a rather gradual transition in flow properties around the ISCO, potentially worsening the leakage of the trapped g-modes. Considering these mode destruction/suppression mechanisms, it is reasonable to suppose that mode survival becomes easier in thinner disks since (1) the ratio of the typical turbulent cell size to the radial extent of the resonant cavity will decrease with disk thickness and (2) there are suggestions that the transition in flow properties around the ISCO is sharper in thinner disks (Reynolds & Fabian 2008; Shafee et al. 2008), thereby producing less mode leakage. This raises the possibility that a g-mode can be sustained against (or even fed by; Nowak & Wagoner 1993) the turbulence in a sufficiently thin disk.

Previously published global MHD disk simulations (e.g., Hawley & Krolik 2001) have modeled flows as thin as h/r ∼ 0.1 and have not reported diskoseismic modes. However, these authors did not conduct a directed search for such modes and hence it is not possible to say that the modes were really not present. Careful examination of local "shearing-box" simulations have indeed failed to find trapped g-modes associated with MHD turbulence (Arras et al. 2006), but this issue has yet to be examined in a global thin-disk setting. Searching for and characterizing these trapped g-modes in global thin-disk simulations will be a major theme of our paper.

From the point of view of explaining HFQPOs, the need for global oscillation modes is diminished if some process does indeed select special radii in the accretion disk. As discussed in Section 1, the discovery of pairs of HFQPOs with small integer ratios has prompted several groups to examine resonance models. In particular, we shall briefly review the parametric resonance model of Abramowicz & Kluźniak (Abramowicz & Kluźniak 2001, 2003; Abramowicz et al. 2002, 2003).

We begin by considering an accretion disk in which the flow deviates only slightly from the Keplerian so that the position of a fluid element (in spherical polar coordinates) is given by

$$\begin{equation} r(t)=r_0+\delta r(t),\quad \theta (t)=\frac{1}{2}\pi +\delta \theta (t),\quad \phi (t)=\Omega t. \end{equation} \tag{ 7 }$$

We have specialized to the case of axisymmetric perturbations. By expanding the resulting equations of motion to third order (and wrapping up the nongravitational forces into two unspecified force functions), one finds a Mathieu-type equation of motion:

$$\begin{equation} \delta \theta _{,\,tt}+\Omega ^2[1+a\cos (\kappa t)]\delta \theta + \lambda \,\delta \theta _{,\,t}=0, \end{equation} \tag{ 8 }$$

where a is a small constant that describes the coupling between the vertical and radial perturbations and λ is a (small) damping constant. Here, we have specialized the equations of Abramowicz et al. (2003; hereafter A2003) to the case where the vertical epicyclic frequency is the same as the orbital frequency. This is valid for a nonspinning black hole and, in particular, the pseudo-Newtonian potential that we use for the simulations in this paper.

One expects a system described by Equation (8) to exhibit a parametric resonance instability when κ = 2Ω/n, where n is a nonzero positive integer. Given that black hole potentials always have κ <  Ω, the smallest value of n for which the resonance condition is obeyed is n = 3, that is, κ = 2Ω/3. This is expected to be the strongest of the set of resonances. It is interesting that the fundamental "test-particle" frequencies at this resonant radius have a 3:2 frequency ratio in agreement with observations of HFQPO pairs. For the Paczynski–Wiita (PW) pseudo-Newtonian potential we use in our simulations (Paczynski & Wiita 1980; hereafter PW),

$$\begin{equation} \Phi =-\frac{GM}{r-2r_g},r_g\equiv \frac{GM}{c^2}, \end{equation} \tag{ 9 }$$

we have

$$\begin{equation} \Omega =\frac{1}{r-2r_g}\sqrt{\frac{GM}{r}} \end{equation} \tag{ 10 }$$

and

$$\begin{equation} \kappa =\sqrt{\frac{GM(r-6r_g)}{r(r-2r_g)^3}}. \end{equation} \tag{ 11 }$$

Hence, the 3:2 resonance occurs at r = 9.2r_g. It must be noted, however, that a full integration of a toy problem by A2003 reveals that higher order effects shift the location of the resonance, making the ratio of the epicyclic frequencies extremely sensitive to the strength of the coupling between the r and θ perturbations. While A2003 suggested that this sensitivity is a strength of the model, allowing application to HFQPO pairs in neutron star systems that do not have simple integer ratios, it inevitably diminishes the power of this model to explain black hole systems. In addition, there is currently no physical model of the coupling between the radial and vertical motions. The simulations described in this paper allow us to assess whether magnetic forces couple these motions in such a way as to drive a parametric resonance instability.

In the remainder of this paper, we construct and analyze numerical simulations of thin disks in order to examine their variability properties, focusing on the presence of local and global modes, as well as parametric resonances. Of course, real accretion disks are believed to require at least an MHD-level description in order to capture the MRI-driven turbulence that transports angular momentum and drives accretion. MHD simulations are addressed in Section 4. However, it is useful to begin with a discussion of ideal hydrodynamic models in order to help isolate the various physical effects present in these complex systems. This will be the focus of this section.

In order to allow us to perform a set of simulations with modest-to-high resolution, we begin with two-dimensional (axisymmetric) hydrodynamic simulations. We expect (and indeed show) that these axisymmetric models are well suited for studying the fundamental m = 0 g-mode. However, once we move to MHD, Cowling's antidynamo theorem (Cowling 1957) leads to qualitatively different behavior in axisymmetric compared with full three-dimensional simulations. Hence, all of the MHD models that we shall describe are performed in three dimensions. Since we will be directly comparing results (e.g., g-mode amplitudes) between the hydrodynamic and MHD simulations, we also need to perform a "bridging" three-dimensional hydrodynamic simulation.

To make the problem tractable despite the severe resolution requirements imposed by the geometrical thinness of the accretion disk, we choose to focus on only the most essential aspects of the physics. From the discussion in Section 2, it is clear that the essential aspect of the relativistic potential that must be captured is the nature of the radial epicyclic frequency (e.g., that it goes to zero at some finite radius and hence produces an ISCO at that radius). In this sense, the PW potential (Equation (9)) is a good approximation for the gravitational field of a nonrotating black hole; its ISCO (at 6r_g) and marginally bound orbit (at 4r_g) are both at the same radius as in the Schwarzschild geometry. We also simplify the simulations by neglecting all radiation processes (radiative heating, radiative cooling, radiative transfer, and the dynamical effects of radiation pressure). In place of a full energy equation, the gas is given an adiabatic equation of state with γ = 5/3.

All simulations are performed in cylindrical polar coordinates (r, z, ϕ). The initial condition consists of a disk with a constant midplane density (ρ₀ = 1) for r > r_isco. The initial density scale height of the disk is assumed to be constant with the radius, implying a sound speed that falls off with radius as approximately r^−3/2. There are two motivations for choosing to model a "constant-h" disk: (1) such a choice is well suited to the cylindrical symmetry of our coordinate grid and (2) according to the standard model of Shakura & Sunyaev (1973), the radiation-pressure-dominated disks of real accreting black holes are likely to maintain an approximately constant scale height in their innermost regions. In more detail, the vertical structure of the disk is given by

$$\begin{equation} \rho (r,z) = \rho _0\,\exp \left(-\frac{z^2}{2h_1^2}\right), \end{equation} \tag{ 12 }$$

$$\begin{equation} p(r,z) = \frac{GMh_2^2}{(R-2r_g)^2R}\,\rho (r,z), \end{equation} \tag{ 13 }$$

where r is the cylindrical radius, z is the vertical height above the disk midplane, and $R=\sqrt{r^2+z^2}$. This corresponds to an isothermal atmosphere, which, when h₁ = h₂, is in vertical hydrostatic equilibrium in the PW potential. In order to give the disk an initial vertical kick, we set h₁ = 1.2h₂ (the values of h₂ for all of our runs are detailed in Table 1). Thus, the initial disk temperature is ∼20% too cold for the density and pressure run, leading to a vertical collapse and bounce of the disk. The initial density is set to zero for r < r_isco. The initial velocity field is everywhere set to

$$\begin{equation} v_\phi =r\Omega =\frac{\sqrt{GMr}}{r-2r_g},v_r=v_z=0, \end{equation} \tag{ 14 }$$

corresponding to rotation on cylinders and pure Keplerian motion for material on the midplane. We impose zero-gradient outflow boundary conditions on both the radial and the vertical boundaries of the simulation, that is, the fluid quantities in the ghost zones are set to the values of the neighboring active zone, and a "diode" condition is imposed on the component of the velocity perpendicular to the boundary that allows outflow but disallows inflow.

Table 1 details our hydrodynamic simulations. We perform four simulations in which strict axisymmetry is imposed at all times (∂/∂ϕ = 0). These axisymmetric simulations were performed using the serial ZEUS-3D MHD code (Stone & Norman 1992a, 1992b) in its pure hydrodynamic axisymmetric mode. In our canonical axisymmetric run (HD2d_1), we set h₂ = 0.3r_g (corresponding to h₂/r = 0.05 at the ISCO). To model the dynamics of the disk in a robust manner, our vertical domain must cover many scale heights; in our canonical run, the vertical domain is z ∈ (−5h₂, + 5h₂). We place 128 uniformly spaced grid cells in this vertical direction, giving 13 cells per scale height. This allows us to resolve hydrodynamic waves with wavelengths of ∼0.5h₂ or greater. The vertical resolution requirements force us to consider a limited radial domain. In the canonical simulation, the radial domain extends from 4r_g (i.e., well within the plunge region) to 16r_g and should correspond to the range of radii where trapped diskoseismic modes or A2003-type resonances occur. Tolerating a grid-cell aspect ratio of ∼2, we place 256 uniformly spaced radial cells in this radial domain. The simulation was evolved for a time 200T_isco, where

$$\begin{equation} T_{\rm isco}\approx 61.6\,GM/c^3, \end{equation} \tag{ 15 }$$

is the orbital period at the ISCO.

In order to test the robustness of the results discussed below to resolution and the limited radial domain, we performed a second simulation (HD2d_1hr) in which the spatial resolution was doubled (i.e., the size of each voxel was halved in both the radial and the vertical dimensions) and the radial domain extended out to 28r_g. We also performed two additional axisymmetric simulations employing the same setup as the canonical simulation but with disks that are half (HD2d_2) and twice (HD2d_3) the thickness (including appropriate modifications to the vertical domain and resolution; see Table 1).

As discussed above, we perform a three-dimensional hydrodynamic simulation in order to aid the later interpretation of the three-dimensional MHD simulations. The setup of our three-dimensional run (HD3d_1) is identical to that for the canonical two-dimensional run except that the computational domain has a ϕ-dimension. To reduce computational expense while capturing the essential physics, we simulate only a Δϕ = 30° wedge of the disk using 32 uniformly spaced grid cells, imposing periodic boundary conditions on the ϕ-boundaries. This three-dimensional simulation was performed using an MPI-parallelized version of ZEUS kindly supplied to us by Eve Ostriker (and similar to the ZEUS-MP code of Vernaleo & Reynolds 2006).

We now discuss the evolution of the axisymmetric hydrodynamic simulations, beginning with the canonical simulation. Starting from the initial condition, the disk undergoes dynamical timescale variability because of pressure gradients, which push matter inside of the ISCO. The strong transients close to the ISCO launch outward-radially directed waves into the disk, which break rapidly to become rolls. This behavior can be seen in Figure 1. These high amplitude transients are short-lived, however, lasting only ∼10T_isco. At subsequent times, the disk settles into a stationary state apart from small amplitude (and decaying) internal oscillations and a very weak accretion stream driven by the numerical viscosity.

Zoom In Zoom Out Reset image size

In order to study the decay of the hydrodynamic fluctuations in a more quantitative manner, we start by computing the quantity

$$\begin{equation} K=\int _{\cal D} \rho v_z^2\,dV. \end{equation} \tag{ 16 }$$

This quantity is particularly well suited to diagnose vertical oscillations of the disk, and will vanish once the hydrodynamic configuration has established a stationary state. In order to diagnose the state of the body of the disk (i.e., side-stepping issues of the outer radial boundary or the plunge region), we choose to compute this integral over a restricted domain ${\cal D}$ consisting of the annulus r ∈ (7r_g, 14r_g). The time dependence of K/max(K) for the canonical axisymmetric simulation (HD2d_1) and its high-resolution counterpart (HD2d_1hr) are shown in Figure 2. The rapid initial decline of K(t) is very similar for these two simulations and corresponds to the strong transients described above. At long times (after about t ∼ 4 × 10³ GM/c³ ≈ 70 T_isco) the behavior of these simulations begins to deviate. HD2d_1 continues to decay in an approximately exponential manner $K(t)\propto e^{-t/\tau _0}$, where τ₀ ≈ 4 × 10³ GM/c³. Superposed on this decay is a distinct oscillation. This corresponds to (twice) the frequency of the trapped g-mode that we shall discuss below. The high-resolution version of this simulation HD2d_1hr shows very similar behavior except that the exponential decay time is longer, τ₀ ≈ 6 × 10³ GM/c³. This suggests that the decay of these small perturbations is due to numerical dissipation, which approximately scales as the square root of the size of the simulation grid cells.

Zoom In Zoom Out Reset image size

We now study the spatio-temporal variability of the disk and, in particular, seek the diskoseismic modes predicted by linear theory. Figure 3 shows the midplane value of v_r on the (r, t)-plane, that is, the value of the function v_r(r, z = 0; t), for run HD2d_1. Note that, outside of the ISCO, the average value is 〈v_r〉 ≪ 0.001c so that this figure essentially plots the fluctuation of v_r from its mean value. At early times, we see strong wave-like disturbances that are generated in the inner parts of the disk (at r ≈ 8r_g) and propagate to both large and small radii. Although the outer radial boundary condition is a zero-gradient outflow, impedance mismatching results in some reflection of these initial strong waves. After these initial transients have died out, it can be seen that the highest amplitude fluctuations are limited to a rather narrow range of radii in the approximate range r=7–9r_g. The fact that these perturbations are essentially vertical on the (t, r)-plane indicates that they are coherent across this radial range, as would be expected if we see a global mode.

Zoom In Zoom Out Reset image size

This temporal variability can be explored in more detail using the power spectral density (PSD), defined as $P(\nu)=\alpha |\tilde{f}(\nu)|^2$, where α is some normalization constant and $\tilde{f}(\nu)$ is the Fourier transform (FT) of the time sequence f(t) under consideration:

$$\begin{equation} \tilde{f}(\nu)=\int f(t)e^{-2\pi i\nu t}\,dt. \end{equation} \tag{ 17 }$$

Note that, for ease of interpretation, most of the PSDs presented in this paper will be in terms of the usual frequency, ν, rather than the angular frequency, ω. Figure 4 shows the PSD of the midplane pressure as a function of r and frequency for HD2d_1 and HD2d_1hr. This is computed using approximately the final half (Δt = 102.4T_isco ≈ 6308 GM/c³) of the simulation in order to avoid the initial strong transients. These PSDs differ in the range r = 12 − 16r_g; run HD2d_1hr has significantly less low-frequency noise than run HD2d_1 in this radial range. We attribute this to effects related to the outer radial boundary, which is at r = 16r_g in run HD2d_1 but substantially further out (r = 28r_g) in run HD2d_1hr.

Zoom In Zoom Out Reset image size

Inside of r = 12r_g, however, the PSDs of these two simulations are very similar and we can trust that neither the resolution nor the outer radial boundary affects the results significantly. We note that simulations in which the initial radial density profile of the disk is truncated before reaching the outer boundary also produce essentially identical results, again giving us confidence that noise infiltration from the outer boundary is not an important issue. In addition to very low frequency noise, the most prominent feature of the inner-disk PSD is a vertical ridge of enhanced power at ν ≈ (4 − 5) × 10⁻³ c³/GM extending from r ≈ 6.5r_g out to r ≈ 9.5r_g. The one-dimensional cuts through this ridge in Figure 5 demonstrate it to have all of the expected properties of a trapped g-mode. First, it exists in a rather narrow range of frequencies, ν ≈ (4 − 5) × 10⁻³ c³/GM, just below the maximum radial epicyclic frequency (κ_max ≈ 5.52 × 10⁻³ c³/GM). Second, it is spatially bounded by the radii at which the mode frequency becomes equal to the radial epicyclic frequency. There is, however, some leaking of the mode down to the ISCO.

Zoom In Zoom Out Reset image size

We can visualize the eigenmode by producing maps of pressure deviations that have been "period-folded" on the period corresponding to the peak power in this mode. More precisely, we use the last half of the simulation to produce maps of the difference between the instantaneous pressure and the time-averaged pressure with a sampling rate of 0.2T_isco. We then sort these maps into 16 phase bins (based on the period corresponding to the peak power in this mode) and average together all maps within a given phase bin. The result is shown in Figure 6. In addition to the g-mode itself, these maps reveal a striking "chevron" pattern at large radii. Time sequences of maps reveal these features to be outward-radially traveling acoustic waves driven by the global g-mode, refracted into the upper layers of the disk atmosphere as they propagate.

Zoom In Zoom Out Reset image size

As a final demonstration that we have properly identified trapped g-modes in our axisymmetric hydrodynamic simulations, we use our set of runs to study the dependence of the mode frequency on the sound speed in the disk. Figure 7 shows the dependence of the mode frequency on the midplane sound speed of the disk (and, hence, disk thickness) at r = 8r_g. As expected from analytical theory (see Appendix A), the difference between the mode frequency and the maximum epicyclic frequency linearly depends on the sound speed. There is, however, a discrepancy in the slope of this linear relationship obtained by the simulations and expected from the analytical theory. Note that the analytical treatment assumes that the gas is strictly isothermal across the whole region of interest, a condition that is clearly violated in the simulated disk. Given the sensitivity of the mode to the radial structure of the disk (as demonstrated by the factor of 2.5 difference in the analytical results between the two pseudo-Newtonian potentials examined in Appendix A), the nonisothermality of the gas in the simulated disk can readily shift the mode frequency away from the analytical value.

Zoom In Zoom Out Reset image size

We begin our analysis of the three-dimensional hydrodynamic simulations (HD3d_1) by examining the decay of hydrodynamic perturbations in our canonical three-dimensional run (HD3d_1) using the quantity K(t) as defined in Equation (16). The K(t) behavior for HD3d_1 is somewhat different from the axisymmetric simulations at late times, reaching a quasi-steady state at K/max(K) ∼ 3 × 10⁻³ with aperiodic fluctuations rather than continuing a ringing exponential decay (Figure 2). Numerical dissipation must be just as effective in HD3d_1 as compared with HD2d_1; hence, these perturbations must be driven by some instability. While it is thought that free Keplerian accretion disks are stable to linear hydrodynamic perturbations (see Balbus & Hawley 1998; Hantao et al. 2006), the presence of the simulation boundaries can introduce true instabilities (associated with reflection from the boundaries) and uncontrolled numerical noise, and these might explain these low-level sustained fluctuations. Visual inspection of density slices in the (r, z)- and (r, ϕ)-planes also suggests that nonaxisymmetric wavelike perturbations interacting with the boundaries are responsible for perturbing the three-dimensional hydrodynamic disk. However, a detailed study of the sustained fluctuations in the three-dimensional hydrodynamic disks is beyond the scope of this paper. Since these perturbations exist at a low level (particularly compared with the MHD turbulence discussed in Section 4; see Figure 2 for a comparison of the sustained fluctuations in the three-dimensional hydrodynamic run compared with the canonical MHD run) and are likely to be driven by our boundaries (hence, are not of astrophysical relevance), they are not of importance for the principal focus of our study.

The principal issue to be addressed here is the impact of the transition to three dimensions on the presence of the trapped m = 0 g-modes in the simulated disks. We might expect the power in these axisymmetric modes to be reduced in the three-dimensional case as the free energy is shared with nonaxisymmetric modes. This is indeed the case. Figure 8 shows the PSD at r = 8r_g of the last Δt = 102.4T_isco of the canonical three-dimensional run HD3d_1. A region of enhanced power is clearly seen in the correct range of frequencies ν ≈ (4 − 5) × 10⁻³ c³/GM to be identified with the axisymmetric g-modes studied in Section 3.3. Furthermore, examination of the radially resolved PSD shows that this region of enhanced power is bounded by the radial epicyclic frequency in precisely the manner expected for trapped g-modes. A comparison of the absolute values of the PSD across the mode does reveal, however, that the mode contains almost an order-of-magnitude less power than in the axisymmetric case.

Zoom In Zoom Out Reset image size

We note the existence of a narrow but large amplitude spike just above a frequency of 7 × 10⁻³ c³/GM in the r = 8r_g PSD of this run. The identification of this feature is not clear; it is not at the frequency of any expected global g- or p-mode. It is, however, confined to a single frequency bin, contains only a small amount of power, and only shows up over a narrow range of radii. It seems likely that this is a noise spike.

We simulate geometrically thin MHD accretion disks by building upon our hydrodynamic computational setup described in Section 3.2. As discussed in Section 3.1, we only consider three-dimensional MHD simulations.

Table 1 details our set of MHD simulations. Most of our discussion will center around run MHD_1 (which we shall refer to as our canonical MHD run). In this run, the hydrodynamic variables are setup as for the canonical hydrodynamic run (see Section 3.2) except that we begin the disk in vertical hydrostatic equilibrium (h₁ = h₂ = 0.3r_g corresponding to h₁/r = h₂/r = 0.05 at the ISCO). An initially weak magnetic field is introduced in the form of poloidal field loops specified in terms of their vector potential A = (A_r, A_z, A_ϕ) in order to ensure that the initial field is divergence free. We choose the explicit form for the vector potential,

$$\begin{eqnarray} &&A_\phi =A_0\,f(r,z)\,p^{1/2}\,\sin \left(\frac{2\pi r}{5h_1}\right),A_r =A_z=0, \end{eqnarray} \tag{ 18 }$$

where A₀ is a normalization constant and f(r, z) is an envelope function that is unity in the body of the disk and then smoothly goes to zero at r = r_isco, r = r_out, and at a location three pressure scale heights away from the midplane of the disk. The use of f(r, z) keeps the initial field configuration well away from either the radial boundaries of the initial disk configuration or the vertical boundaries of the calculation domain. The final multiplicative term produces field reversals with a radial wavelength of 5h. This results in a number of distinct poloidal loops throughout the disk. The normalization constant A₀ is set to give an initial ratio of gas-to-magnetic pressure of β = 10³ in the inner disk. In our canonical MHD simulation, this initial condition is evolved for a duration of 630T_isco (38800 GM/c³) using the MPI-parallelized version of ZEUS described above.

We supplement the ideal MHD algorithms of ZEUS in two ways. First, it is necessary to impose a floor on the density field of 10⁻⁵ times the initial maximum density in order to prevent the numerical integration from producing negative densities. This essentially amounts to a subtle distributed mass source. The density only reaches this floor close to the z-boundary (i.e., many scale heights above and below the disk plane). Second, we implement the prescription of Miller & Stone (2000) to include some effects of the displacement current, principally forcing the Alfvén speed to correctly limit to the speed of light as the magnetic fields become strong. We note that this "Alfvén speed limiter" only plays a role within small patches of the tenuous magnetized atmosphere that forms at large vertical distances above and below the disk; it never plays a significant role in the body of the accretion disk.

Periodic boundary conditions were imposed on the ϕ-boundaries, and zero-gradient outflow boundary conditions were imposed at both the inner and outer radial boundaries (Stone & Norman 1992a, 1992b). However, the choice of the z-boundary condition for this kind of simulation is notoriously problematic. The most physically motivated choice would be a free outflow boundary. However, as described in Stone et al. (1996), field-line "snapping" at these free boundaries can halt such a simulation. Indeed, our own test simulations employing zero-gradient outflow boundary conditions on the z-boundaries were found to be subject to these difficulties, as well as occasional numerical instabilities appearing to result from an interplay of the imposition of the density floor and the free boundary. Furthermore, these tests showed that the tenuous matter high above the disk midplane generally flows slowly across these boundaries at very subsonic and sub-Alfvénic speeds; strictly, this anyways invalidates the use of such boundary conditions (since the flow on the other side of the boundary should be able to act back on the simulation domain).

We adopt the solution of Stone et al. (1996) and choose to employ periodic boundary conditions in the z-directions. While this is obviously unphysical in the sense that matter cannot leave the simulation domain in the vertical direction, it does guarantee mathematically reasonable behavior at the boundary (eliminating numerical instabilities) and, more importantly, appears to have no effect on the dynamics of the accretion disk itself (as diagnosed through comparisons with our vertical-outflow test runs). In order to further isolate the simulated disk from the vertical boundaries, we expand the vertical domain (compared with the canonical hydrodynamic simulation) to z ∈ (−3, 3) (i.e., ±10h₁).

We perform four additional simulations aimed at demonstrating the robustness of the canonical simulation. In run MHD_2, we restrict the vertical domain back to z ∈ (−1.5, 1.5), allowing us to gauge the importance of the location of the z-boundaries. Run MHD_2hr has a setup identical to MHD_2 except that the radial and vertical resolutions are doubled compared with the canonical run (i.e., the voxel size is reduced by a factor of 2 in each of the radial and vertical directions), and the azimuthal domain is doubled to Δϕ = 60° at fixed resolution. Due to the factor of 8 increase in the number of computational cells (and the decrease in the timestep), this run was only integrated for a duration of 85T_isco (5236 GM/c³). While the run time is insufficient to conduct a detailed temporal study, a comparison with run MHD_2 does allow us to investigate the effect of both the radial/vertical resolution and the extent of the ϕ-domain on the turbulent state (see below). As discussed below, we find that the properties of the simulated disks are very similar to these two runs.

Run MHD_3 is similar to MHD_2 except that the outer radius is pushed to r = 28r_g (at fixed resolution), doubling the size of the radial domain. Comparing MHD_2 with MHD_3, we find no evidence that the outer radial boundary affects the inner disk (r < 12r_g) in any way.

We now discuss the evolution and general properties of these MHD disks, centering our discussion around run MHD_1. At very early times (t < 5T_isco), strong hydrodynamic transients dominate the evolution as radial pressure forces drive mass into the region within the ISCO. Similar to the hydrodynamic cases, these transients launch outwardly directed axisymmetric waves that break to form rolls. These strong hydrodynamic transients largely damp away in all but the outermost parts of the disk within 10T_isco. Concurrently, the MRI amplifies the initial magnetic field until the (domain wide) volume-averaged ratio of gas-to-magnetic pressure peaks at 〈β〉 ∼ 5 (at t ≈ 10T_isco), at which point most of the flow has become turbulent. The entire flow (outside of the ISCO) becomes turbulent by t ≈ 20T_isco. Figure 9 displays the density field across a vertical slice through the accretion disk at various times.

Zoom In Zoom Out Reset image size

Between t = 10T_isco and t = 20T_isco, the total magnetic energy in the computational domain declines until 〈β〉 ∼ 20. After this, a quasi-steady state seems to be achieved where the magnetic field generation by the MRI-driven MHD turbulence is balanced by the removal of magnetic field energy due to numerical reconnection and magnetic buoyancy. The fact that buoyancy is playing an important role is revealed by examining time sequences of the strengths of B-field components in (r, z) cross-sections of the disk. One clearly sees highly magnetized structures being generated close to the midplane of the disk, which then propagate vertically away from the midplane.² From this time until the end of the simulation at t = 630T_isco, MHD turbulence and hence accretion is sustained. Over the course of the simulation, the disk loses a little more than 60% of its mass (Figure 10).

Zoom In Zoom Out Reset image size

The total magnetic energy and thermal energy undergo a slow decline as mass is drained out of the simulated disk. During this decline, the volume-averaged plasma-β parameter within the domain remains in the range 〈β〉 ∼ 20–30. However, as expected, the β parameter within the high-density body of the disk is appreciably higher, reaching values of β ≈ 70–100. While this is significantly larger than the β found in global simulations of thicker disks (e.g., Hawley & Krolik 2001, 2002), it is in line with what might be expected for thin accretion disks with a zero net field as diagnosed through local shearing-box simulations both with and without vertical stratification (Stone et al. 1996; Hawley et al. 1996; Miller & Stone 2000).

A generic concern in this class of simulation is the effect of the z-boundaries. Once it achieves its quasi-steady state, our canonical MHD simulation displays a 2–3 order-of-magnitude drop in magnetic pressure, and a 5 order-of-magnitude drop in gas pressure, between the disk and the z-boundary. Thus, the disk boundary seems to be well isolated from the boundary. Further confidence is gained from an examination of run MHD_2 in which the z-boundaries have been brought in from z = ± 3 to z = ± 1.5. Despite the fact that the magnetic pressure now only drops by 1 order of magnitude from the disk to the boundary, all of the results from the canonical MHD run discussed in this paper are reproduced by MHD_2. More precisely, (1) the PSD of the fluid variables (Section 4.3.2) are essentially indistinguishable, failing to show any evidence of global modes but displaying prominent local p-modes, (2) the PSD of the mass accretion rate has a broken power-law form with indices that differ from those found in MHD_1 by Δγ ≈ 0.1 (comparable with the 1σ error bars) and break frequencies that differ by Δ(log ν_break) ≈ 0.05 (again, comparable with the 1σ error bars).

Another generic concern with thin disk simulations is whether the resolution in the vertical direction is adequate. The canonical MHD simulation (MHD_1) has approximately 13 grid cells spanning a vertical range Δz = h₂, implying that we cannot follow any modes with a wavelength smaller than λ ∼ h₂/2. This is just sufficient to follow the fastest growing MRI mode (with wavelength λ ∼ 2πh/β^1/2) if β ≲ 100. To ensure that we have, in fact, achieved adequate resolution, we compare run MHD_2 and its high resolution counterpart, run MHD_2hr. For example, Figure 11 compares the time dependence of the thermal and magnetic energies, as well as the height dependence of the plasma-β parameter. While there is some indication of increased buoyancy-driven escape of magnetic fields from the high-resolution simulation (as is apparent from the higher value of β at intermediate heights), the two simulations generally compare very well. Thus, we conclude that we have achieved adequate numerical resolution.

Zoom In Zoom Out Reset image size

4.3.1. The Importance of Correcting for Secular Changes During the Simulation

The length of our canonical MHD run makes it well suited for a detailed study of temporal variability. In particular, the long stream of simulation data facilitates the construction of PSDs. However, as we now discuss, significant complications arise in the analysis of these MHD simulations as compared with the pure hydrodynamic simulations.

In the case of the hydrodynamic simulations, the background (i.e., unperturbed) flow achieves almost a stationary state once the large transients caused by the initial conditions have died out. In particular, the lack of angular momentum transport within the hydrodynamic disk (other than that due to the very small numerical viscosity) allows the disk to achieve a nonaccreting state. Density and pressure fluctuations about this background state can then be studied via the straightforward construction of the PSD. The constant background level does not contribute to the power spectrum and, hence, the PSD faithfully characterizes the fluctuations of interest.

However, even once the initial transients have been dissipated, MHD disks never achieve this kind of stationary background flow. MHD turbulence leads to continued accretion that depletes mass from the simulated disks. This decline in total mass leads to secular changes (with an approximately exponential form) in the density and pressure of the background flow. Unless corrected for, even a rather slow exponential decay can have a significant influence on the PSD of the pressure or density fluctuations, severely affecting attempts to characterize the properties of the astrophysically relevant fluctuations (i.e., the fluctuations that would be present in the ideal case of a steady-state disk in which the mass was replaced from a large reservoir).

To see this, consider some quantity whose time series f(t) we extract from the simulation (e.g., this could be the midplane density or pressure at some given radius). Let us assume that this can be decomposed as

$$\begin{equation} f(t)=\epsilon (t)[1+g(t)], \end{equation} \tag{ 19 }$$

where g(t) is the time series of the astrophysically interesting fluctuations about some mean state (i.e., 〈g〉 = 0) and (t) is a decay function that describes the secular change in the background state due to the draining of mass from the simulation. It must be noted that the decomposition given in Equation (19) is not completely general; this assumes that the amplitude of the "real" fluctuations is proportional to the mean background value, and that the properties of the fluctuations otherwise remain invariant as the background state decays. This decomposition would be valid if the decay of the simulated disk simply amounted to a gradual decline in the density scale of the simulation while the temperature and velocity structures remained unchanged. We shall refer to such (simulated) disks as density-invariant disks. This does indeed appear to describe our simulated MHD disks (i.e., there is little or no corresponding secular change in disk thickness or characteristic turbulent velocities). In this case, the mass accretion rate will be proportional to the density and the decay will hence have an exponential form, $\epsilon =\epsilon _0\,e^{-t/t_0}$.

In the uninteresting case where there are no fluctuations (g(t) = 0, ∀t), the observed signal is just f(t) = (t), leading to a FT and PSD given by

$$\begin{equation} \tilde{\epsilon }(\omega)=\frac{A_1}{(1/t_0)+i\omega }\,\, {\rm and}\,\,P_\epsilon (\omega)=\frac{A_2}{(1/t_0)^2+\omega ^2}, \end{equation} \tag{ 20 }$$

where A₁ and A₂ are uninteresting normalization constants. Thus, for ω ≫ 1/t₀, the PSD of the exponential decay goes as P(ω) ∝ ω⁻².

In the more interesting case of nonzero fluctuations, the FT of the observed signal is

$$\begin{equation} \tilde{f}(\omega)=\tilde{\epsilon }(\omega)+\int \tilde{\epsilon }(\omega ^\prime)\tilde{g}(\omega -\omega ^\prime)\,d\omega ^\prime, \end{equation} \tag{ 21 }$$

that is, the sum of the FT of the exponential decay with the convolution of the FTs of the decay and the interesting fluctuations. When the fluctuations are small compared with the (decaying) background state, as is the case for the density and pressure, one can see that the PSD of the observed signal will be dominated by 1/ω² associated with the decay. Even in the case where the fluctuations are large compared with the decay, the PSD will still be influenced by the exponential decay due to the convolution term in Equation (21). In particular, regions of the "real" PSD that are steeper than ω⁻² (including the high-frequency wing of any QPO or regions above a high-frequency break) will tend to get filled.

Clearly, we must correct for this decay of the background state, and be cognizant of manifestations of any remaining uncorrected effects of this decay. This procedure plays the same role as the "prewhitening" employed by Schnittman et al. (2006; hereafter SKH). We proceed by dividing the observed time series by a "best-fitting" exponential decay function. Here, the time constant of the exponential decay t₀ is estimated via two methods. First, we can set t₀ by the requirement that the starting and final values of the observed series are equal (f(t = 0) = f(t = T), where t = T corresponds to the end of the time series). This is the procedure adopted by SKH except that they employed a linear form for the decay function. Secondly, we can set t₀ via a least-square fit of an exponential form to the observed time series. These methods give very similar values of t₀ and similar corrected PSDs when applied to density or pressure fluctuations, as we shall now see.

4.3.2. Power Spectra of Fluid Variables in Disk Midplane

We now examine the PSD of the azimuthally averaged fluid variables (velocities, pressure, and density) for the simulated MHD disks, as we did with the hydrodynamic disks in Section 3.3. The top panels of Figure 12 show the PSD for the radial and vertical components of velocity in the midplane of the disk, v^mid_r(r, t) and v^mid_z(r, t), respectively, for a duration lasting Δt = 409.6T_isco starting at t = 100T_isco (i.e., well after all of the initial transients have dissipated and a quasi-steady turbulent state has been established). Note that, outside of the ISCO, these velocity components are themselves first-order fluctuating quantities (i.e., the radial and vertical velocities of the background state are approximately zero) and, within the density-invariant disk assumption, will have characteristic fluctuation amplitudes that remain constant throughout the simulation. Hence, we do not need to correct these PSDs for the effect of the density decay in the simulation.

Zoom In Zoom Out Reset image size

It is readily seen that the PSD of the midplane radial velocity v^mid_r is dominated by the radial epicyclic frequency from r ∼ 7–8r_g out to the outer radial boundary. Inside of r ≈ 7 − 8r_g, the PSD shows broadband power associated with the transition from turbulent to plunging flow.

The PSD of the midplane vertical velocity v^mid_z is quite different and can be interpreted in the light of the local fluid oscillations discussed in Section 2.1. Below the radial epicyclic frequency, there is a broad spectrum of g-modes. Above the orbital frequency, the distinct modes described by Equation (4) become apparent; the n = 0 mode (ω>Ω) and n = 2 mode (ω>2.08Ω) can be picked out as distinct tracks, and there are hints of higher frequency modes as well. The fact that only even-n modes are seen is readily understood given that the odd-n modes have v_z-nodes at the midplane, as is evident from the power deficit centered on the n = 1 frequency in the v_z PSD.

The bottom panels of Figure 12 show the corresponding PSD for the density and pressure from the canonical simulation. The displayed PSDs have been corrected for the decay of the background state by dividing through by an exponential curve that best fits (in a least-square sense) the domain-averaged density or pressure. Very similar results are obtained using an exponential determined from just the end points. The radial banding of the low-frequency noise is largely the effect of the prewhitening procedure; the full (noncorrected) low-frequency power is somewhat greater, and the bands correspond to the residuals between the real decay and the simple exponential model. In addition to this low-frequency noise (which is of secondary interest to the investigation presented in this paper), both PSDs show an enhancement corresponding to the n = 1 vertical p-mode (with ω>1.63Ω). As expected, there is no enhancement in the density or pressure PSDs corresponding to the n = 0 or n = 2 p-modes; these have pressure and density nodes at the midplane.

In stark contrast to the hydrodynamic simulation, the (decay-corrected) midplane pressure PSD does not show the κ-bounded vertical "ridge" on the frequency–radius plane that is characteristic of trapped g-modes. The absence of an excited trapped g-mode is also confirmed by examining the PSD at r = 8r_g (Figure 13, normalized such that a direct comparison can be made with the three-dimensional hydrodynamic results of Figure 8). It is important to note, however, that a trapped g-mode of the strength seen in our three-dimensional hydrodynamic simulation would not stand out from the background turbulence. Thus, while it is apparent that the fundamental trapped g-mode is not strongly excited by the turbulence (e.g., in the way that the local p-modes are), we cannot say whether the g-mode is actively damped by the turbulence.

Zoom In Zoom Out Reset image size

Finally, we note that there is no indication that the parametric resonance instability of Abramowicz & Kluzniak is at work, at least in its simplest form. As discussed in Section 2, this model predicts its strongest resonance at the location where the radial and vertical epicyclic frequencies are in a 2:3 ratio; this occurs at r = 9.2r_g in the PW potential. As shown in Figure 13, the PSD of the pressure fluctuations at r = 9.2r_g shows no structure associated with the local epicyclic frequencies. We have verified that a similar conclusion holds true for the PSD of the other variables. While it is possible that nonlinear effects have shifted the location of the resonance inward from r = 9.2r_g (A2003), we note that the PSDs of Figure 12 show no obvious radius at which the power at the radial epicyclic and the orbital frequencies appears to be locally enhanced.

The previous section addressed the temporal properties of the fundamental fluid variables through the body of the disk. Of course, real observations of accretion disks measure the electromagnetic radiation generated by the accretion flow. While our simulations miss all of the physics relevant for making a first-principles prediction of the observed light curve, it is still instructive to analyze a simple scalar quantity that can be generated from the simulation and may be related to the observed hard X-ray radiation (since it is the hard X-rays that carry the HFQPO signal).

Here we consider one such proxy for the observed light curve, the instantaneous mass accretion rate into the black hole,

$$\begin{equation} \dot{M}=\int _{\partial {\cal R}\,_i}r(-v_r)\rho \,dS, \end{equation} \tag{ 22 }$$

where the integral is performed over the surface $\partial {\cal R}\,_i$ defining the inner radial boundary of the computational domain. Figure 14 shows the raw light curve (i.e., not corrected for the exponential decay of the disk). In the remainder of this section, we address the temporal properties of this light curve. Of particular interest is the presence of breaks or QPOs in the PSD of this light curve. Hence, we need a quantitative approach by which the significance of such features in the PSD can be assessed. We begin by discussing our general statistical approach, which differs from that advocated by SKH.

Zoom In Zoom Out Reset image size

4.4.1. Analysis Method

Our approach to PSD fitting is predicated on the assumption (also made by SKH) that at a given frequency, the power density has an exponential probability distribution; if the mean is p₀, then the probability of measuring a power between p and p + dp is

$$\begin{equation} P(p)dp={1\over p_0}e^{-p/p_0}dp\;. \end{equation} \tag{ 23 }$$

Suppose that we have a model that predicts a power density p_mod,i for frequency bin i, and that in our MHD simulation, we actually observe a power density p_obs,i in that bin. Then the likelihood of the data given the model in that bin is

$$\begin{equation} {\cal L}_i=(1/p_{{\rm mod},i})\exp (-p_{{\rm obs},i}/p_{{\rm mod},i}). \end{equation} \tag{ 24 }$$

The likelihood of the whole power density spectrum given the model is the product of the individual likelihoods, but the log likelihood is typically more useful:

$$\begin{equation} \ln {\cal L}=\sum _i \left[-\ln p_{{\rm mod},i}-p_{{\rm obs},i}/p_{{\rm mod},i} \right]. \end{equation} \tag{ 25 }$$

This is the figure of merit for a given model. It can, therefore, be used for standard tasks, such as parameter estimation and model comparison (e.g., determining if a QPO or break is required by the data).

Note that to maximize the information content, the bin size should be the smallest possible, in this case the frequency resolution of the raw PSD. As a result, this method does not require rebinning to coarser resolution. Our approach, therefore, yields an accurate evaluation of one or more precisely specified models, as opposed to the broader but less sensitive method of trying to detect a signal in a model-independent way.

We apply a Markov Chain Monte Carlo method to search for best fits and establish confidence regions. As discussed in Section 4.3.1, we correct for secular changes. Given the large amplitude fluctuations in the $\dot{M}$ and S_tot curves, the end-point method discussed in Section 4.3.1 is not appropriate. Hence, we employ the least-square method described in Section 4.3.1 in order to determine and then divide out the best-fitting exponential decay.

4.4.2. Results

The mass accretion rate $\dot{M}$ has sufficient intrinsic variability that secular changes in the disk properties have only minor effects on the PSD. This is evident from Figure 15, which compares the PSD and Fourier phase of the raw $\dot{M}$ curve as a function of frequency (left panel) with these after multiplication by the exponential in time that minimizes the overall rms amplitude (right panel). The lack of significant differences suggests that inferences drawn from the PSD are robust. The comparison shown in Figure 15 is for the second quarter of run MHD_1 (approximately t = 150T_isco to t = 310T_isco); analyses of the third and fourth quarters also reveal a lack of sensitivity to the secular decay of the disk.

Zoom In Zoom Out Reset image size

Having established that the secular decay of the disk is unimportant for the PSD of ${\dot{M}}$, we will directly work with the raw time series from run MHD_1, without multiplication by an exponential function. We separately analyze the second, third, and fourth quarters of the data (each encompassing a period of Δt ≈ 160T_isco) to look for trends or stability in the PSD, while discarding the first quarter as potentially biased by initial conditions. The results are summarized in Table 2, where we compare single power-law models for the PSD with models involving a broken power law. We use the method described in the previous section; note that only differences in the log likelihood $\ln {\cal L}$ are important rather than the absolute magnitude of $\ln {\cal L}$.

Table 2. Model Comparisons for PSD of Mass Accretion Rate $\dot{M}$

Data segment	Model	Parameters and Uncertainties	Maximum $\ln {\cal L}$
1	Single PL	Γ = 2.50 ± 0.05	9131
	Broken PL	Γ1 = 0.82 ± 0.09, Γ2 = 2.96 ± 0.07,	9238
		log10(νbreak) = −1.94 ± 0.04
2	Single PL	Γ = 2.70 ± 0.08	6508
	Broken PL	Γ1 = 1.45 ± 0.11, Γ2 = 2.91 ± 0.10,	6609
		log10(νbreak) = −1.65 ± 0.07
3	Single PL	Γ = 2.16 ± 0.04	9561
	Broken PL	Γ1 = 1.47 ± 0.07, Γ2 = 2.89 ± 0.09,	9630
		log10(νbreak) = −1.60 ± 0.04

Note. Γ is the power-law index of a single power-law model for the PSD; P(ν) ∝ ν^−Γ. Γ₁ and Γ₂ are the two indices obtained from a broken power law; $P\propto \nu ^{-\Gamma _1}$ for ν < ν_break and $P\propto \nu ^{-\Gamma _2}$ for ν>ν_break. All error bars are 1 standard deviation.

Download table as:  ASCIITypeset image

From this analysis, we find compelling evidence of a break in the power law characterizing the PSD of ${\dot{M}}$. Compared with a single power law, the broken power fit is better by $\Delta \ln {\cal L}=76-90$; hence, the maximum likelihood ratio is at least exp(76) = 10³³ in all three data segments independently. The break frequency and power-law slopes are consistent between the second and third data segments, but these do not match the first data segment. The break frequencies are within a factor of 2 of, but not consistent with, the orbital frequency at the ISCO, log₁₀(ν_ISCO) = −1.79. Therefore, although there is a clear steepening in the power density spectrum, it is not possible at this point to assign a specific physical meaning to the break. We note that this general form of the PSD, that is, an approximate power law with curvature or a break at frequencies close to the ISCO orbital frequency, has been previously seen in the mass accretion rate of global disk simulations (Hawley & Krolik 2001, 2002).

We see no indications of QPOs in the $\dot{M}$ PSD. Quantitatively, we add a Lorentzian QPO to the broken power-law PSD model in which the QPO centroid frequency, FWHM, and amplitude are allowed to be free parameters, with the one restriction that the quality factor of the QPO must exceed 2. We find that the peak power of any QPO cannot exceed 2% of the continuum power measured at the centroid of the QPO at the 99% confidence level.

In recent years, several groups have reported temporal analyses of MHD disk simulations. Here, we briefly compare our work with some of these published results.

Probably the most relevant previous work is that of Arras et al. (2006; hereafter ABT). These authors perform a local, shearing box MHD simulation of a patch of an accretion disk; this provides a controlled environment in which fluid modes can be characterized. ABT found that the MHD turbulence excites a spectrum of distinct acoustic modes as well as radial epicyclic motions. However, they noted a lack of distinct inertial modes (g-modes) and used this fact to argue against the excitation of trapped g-modes in global accretion disks. Our findings are completely in line with those of ABT, and represent an extension of ABT's conclusions to global simulations of thin accretion disks.

There have been QPOs reported from global simulations. Kato (2004d) performed a global MHD disk simulation in a PW potential and presented an analysis of quantities derived from the mass accretion rate. Through visual inspection of the resulting PSDs from four periods of the (long) simulation, he reported a pair of transient QPOs and a pair of QPOs that are labeled as persistent (although it is not clear that they are present in the PSD of all data segments, and the statistical significance of the features is unclear). The QPOs are attributed to resonances between the vertical and radial epicyclic frequencies, and it is found that these QPO pairs have frequency ratios that are approximately 3:2. We find no evidence of these resonances in our simulations. In another interesting difference, an inspection of the radially resolved PSDs in Kato (2004d) reveal no signs of the local p-modes that seem to feature so prominently in our PSDs. While the reason for these discrepancies is unclear, we do note that the Kato (2004d) simulations have an order-of-magnitude less resolution in both the azimuthal and (more importantly) the vertical directions, although his simulations do have a significantly larger computational domain. It is possible that the Kato (2004d) simulations have failed to adequately resolve the vertical dynamics of the thin disk.

Chan et al. (2006) have also performed global MHD disk simulations in a PW potential using the pseudo-spectral algorithm of Chan et al. (2005, 2006). Their study included a postprocessing of the simulation to include detailed radiative transfer, and was explicitly targeted at understanding the variability (including the large amplitude flaring) of the hot accretion flow around the black hole at the center of the Galaxy. They found that the turbulence of the quiescent flow could only produce a factor of 2 modulations in the observed luminosity. To model the large amplitude flares, they introduced large density perturbations into the flow. After being perturbed, the disk displays a QPO with a frequency equal to the orbital frequency at the magnetosonic point. Given that the Chan et al. study explores a rather different regime of accretion than our present study (i.e., hot, thick accretion flows versus thin, cold accretion flows), it is hard to make a direct comparison.

Finally, SKH have performed detailed analyses of General Relativistic MHD (GRMHD) simulations of disks performed using the code of De Villiers & Hawley (2003). In particular, they have studied a long (6000 GM/c³) simulation of a disk around a Schwarzschild black hole, focusing on the temporal behavior of proxy light curves (rather than the underlying fluid properties discussed in this paper). It is unclear from their analysis whether their simulated accretion disk has excited local p-modes of the type that we find in this current work. SKH did find, though, that a proxy light curve based on radiative transfer through the disk assuming black body emission and free–free absorption displays QPOs with an approximate 3:2 ratio. However, these QPOs are transient, only appearing at certain times and certain viewing inclinations.

As described above, several authors have reported transient QPO pairs from MHD simulations with a frequency ratio of 3:2. It is tempting to interpret these as resonance phenomena. Here, we note that there are several ways that approximate 3:2 ratios can be generated that do not necessarily involve resonances and, hence, one should guard against overinterpreting QPO pairs. Indeed, we must remember that some sources have QPO frequency ratios that are inconsistent with 3:2 (e.g., the 67 Hz and 41 Hz QPOs from GRS 1915+105; see Strohmayer 2001).

For example, consider the local vertical p-modes of the disk (discussed in Section 2.1). At a given radius, the n = 1 vertical pressure mode has a frequency of $\omega _{\rm vert,1}=\sqrt{\gamma +1}\,\Omega$, where Ω is the orbital frequency. For a gas-pressure-dominated disk in which γ = 5/3, ω_vert,1/Ω = 1.63, and for a radiation pressure dominated disk in which γ = 4/3, ω_vert,1/Ω = 1.53. If some unspecified physical process enhances emission from a given ring of the disk, one can imagine a situation where QPOs are generated at the orbital frequency and the n = 1 vertical pressure mode, thereby giving frequency ratios compatible with the measured ratios in several sources. Alternatively, the next lowest vertical mode that has a vertical velocity node in the midplane (and thus the maximum variation of pressure and density there) has a frequency of $\omega _{\rm vert,3}=\sqrt{3\gamma +1}\,\Omega$. As a result, when γ = 5/3, we have ω_vert,3/ω_vert,1 = 1.5, and when γ = 4/3, the ratio is 1.46. Once again, if nature picks out a specific radius and the emission is modulated by the n = 1 and n = 3 modes, we would see QPOs with frequency ratios entirely consistent with the observations.

As another example, suppose that the disk emission is modulated at the vertical and radial epicyclic frequencies and that the emission is distributed in radius according to a standard Page & Thorne (1974) disk. The radial distribution of the emission then peaks close to the radius where the radial epicyclic frequency is a maximum (and, hence, is slowly changing with radius). One might then expect to see a pair of QPOs corresponding to these two epicyclic frequencies. The frequency ratio only weakly depends on the spin parameter, ranging from 1.46 at a/M = 0 to 1.7 at a/M = 0.9.