THE TORQUING OF CIRCUMNUCLEAR ACCRETION DISKS BY STARS AND THE EVOLUTION OF MASSIVE BLACK HOLES

Maser-emitting nuclear accretion disks are unique, clean probes of the environment near MBHs. The first discovered and best-studied maser disk, NGC 4258, is also the thinnest and most Keplerian (to better than 1%; Maloney 2002) nuclear disk yet discovered. Radio observations of H₂O maser emission from the edge-on disk (Greenhill et al. 1995) reveal the disk morphology, gas velocities, and accelerations, and allow accurate measurements of the MBH mass (M_• = 3.7 × 10⁷ M_☉), the distance to the host galaxy (D = 7.2 ± 0.3 Mpc), and the spatial extent of the maser region (R_a = 0.13 pc to R_b = 0.26 pc) (Herrnstein et al. 1996, 1999). Together with optical and X-ray observations, these also constrain the disk's accretion rate ($10^{-6}<\dot{M}<10^{-4}\,M_{\odot }\,\mathrm{yr^{-1}}$; Neufeld & Maloney 1995). The NGC 4258 maser disk shows a clear O(10°) warp (Herrnstein et al. 1996). It is noteworthy that observations also indicate possible warps in the maser disks of Circinus (Greenhill et al. 2003) and NGC 3393 (Kondratko et al. 2008).

Here, the stellar cusp around the MBH is modeled as a single mass power-law cusp of stars, N_⋆(<  r) = N_h(r/r_h)^{3 − γ}, which extends inward from the radius of influence at r_h = GM_•/σ²_b, where the velocity dispersion of the galactic bulge σ_b can be estimated from the empirical M_•/σ relation (Ferrarese & Merritt 2000; Gebhardt et al. 2003; Shields et al. 2003), and contains N_h = μ_hM_•/M_⋆ stars of mass M_⋆ each, where μ_h ∼ O(1). We model NGC 4258 with r_h = 7 pc and assume μ_h = 2 (the formal value for an MBH-less singular isothermal distribution), M_⋆ = 1 M_☉, and unless stated otherwise, γ = 1.75 (a dynamically relaxed cusp; Bahcall & Wolf 1976). For the purpose of simulating the torques on the maser disk of NGC 4258, only stars in the radial range 0.01 pc to r_h are considered. The maser disk model is described in detail in Section 2.2.

Resonant relaxation (Rauch & Tremaine 1996; Rauch & Ingalls 1998; Hopman & Alexander 2006) is a rapid angular momentum relaxation mechanism that operates in potentials with a high degree of approximate symmetry, which restricts orbital evolution (e.g., fixed ellipses in a nearly Keplerian potential or fixed planar rosettes in a nearly spherical potential). In such potentials the fixed, orbit-averaged, stellar mass distribution exerts a constant residual torque on a test mass, which persists over a coherence time t₀ as long as perturbations due to deviations from the perfect symmetry remain small. The accumulated change in angular momentum J over t₀ then becomes the "mean free path" in J-space for the noncoherent random walk phase on timescales longer than t₀. When t₀ is long, the mean free path is large and the random walk is rapid. The efficiency of RR is determined by the nature of the physical process that perturbs the symmetry and limits t₀ (e.g., the Keplerian symmetry is perturbed far from the MBH by the potential of the stars and near it by relativistic precession). For circular orbits in a near spherical potential, such as those of gas streams in an accretion disk, the RR torques can only change the orbital orientation, but not the eccentricity, i.e., J(r) = J_c(r) = const (Gürkan & Hopman 2007), where $J_{c}=\sqrt{GM_{\bullet }r}$ is the maximal (circular) angular momentum at radius r. On timescales t < t₀, the direction of the angular momentum vector of a circular orbit of radius r changes coherently due to the residual forces by stars on the same scale,

$$\begin{eqnarray} &&w(r)\equiv \left|\Delta J_{\perp }(r)\right|/J_{c}(r)=\beta _{\perp }\sqrt{N_{\star }(r)}(M_{\star }/M_{\bullet })t/P(r), \nonumber \\ &&\vspace*{-4pt} \end{eqnarray} \tag{ 1 }$$

where N_⋆ is the number of stars within r and P is the orbital period. The numeric prefactor β_⊥ ∼ O(1) may depend somewhat on the parameters of the cusp and can be determined from N-body simulations, where it is measured to be $\beta _{\perp }\gtrsim \sqrt{2}$ (Eilon et al. 2009; G. Kupi & T. Alexander 2012, in preparation); $\beta _{\perp }=\sqrt{2}$ is adopted here. The "warp factor" w (0 ⩽ w ⩽ 2) corresponds to a tilt in $\mathbf {J}$ by cos i = 1 − w²/2.

In the limit where the disk mass M_d is negligible compared to the MBH mass M_•, the coherence time is set by the randomizing effect of RR itself on the torquing stars ("self-quenching"). In this case, the coherence time for the torque on a circular orbit of radius r by stars on the same scale is

$$\begin{equation} t_{{\rm sq}}(r)\sim A_{{\rm sq}}(M_{\bullet }/M_{\star })P(r)/\sqrt{N_{\star }(r)},\vspace*{-4pt} \end{equation} \tag{ 2 }$$

where A_sq is an O(1) numeric prefactor. Preliminary analysis of N-body simulations indicates that A_sq = 1.0 ± 0.1 (G. Kupi & T. Alexander 2012, in preparation); A_sq = 1 is adopted here. The warp factor grows over the self-quenching time to $w(t_{{\rm sq}})=\beta _{\perp }A_{{\rm sq}}\simeq \sqrt{2}$, that is, quenching occurs when the accumulated tilt in the orbital direction grows to i ≃ π/2. An approximate correction for the coherence time when the back reaction of the disk on the stars cannot be neglected is introduced in Section 1.3.

We argue here, based on an approximate analysis of scales, that purely gravitational interactions between the stars and the disk, rather than hydrodynamic ones, typically dominate the torquing of the disk; that around lower mass MBHs (≲  10⁷ M_☉) the disk warps in response to the stellar torques, while around more massive MBHs it changes its orientation as a nearly rigid body; and that stellar torques are expected to warp the maser disk of NGC 4258 by O(10°).

Direct stars/disk interactions. Stars crossing the disk exert a torque on it by direct hydrodynamic interaction. To assess whether this effect competes with the purely gravitational RR torques, we estimate the torque exerted on the disk in reaction to the ram pressure the disk exerts on the stars. The residual torque by N_⋆ stars of mass M_⋆ and radius R_⋆, on randomly oriented orbits, in a volume of size R enclosing a disk of mass M_d, is

$$\begin{equation} T_{\mathrm{ram}}\sim RF_{\mathrm{ram}}\sim R\sqrt{N_{\star }}\rho v^{2}\pi r_{X}^{2}\,, \end{equation} \tag{ 3 }$$

where ρ is the disk density, v² ∼ GM_•/R, r_X = max (R_⋆, r_B) is the effective cross-section radius of the star, where r_B = 2GM_⋆/(v² + c²_s) ≃ 2(M_⋆/M_•)R is the Bondi radius in the limit c_s ≲ O(10 km s⁻¹) ≪ v ∼ O(10³ km s⁻¹), which holds for a cold accretion disk close to the MBH. For M_⋆ = 1 M_☉, R_⋆ = 1 R_☉, and for M_• = 3.7 × 10⁷ M_☉, R = 0.26 pc and H/R = 0.002 (the values for NGC 4258; Herrnstein et al. 2005), r_X = R_⋆. The gas density can be estimated by ρ ∼ M_d/(πR²2H), where H is the disk's exponential scale height. We further assume that the cold thin disk is close to its self-gravity stability limit (as implied by the non-smooth morphology of many observed galactic maser disks, Maloney 2002; Braatz & Gugliucci 2008, and is predicted to be the general case in active galactic nucleus (AGN) disks; Goodman 2003), M_d/M_• ∼ H/R. The ram torque is then

$$\begin{equation} T_{\mathrm{ram}}\sim R\sqrt{N_{\star }}\frac{M_{d}}{2\pi R^{2}H}\frac{GM_{\bullet }}{R}\pi R_{\star }^{2}=\frac{1}{2}\sqrt{N_{\star }}\frac{GM_{\bullet }^{2}}{R}\left(\frac{R_{\star }}{R}\right)^{2}\,, \end{equation} \tag{ 4 }$$

while the residual RR torque on the disk by the stars is

$$\begin{equation} T_{{\rm RR}}\sim \sqrt{N_{\star }}\frac{GM_{\star }M_{d}}{R}\,. \end{equation} \tag{ 5 }$$

The RR torques affect the disk continuously over an orbital time, whereas the ram torque is applied only twice per orbit, when the star crosses the disk. The ratio between the durations the torques are active is Δt_ram/Δt_RR ∼ 4H/(2πR), and therefore the ratio of accumulated changes in the disk angular momentum per orbital time is

$$\begin{equation} \frac{\Delta J_{\mathrm{ram}}}{\Delta J_{{\rm RR}}}=\frac{T_{\mathrm{ram}}}{T_{{\rm RR}}}\frac{\Delta t_{\mathrm{ram}}}{\Delta t_{{\rm RR}}}=\frac{1}{\pi }\frac{M_{\bullet }}{M_{\star }}\left(\frac{R_{\star }}{R}\right)^{2}\,. \end{equation} \tag{ 6 }$$

For the parameters of NGC 4258, ΔJ_ram/ΔJ_RR ∼ 10⁻⁷. Even if AGN disks are limited by gravitational instability to radii as small as R ∼ 2000r_g (r_g = GM_•/c²), as argued by Goodman (2003), then ΔJ_ram/ΔJ_RR ∼ 0.02(M_•/10⁶ M_☉)⁻¹, and RR torquing dominates for all MBHs with mass M_• ≳ few × 10⁴ M_☉. The torquing of a thin disk by ram pressure is therefore negligible, and the torquing effects of the stars on the disk are well approximated by pure gravitational interactions.

Stellar torques versus disk angular momentum transport. The effect of the RR torques on the disk geometry depends on the ratio between the stellar torques and the rate of angular momentum diffusion in the disk. When the internal viscous torques are small, differential torquing across the disk results in warps, which persist on the warp diffusion timescale t_warp ∼ (R/H)²P/2πα₂, where α₂ ∼ O(1) is the dimensionless viscosity parameter associated with the vertical viscosity (Equation (12)). The necessary condition for the disk to be Keplerian, N_⋆ ≪ M_•/M_⋆, and the condition for the warp to persist against diffusion, t₀ < t_warp, provide together a soft upper limit to the mass range of MBHs where disk warping is efficient:

$$\begin{equation} M_{\bullet }<M_{\star }\left(\frac{R}{H}\right)^{4}\frac{f_{N}}{4\pi ^{2}\alpha _{2}^{2}}\,, \end{equation} \tag{ 7 }$$

where N_⋆ = f_NM_•/M_⋆ is the number of stars out of the N_h ∼ O(M_•/M_⋆) in the MBH radius of influence that are close enough to the disk to affect it appreciably. Conversely, when t₀ > t_warp, the disk flattens out rapidly, and the stellar torquing results in an overall evolution in the inclination of the disk as a nearly rigid body. The upper limit for warping depends quite sensitively on the assumed parameters, especially the disk's aspect ratio. For R/H = 250, f_N = 0.1, M_⋆ = 1 M_☉, and α₂ = 1, it is M_• ∼ O(10⁷ M_☉).

Residual stellar angular momentum. A random density perturbation on the spatial scale r has a Poisson magnitude $\sqrt{N_{\star }(r)}$ and carries angular momentum of the order of $J_{N}\sim \sqrt{N_{\star }(r)}M_{\star }\sqrt{GM_{\bullet }r}$. This is the maximal angular momentum that can be transferred to the disk from that scale. The torque exerted by the stars on the disk will lead to an equal and opposite torque by the disk on the stars. By the time the disk has been torqued by J_N, the reaction force on the stars will disperse the perturbation, and another, uncorrelated one, will take its place. The back-reaction thus introduces a new coherence timescale for the perturbation, t_react(r) such that $\left|\Delta \mathbf {J}(t_{\mathrm{react}})\right|=J_{N}(r)$, which may be shorter than the self-quenching coherence timescale, t_sq, and thus limit the action of the torques on the disk.

We estimate t_react here by approximating the disk, which extends between R_α and R_β, as a thin ring of mass M_d and mass-weighted mean radius $R_{d}=2\pi M_{d}^{-1}\int _{R_{\alpha }}^{R_{\beta }}R^{2}\Sigma dR$, and evaluating the stellar torque on the disk in the limits of small and large spatial scales. The lever arm on the gas disk is R_d. When r ≫ R_d, the magnitude of the force by the stars on the disk is $\sqrt{N_{\star }(r)}GM_{\star }/r^{2}$, and so ΔJ(t)/J_c(R_d) = β_⊥A_sq(R_d/r)^1/2t/t_sq(r). Conversely, when r ≪ R_d, the force on the disk is $\sqrt{N_{\star }(r)}GM_{\star }/R_{d}^{2}$, and so ΔJ(t)/J_c(R_d) = β_⊥A_sq(r/R_d)^3/2t/t_sq(r). In shorthand notation,

$$\begin{equation} \Delta J(t)/J_{c}(R_{d})=\beta _{\perp }A_{{\rm sq}}(r/R_{d})^{1/2-\Theta }t/t_{{\rm sq}}(r)\,, \end{equation} \tag{ 8 }$$

where Θ = sign(r − R_d). The back-reaction timescale is then

$$\begin{eqnarray} t_{\mathrm{react}}(r) &=& \frac{\sqrt{N_{\star }(r)}}{\beta _{\perp }A_{{\rm sq}}}\frac{M_{\star }}{M_{d}}\left(\frac{r}{R_{d}}\right)^{\Theta }t_{{\rm sq}}(r)\nonumber\\ &=& \frac{1}{\beta _{\perp }}\frac{M_{\bullet }}{M_{d}}\left(\frac{r}{R_{d}}\right)^{\Theta }P(r), \end{eqnarray} \tag{ 9 }$$

which increases monotonically with r. The coherence time is t₀(r, R_d) = min (t_sq, t_react). Since typically $t_{\mathrm{react}}(r,R_{d})/ t_{{\rm sq}}(r,R_{d}) = \sqrt{N_{\star }(r)}(M_{\star }/M_{d})/\sqrt{2}<1$ for r ∼ R_d, where the stellar torques on the disk are most effective, the back-reaction time is what sets the limit on the coherence time. The change in the angular momentum of the disk over the back-reaction time is

$$\begin{equation} w_{\mathrm{react}}(r)=\sqrt{N_{\star }(r)}\frac{M_{\star }}{M_{d}}\sqrt{\frac{r}{R_{d}}}\,. \end{equation} \tag{ 10 }$$

For the NGC 4258 cusp parameters (Section 1.1) and a disk mass of M_d ∼ 3000 M_☉ (Martin 2008; BA09) with r ∼ R_d ≃ 0.16 pc, the predicted overall change in disk orientation is ∼27°, and the differential change across the inner and outer edges of the masing region, the observed warp angle ω, is an order unity fraction of that (BA09). We therefore anticipate that RR warping can lead to the O(10°) warp that is observed in NGC 4258.

The equation governing the evolution of the angular momentum surface density, L, of a thin accretion disk under the influence of internal and external torques is given in the limits of a Keplerian velocity field, no azimuthal modes, and diffusive warp propagation (α₁ > H/R) by (Papaloizou & Pringle 1983; Pringle 1992; Ogilvie 1999)

$$\begin{eqnarray} \frac{\partial \mathbf {L}}{\partial t} & = & \frac{1}{R}\frac{\partial }{\partial R}\left[3R\frac{\partial }{\partial R}\left(\nu _{1}L\right)\boldsymbol{\ell }+\frac{1}{2}\nu _{2}RL\frac{\partial \boldsymbol{\ell }}{\partial R}\right]\nonumber \\ & & +\frac{1}{R}\frac{\partial }{\partial R}\left[\left(\nu _{2}R^{2}\left|\frac{\partial \boldsymbol{\ell }}{\partial R}\right|^{2}-\frac{3}{2}\nu _{1}\right)\mathbf {\mathbf {L}}\right]\nonumber \\ & & +\frac{1}{R}\frac{\partial }{\partial R}\left[\nu _{3}R\mathbf {\mathbf {L}}\times \frac{\partial \boldsymbol{\ell }}{\partial R}\right]+\mathbf {T}_{\mathrm{ext}}+ {\bf T}_{\rm src}, \end{eqnarray} \tag{ 11 }$$

where $L=\Sigma \sqrt{GM_{\bullet }R}$, $\boldsymbol{\ell }=\mathbf {L}/L$, and Σ is the disk's mass surface density. The disk is thus effectively modeled as a set of concentric rigid thin annuli. The first term on the right-hand side of Equation (11) describes the angular momentum that is carried to the central sink with the inflowing mass, the second term the angular momentum that is advected outward by the disk's internal viscous torques, which are expressed by the azimuthal kinematic viscosity ν₁ (responsible for the mass inflow) and the vertical kinematic viscosity ν₂ (responsible for the unwarping of the disk). The third term describes the precessional torque that accompanies large amplitude warps (Ogilvie 1999). We assume isotropic viscosities, $\nu _{n}=\widetilde{\alpha }_{n}c_{i}H$ (n = 1, 2, 3), where c_i is the mid-plane isothermal sound speed and use the second-order expansion of the three dimensionless viscosity parameters

$$\begin{equation} \widetilde{\alpha }_{n}=A_{n}+B_{n}\psi ^{2}+O(\psi ^{4})\quad (n=1,2,3) \end{equation} \tag{ 12 }$$

in terms of the local dimensionless local warp amplitude $\psi =R\left|\partial \boldsymbol{\ell }/\partial R\right|$ (Ogilvie 1999; Lodato & Price 2010). The first-order coefficients are

$$\begin{equation} A_{1}=\alpha _{1},\,\,\, A_{2}=\frac{2\big(1+7\alpha _{1}^{2}\big)}{\alpha _{1}\big(4+\alpha _{1}^{2}\big)}\,,\,\,\, A_{3}=\frac{3\big(1\hbox{--}2\alpha _{1}^{2}\big)}{2\big(4+\alpha _{1}^{2}\big)}. \end{equation} \tag{ 13 }$$

We find that the second-order terms do not change the results significantly in our simulations. The second-order coefficients are listed for completeness in Appendix B. For typical values α₁ < 1, ν₂ ≫ ν₁, that is, the disk unwarps faster than it flows into the MBH.

$\mathbf {T}_{\mathrm{ext}}=\mathbf {T}_{\rm BP}+\mathbf {T}_{{\rm RR}}$ expresses the external torques per unit area applied to the disk, which here are those due to relativistic frame dragging and due to the stellar torques. $\mathbf {T}_{\mathrm{SRC}}$ expresses the torque per unit area due to a source term at the disk's outer edge (Section 2.3.1). The effect of frame dragging on the disk, the Bardeen–Petterson (BP) effect, is given to lowest post-Newtonian order by¹ (Bardeen & Petterson 1975)

$$\begin{equation} \mathbf {T}_{\mathrm{BP}}=\boldsymbol{\Omega }_{\mathrm{LT}}\times \mathbf {L}, \end{equation} \tag{ 14 }$$

where $\boldsymbol{\Omega }_{\mathrm{LT}}=(2G/c^{2}R^{3})\mathbf {J}_{\bullet }$ is the Lense–Thirring precession angular frequency and $\mathbf {J}_{\bullet }=\boldsymbol{\chi }GM_{\bullet }^{2}/c$ is the spin of the MBH, with $0\le |\boldsymbol{\chi }|\le 1$ the dimensionless spin parameter. The spin of the MBH evolves correspondingly by (King et al. 2005)

$$\begin{equation} d\mathbf {J}_{\bullet }/dt=-2\pi \int _{R_{\alpha }}^{R_{\beta }}\mathbf {T}_{{\rm BP}}RdR. \end{equation} \tag{ 15 }$$

Since the torque is perpendicular to the MBH spin, it changes only its direction, not its magnitude. The BP torques work to align $\mathbf {L}$ and $\mathbf {J}_{\bullet }$. The frame-dragging torques dominate over the viscous torques out to the distance where the warp diffusion rate becomes faster than the Lense–Thirring precession rate, $R_{\mathrm{BP}}=\sqrt{\nu_{2}(R_{\mathrm{BP}})/\Omega _{\mathrm{LT}}(R_{\mathrm{BP}})}$. Typically R_BP ≫ R_α. The integrand in Equation (15) scales as $\Sigma (R)R^{^{-3/2}}$, so the torque on the MBH is dominated by contributions from the innermost region of the disk that remains non-aligned, and can be estimated by $\dot{J}_{\rm BP}\sim \pi R_{\rm BP}^{2}T_{\rm BP}(R_{\rm BP})=\pi \nu _{2}(R_{\rm BP})\Sigma (R_{\rm BP})\sqrt{GM_{\bullet }R_{\rm BP}}$. The timescale for the MBH spin to align itself with the disk is then $t_{\parallel }\sim J_{\bullet }/\dot{J}_{\rm BP}$, which in steady-state ($\dot{M}\simeq 3\pi \nu _{1}\Sigma$, ν₂ ≃ ν₁/2α²₁) can also expressed as

$$\begin{equation} t_{\parallel }=6\chi \alpha _{1}^{2}(M_{\bullet }/\dot{M})\sqrt{r_{g}/R_{\rm BP}}\,. \end{equation} \tag{ 16 }$$

The torque on the MBH spin by direct accretion of matter, $\dot{J}_{{\rm acc}}\sim \dot{M}\sqrt{GM_{\bullet }R_{\alpha }}$, can be neglected relative to that by the BP effect, $\dot{J}_{\mathrm{BP}}$, since $\dot{J}_{{\rm acc}}/\dot{J}_{\mathrm{BP}}\sim 6\alpha _{1}^{2}\sqrt{R_{\alpha }/R_{\mathrm{BP}}}\ll 1$. Accretion can no longer be ignored on timescales long enough for it to change the MBH mass appreciably, t ≳ t_E = ηcσ_T/4πGm_p = η4.5 × 10⁸ yr, where t_E is the e-folding, or Salpeter, timescale for growth by Eddington-limited accretion, σ_T is the Thomson cross-section, m_p is the proton mass, and η is the radiative efficiency of the accretion.

The internal structure of the disk needs to be specified to determine the viscosity. In the α-disk prescription, the mid-plane temperature T reflects the azimuthal kinetic viscosity of the disk via ν₁ = α₁c²_i/Ω_K, where c²_i = kT/μ is the isothermal sound speed, μ is the mean molecular weight, and $\Omega _{K}=\sqrt{GM_{\bullet }/R^{3}}$ is the Keplerian angular frequency. The sound speed, in turn, determines the disk's exponential scale height H = c_i/Ω_K.

The physical parameters of the masing region in the accretion disk of NGC 4258 can be inferred from the observed aspect ratio H/R ∼ 0.002, which suggests a temperature of T ∼ 600 K (Argon et al. 2007) for a gas pressure supported disk (observations disfavor magnetic support; see review by Lo 2005), and from the fact that H₂O maser emission is possible there. This requires the mid-plane molecular hydrogen density $n_{\mathrm{H}_{2}}$ to be in the range 10⁷ cm⁻³ to 10¹¹ cm⁻³ (all the hydrogen is assumed to be molecular), the mid-plane gas temperature T to be in the range 300–400 K to 1000 K, and the mid-plane pressure p/k to be in the range 10¹⁰ K cm⁻³ to 10¹³ K cm⁻³ (Maloney 2002). In addition, stability against fragmentation requires that Toomre's criterion holds locally,² Q = c_iΩ_K/πGΣ > 1 (Toomre 1964), which further constrains the disk mass and temperature.

It is unclear whether the accretion disk's internal viscosity alone can provide enough heat to warm the molecular gas to masing temperatures. Neufeld & Maloney (1995) argue that heating by X-ray irradiation by the central source is essential for creating the required conditions in the masing region. Outside the masing region, the disk's properties are not usefully constrained by current observations.

Here, we make a choice of convenience to adopt the known solution³ of a gas pressure-dominated thin accretion disk, where the temperature is completely determined by the internal viscosity (that is, no external heating) and free–free (Kramer's law) opacity (Shakura & Sunyaev 1973). The prefactor κ_a of Kramer's law, κ = κ_a(ρ/ρ_a)(T/T_a)^−7/2, where $\rho =\Sigma /(\sqrt{2\pi }H)$ is the total mid-plane mass density (the subscript a denotes values at R_a) is then adjusted to be high enough to maintain the required temperature at the masing region.⁴ We assume solar abundances (X = 0.7057; Arnett 1996) and that all the hydrogen is molecular, which corresponds to a mean molecular weight of μ = 2.358m_p. The adopted value [H/R]_a = 0.002 fixes the surface density $\Sigma {}_{a}=\rho _{a}\sqrt{2\pi }[H/R]_{a}R_{a}$, where $\rho _{a}=n_{\mathrm{H}_{2}}(R_{a})2m_{p}/X$. We find that in order to have masing conditions between R_a and R_b, it is necessary to assume the highest temperature value at the inner edge, T_a = 10³ K. The assumed density there, $n_{H_{2}}(R_{a})=3\times 10^{8}\,\mathrm{cm^{-3}}$, which is within the masing range and consistent with the limits on p_a/k = (ρ_a/μ)T_a, is chosen so the total disk mass is 3 × 10³ M_☉. The mid-plane temperature is given by σ_SBT⁴ = (27/32)Σ²Ω²_Kν₁κ (e.g., Frank et al. 2002, Equation (5.76)), which here yields

$$\begin{equation} T(R,t)=\left[\frac{27k_{B}^{1/2}}{\sqrt{2\pi }32\sigma _{{\rm SB}}}\frac{\alpha _{1}\kappa _{a}}{\mu {}^{1/2}}\frac{T_{a}^{7/2}}{\rho _{a}} \Sigma ^{3}(R,t)\Omega _{K}^{2}(R)\right]^{1/7}\!. \end{equation} \tag{ 17 }$$

For the values of T_a and ρ_a above, Equation (17) requires for self-consistency κ_a = 7.195 × 10⁵/α₁ cm² g⁻¹. This is an extremely high near/mid IR opacity (λ_max ∼ 3–10 μm for T ∼ 300–1000 K), as compared, for example, to the mean theoretical opacity of 〈κ(2.2 μm)〉 = 3800 ± 700 cm² g⁻¹ used to model dusty star-forming cores (Shirley et al. 2011). This is in line with the conclusion that an external source of heating is required to explain the gas temperature.

The finite difference method is used to integrate the evolution equations (Equation (11) with G = c = M_• = 1) for the angular momentum surface densities {L^j_i}_{i = 0}^{N + 1} at time t_j, at positions {R_i}^{N + 1}_{i = 0} on a logarithmically spaced grid between R₁ = R_α at the innermost stable circular orbit (ISCO) and R_N = R_β ≫ R₁. The edge points R₀ and R_{N + 1} are used to enforce the boundary conditions. A brief description of the numerical scheme (Papaloizou & Pringle 1983; Pringle 1992) is given in Appendix A.

2.3.1. Initial Conditions and Boundary Conditions

The initial conditions are a flat disk with surface density Σ = Σ_a(R/R_a)^−3/4 and normal $\boldsymbol{\ell }_{0}$, where Σ is chosen so as to satisfy the maser conditions between the masing region limits R_a and R_b. The MBH spin is initially along the z-axis. We assume the boundary conditions L(R₁) = 0 and $\partial \boldsymbol{\ell }/\partial R|_{R_{1}}=0$ (central mass sink and no warp), and $\partial (\nu _{1}\mathbf {L})/\partial R|_{R_{2}}=0$ (no torque). These boundary conditions allow mass through the inner boundary into the MBH. To prevent the disk from draining on the viscous timescale and from drifting away from the masing conditions much before that, a source term is added at the outer edge. It is adjusted every time-step t_j to compensate for the mass lost from the disk, (ΔM)^j = ∑^N_{i = 1}(M_i^j − M^{j − 1}_i), where M^j_i = 2πR_iΔR_iΣ^j_i is the mass in a disk ring of width ΔR_i = R_i − R_{i − 1} ≪ R_i. The angular momentum density at the outermost grid point is then augmented by an amount $(\Delta L_{N}){}^{j}\boldsymbol{\ell }_{0}$, such that $\left|\mathbf {L}_{N}^{j}+(\Delta L_{N})^{j}\boldsymbol{\ell }_{0}\right|=\left[\Sigma _{N}^{j}+(\Delta \Sigma _{N})^{j}\right]\sqrt{R_{N}}$, where (ΔΣ_N)^j = −X_M(ΔM)^j/2πR_NΔR_N (X_M is an order unity stabilization factor; see below). The magnitude of the angular momentum needed to compensate for the mass loss is therefore

$$\begin{eqnarray} (\Delta L_{N})^{j} & = & -\mathbf {L}_{N}^{j}\cdot \boldsymbol{\ell }_{0}\nonumber \\ & & +\sqrt{\big(\mathbf {L}_{N}^{j}\cdot \boldsymbol{\ell }_{0}\big)^{2}+\big[\Sigma _{N}^{j}+(\Delta \Sigma _{N})^{j}\big]^{2}R_{N}-\big(L_{N}^{j}\big)^{2}}\!.\nonumber \\ \end{eqnarray} \tag{ 18 }$$

The fluctuating factor X_M = 1 ± stabilizes the disk against the accumulation of numeric integration errors by allowing small overcorrections or undercorrections, as needed.⁵ When (ΔM)^j < 0, X_M is set to 1 + , and conversely, when (ΔM)^j > 0, it is set to 1 − . Experimentation indicates that = 0.1 is a good choice.

We find that with the source term thus defined, a flat disk with the inner boundary conditions $\nu _{1}\Sigma |_{R_{1}}=0$ rapidly converges to a close approximation of its theoretically expected steady state solution with the expected mass-loss rate of $\dot{M}\simeq 3\pi \nu _{1}\Sigma$ (for R ≫ R_α). Likewise, the mass of a non-stationary disk rapidly converges to a steady state value that is typically within O(10⁻³) of its initial mass, even as its geometry continuously changes. The mean accretion rate over some period $t_{j_{1}}$ to t_j2 can then be estimated by $\langle \dot{M}\rangle =(\sum _{j=j_{1}}^{j_{2}}X_{M}(\Delta M)^{j})/(t_{j2}-t_{j1})$.

This mass replenishment scheme represents an idealized case of a self-regulating mass supply, which prevents both the draining and overloading and fragmentation of the disk. It should be emphasized that this assumption is introduced here for convenience only, to stabilize the disk long enough to collect robust statistics. Our statistical conclusions about the dynamical response of the disk to the external stellar torques do not depend strongly on the existence of such a mass source, as long as the disk mass is in quasi-equilibrium for at least the relatively short RR timescale.

2.3.2. The Stellar Resonant Relaxation Torques

The RR torques due to the O(108) stars in the radius of influence cannot be modeled directly. Instead, they are approximated here by representing the stellar cusp by a small number of concentric spherical shells delimited by the radii {rk}Nsk = 0 (where r0 = 0). The spacing between consecutive shell radii is chosen to be large enough so that the residual forces are approximately independent of each other, rk + 1/rk ⩾ 22/(3 − γ), which corresponds to the requirement that $\sqrt{N_{\star }({<}\,r_{k})}\ge 2\sqrt{N_{\star }({<}\,r_{k-1})}$ (BA09) and is broadly consistent with the rigorous derivation of Kocsis & Tremaine (2011). The maximal number of such logarithmically spaced independent star shells in the simulations of NGC 4258 that are presented below is typically Ns = 5. The vector RR torques from each shell are represented by the gravitational field of a thin ring of mass $M_{k}=\sqrt{N_{\star }(<r_{k})-N_{\star }(<r_{k-1})}M_{\star }$ with radius $\bar{r}_{k}=(r_{k}+r_{k-1})/2$ and normal $\mathbf {n}_{k}(t)$, which is interpolated smoothly in time by cubic spline from a sequence of isotropic random normals $\lbrace \mathbf {n}_{k}^{j}\rbrace _{j=0}^{[t_{\mathrm{sim}}/t_{0}]+1}$ at times tj = jt0(rk), where tsim is the duration of the simulation and t0(r) = min (tsq(r), treact(r)) (Equations (2) and (9)) is the coherence time. This simulates the random reorientation of the local residual force after a vector RR timescale. Finally, the torque exerted by a star ring at $\bar{r}_{k}$ on a disk gas ring at Ri, which is directed along the line of intersection of the two rings, is calculated numerically by integrating over the gravitational force between all mass elements in the rings (Nayakshin 2005),

$$\begin{eqnarray} \mathbf {T}_{\mathrm{RR}}(k & \rightarrow & i)=\frac{M_{i}M_{k}R_{i}\bar{r}_{k}}{\big(R_{i}^{2}+\bar{r}_{k}^{2}\big)^{3/2}}\frac{\sin \beta }{4\pi ^{2}} \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} & &\times \int _{0}^{2\pi }d\phi _{1}\int _{0}^{2\pi }d\phi _{2}\frac{\sin \phi _{1}\sin \phi _{2}}{[1-\delta \cos \lambda ]^{3/2}}(\boldsymbol{\ell }_{i}\times \mathbf {n}_{k})\,, \end{eqnarray} \tag{ 20 }$$

where M_i is the mass of the disk annulus, $\beta =\cos ^{-1}(\boldsymbol{\ell }_{i}\cdot \mathbf {n}_{k})$, $\delta =1-(\bar{r}_{k}^{2}-R_{i}^{2})/(\bar{r}_{k}^{2}+R_{i}^{2}),$ and⁶ cos λ = cos βsin ϕ₁sin ϕ₂ + cos ϕ₁cos ϕ₂.

2.3.3. Validation of the Code

Figure 2 shows that RR-induced warping can easily reproduce the observed 8° warp angle across the maser region of NGC 4258 (Herrnstein et al. 2005), without requiring any particular initial conditions. The three scenarios explored here, RR with maximal (χ = 1) frame dragging and a large (i₀ = 3/4 rad) initial tilt, or with a small (i₀ = 0 rad) initial tilt, or without frame dragging (χ = 0), all exhibit a probability P₈ ∼ O(1) to have ω ⩾ 8° (see probabilities quoted in Figure 2). This should be contrasted with the scenario where only frame dragging is operating (Caproni et al. 2007). In the example shown here, the warp angle never exceeds 5°. By careful choice of the initial parameters, it may be possible to obtain a large enough warp angle and satisfy the constraints set by the direction of the radio jet with frame dragging only. However, this requires considerable fine tuning and may require a very long-lived disk for the warp to grow (Martin 2008).

Zoom In Zoom Out Reset image size

The warping probability increases with the steepness of the stellar cusp, as more stars are available near the disk to torque it. For example, P₈ = 0.02, 0.11, and 0.29 for i₀ = 0 rad χ = 1 in a sequence of models with γ = 3/2, 7/4, and 2, respectively.

We find that for the i₀ = 0 rad models of NGC 4258, the rms RR-induced misalignment between the MBH spin axis and the disk's total angular momentum is rms(i) = 15° after 10⁸ yr for the case χ = 1, and as large as rms(i) = 44° for the χ = 0 case. This is an example of the stabilizing effect of the BP torques, which introduce a preferred plane for the disk and suppress warps. It is interesting to note that the observed maser disk in NGC 4258 is tilted by ∼ 30° to the present jet axis, as identified by the north and south hotpots (Cecil et al. 2000; Wilson et al. 2001). Assuming that the jet is aligned with the MBH spin axis, then such a tilt is consistent with RR torquing of a disk (see Figure 1(d)), possibly around a sub-maximal spinning MBH.

The warping of the disk exposes it to ionizing radiation from its innermost parts, which may play a central role in determining the physical conditions in the outer regions of the disk. The fraction of the central luminosity that falls on the disk (assuming isotropic emission) is expressed by the disk's covering fraction, C_f = ∫^2π₀dϕ|cos θ_max(ϕ) − cos θ_min(ϕ)|/4π, where θ_min and θ_max are the minimal and maximal inclinations above and below the disk's mean plane along R in azimuthal direction ϕ. Figure 3 shows that RR-induced warping gives the disk a minimal varying covering factor of C_f ∼ 0.1–0.3, on top of any contribution from a preexisting large-scale warp (e.g., by the BP effect).

Zoom In Zoom Out Reset image size

As noted by Lodato & Pringle (2006), warping increases the mass accretion rate over that in a flat disk because of the additional dissipation due to the vertical shear. We confirm here that the warps are associated with mass accretion fluctuations of up to factors of 3–4 over that in a flat disk (Figure 2), resulting in an overall increase by a factor of up to $f_{\dot{M}}\sim 1.4$ in the average mass accretion rate. The power spectrum of the accretion rate fluctuations in the maser model is very broad (Figure 4) and can be roughly approximated by a broken power law, which flattens beyond ≳ 10⁷ yr. This timescale corresponds to the RR coherence time in the outer half of the disk, r ≲ R_d = 0.16 pc to r ∼ R₂ = 0.27 pc. A rough estimate of this variability timescale can be obtained by evaluating the coherence time at the disk's mass-weighted mean radius, since that is where the gravitational coupling with the stars is maximal. That length scale is typically in the regime where the back-reaction time is shorter than the self-quenching time (Equation (9)), and so

$$\begin{equation} t_{{\rm var}}\sim t_{\mathrm{react}}(R_{d})\sim (M_{\bullet }/M_{d})P(R_{d})\simeq 1.1\times 10^{7}\,. \end{equation} \tag{ 21 }$$

A similar correspondence is seen between the variability timescale of the AGN model (Section 3.3), t_var ≃ 1240, and the power-law break in its power spectrum (Figure 4). Thus, t_var can be interpreted as the shortest timescale for which there is substantial variability power.

Zoom In Zoom Out Reset image size

The enhanced mass accretion rate translates directly to a faster rate of change in the magnitude of the MBH spin, $\dot{J}_{\bullet }\simeq \dot{M}\sqrt{GM_{\bullet }R_{\alpha }}$, as the MBH spin vector remains aligned with the inner regions of the disk due to the BP torques. The BP torques also couple the fluctuating RR torques on the disk to the MBH spin orientation (Figure 5, Equation (15)). For the effect to be significant, the disk must be massive enough to carry a substantial amount of angular momentum (Equations (14) and (15)), but not so massive as to remain effectively immobile against the RR torques.

Zoom In Zoom Out Reset image size

We find that the response of the MBH spin in NGC 4258 to the BP torques is very small because of the very low mass of the disk. However, the effect can be larger in systems with a more massive accretion disk. Figure 5 shows as an example the evolution of the MBH spin axis relative to its initial direction (we omit the much smaller spin evolution due to accretion), and that of the total angular momentum of the disk, for a low-mass luminous AGN model that corresponds to a Seyfert galaxy (L ∼ 7 × 10⁴³ erg s⁻¹ for a radiative efficiency of η = 0.1). The angular momenta of the MBH and the disk both execute a random walk around their initial orientation. The amplitude of the disk angular momentum jitter is larger than that of the MBH spin because the RR coherence time, t_react = 620 yr, at R_d = 2.2 × 10⁻³ pc is much shorter than the BP timescale t_∥ = 4.1 × 10⁵ yr, and so the MBH spin cannot "catch up" with the faster changes in the direction of the stellar torques. We find that for this specific model, the MBH spin scatters across ∼ 5° over an e-folding time (4.5 × 10⁷ years). Another consequence of the short coherence time (compared to t_warp(R_d)  ≃  3.5 × 10⁴ yr) is that the RR-induced warps in this disk model have angles of only ∼ 1°. Their effect on the mass accretion rate is correspondingly small, only a ∼ 2% increase over that in a flat disk.

A circumnuclear accretion disk does not exist in isolation around an MBH. Rather, it shares the volume with the high-density nuclear cluster that is predicted to form there while the nucleus evolves. The cluster is expected to be on average spatially isotropic near the MBH, but the Poisson fluctuations in the distribution of orbital inclinations lead to a slowly varying residual force that exerts coherent torques on the disk (RR). We argue (Section 1.3) that a thin Keplerian disk is primarily torqued by these purely gravitational interactions with the stars, rather than by the hydrodynamical ones that occur when stars plunge through the disk. The RR torques warp the disk and can lead to order unity variability in the disk geometry (Figures 1 and 2), its mass accretion rate (Figures 2 and 4), and covering factor (Figure 3). The strong coupling between the stellar potential fluctuations and the disk via RR is then further extended to the MBH spin via the frame-dragging BP effect, and this allows angular momentum to be transferred from the nuclear cluster to the MBH. The combined effects of the perpendicular RR and BP torques excite a jitter in the MBH spin direction (Figure 5). In addition, the increased mass accretion rate due to the RR-induced warping, together with the disk/MBH spin alignment due to the BP effect, leads to an increased growth rate of the MBH spin magnitude. We conclude that gravitational interactions between the stars and the disk excite a substantial level of irreducible variability in the disk properties—stationary accretion in a circumnuclear disk is merely an idealization.

The large-scale warping of the disk (for NGC 4258, this coincides with the maser region warp) reflects the RR coherence timescale at the mass-weighted mean disk radius, t_var ∼ (M_•/M_d)P(R_d) yr (∼10⁷ yr for NGC 4258). Generally, the effect of RR on the disk and the MBH will be substantial when the disk is long-lived, t_disk > t_var. There is a large spread in the estimates of AGN and QSO lifespans in the luminous phase, few ×10⁶–10⁸ yr (e.g., Grazian et al. 2004; Porciani et al. 2004), with some recent analyses suggesting an even longer lifespan of ≲  10⁹ yr (Gilli et al. 2009; Ross et al. 2009; Schawinski et al. 2009). For long-lived disk systems, the cumulative effects of RR-torquing can be large.

One implication of the factor $f_{\dot{M}}\sim 1.2\hbox{--}1.4$ increase in the mean mass accretion rate over that in a flat stationary disk (Section 3.2) is that an MBH fed by an RR-torqued disk, whose mass accretion rate is raised to the Eddington limit by warping, will be $\exp [(f_{\dot{M}}-1)(t_{{\rm disk}}/t_{E})]$ more massive than one that is fed by an otherwise identical flat disk, which is accreting at only $1/f_{\dot{M}}$ of the Eddington limit. For example, RR-torquing can accelerate MBH growth over a time t_disk = 10t_E = 4.5 × 10⁸ yr (assuming η = 0.1) by a factor of up to ∼ 50 (for $f_{\dot{M}}=1.4$). Alternatively, if radiation pressure prevents the increased RR-induced mass flow from accreting on the MBH, the accumulating excess mass may trigger outflows or disk fragmentation and star formation.

Another cumulative effect of RR-torquing of the disk is the displacement of the MBH spin from its initial orientation due to the BP-induced random jitter. It is interesting to note that the few degrees amplitude of the displacement is still consistent with the typical observed ∼ 5° opening angle of AGN radio jets on large scales (Oppenheimer & Biretta 1994), assuming the jet direction reflects the spin direction.

A direct geometric consequence of RR-induced warping is the substantial time-averaged covering factor the disk acquires relative to the central continuum source. This allows the disk to intercept the central radiation and to be heated and ionized by it. Such X-ray irradiation may in fact be essential for raising the gas temperature to the range necessary for maser emission (Neufeld & Maloney 1995), for the production of the observed emission lines in AGN disks in general (Collin-Souffrin 1987), for driving winds and outflows from disks (Proga & Kallman 2004), or for explaining narrow line Seyfert 1 galaxies and ultrasoft AGNs (Puchnarewicz & Soria 2002).

The magnitude of the effects of the stars on the disk depends on the central concentration of the stellar cluster and is therefore sensitive to the degree of mass segregation in the system and to the fraction of massive compact remnants in the population. The more centrally concentrated the stellar density profile, and the more massive the stars (for RR, the relevant quantity is 〈M²_⋆〉/〈M_⋆〉; Rauch & Tremaine 1996), the larger are the effects on the disk. It is worth noting that strong mass segregation with a power-law cusp profiles of γ ≳ 2 is expected to occur in nuclei with old stellar populations (Alexander & Hopman 2009; Keshet et al. 2009; Preto & Amaro-Seoane 2010), where the inner parts of the cusp are dominated by stellar black holes of mass O(10 M_☉). It is thus plausible that a substantial fraction of galactic nuclei have conditions that are conducive to RR-torquing of an accretion disk.

In this study we presented numerical results for galactic nuclei with MBH masses in the range ∼4 × 10⁶ M_☉ (with a massive disk close to $\dot{M}_{E}$) to ∼ 4 × 10⁷ M_☉ (with a low-mass, low-accretion rate disk). In order to scale the results to other systems, it is necessary to specify the scaling of the disk properties with MBH mass. We defer this to future work and present here only a simple preliminary analysis that suggests that the effects are generally more important in lower-mass MBHs. Assuming that all AGN disks extend up to some universal large multiple of r_g (e.g., ∼ 2000r_g for α-disks limited by gravitational instability; Goodman 2003) and assuming they all have the same aspect ratio, then R_d∝M_• and M_d/M_• ≃ H/R, so that M_d∝M_•. The angular momentum in the disk therefore scales as $J_{d}\propto M_{d}\sqrt{M_{\bullet }R_{d}}\propto M_{\bullet }^{2}$. The M_•/σ relation indicates that the MBH radius of influence scales as r_h∝M^1/2_• (assuming M_•∝σ⁴), so that the number of stars in the MBH cusp inside R_d scales as N_⋆(R_d)∝M_•(R_d/r_h)^{3 − γ}∝M^{(5 − γ)/2}_•. The residual angular momentum in the stars that are efficiently coupled to the disk then scales as $J_{N}\propto \sqrt{N_{\star }(R_{d})}\sqrt{M_{\bullet }R_{d}}\propto M_{\bullet }^{(9-\gamma)/4}$, and it then follows that J_d/J_N∝M^{(γ − 1)/4}_•. The ratio of the disk and stellar angular momenta is a slowly rising function of M_• for γ > 1. This implies that, all other parameters being equal, it should become progressively more difficult for the stars to torque the disk as the MBH mass increases. In addition, as argued in Section (1.3), the larger M_•, the more the disk responds to the RR torques as a rigid body, rather than by growing warps.

The calculations and conclusions presented here are limited by various assumptions and approximations. They are listed here briefly. Some of these issues will be addressed in future work.

Limitations of the galactic nucleus model. The RR torques by the stellar cluster are approximated by the gravitational field of a small number of thin rings (Section 2.3.2), which provides only a rough approximation of the true power spectrum of the stellar perturbations. The back-reaction of the disk on the stars is approximated simplistically by limiting the RR coherence time so it does not exceed the back-reaction time. A more realistic fluctuation spectrum will allow, among other issues, to better model the very short timescale fluctuations and explore whether these have any relation to the observed optical/UV variability of AGNs. Our conclusions on mass accretion rates and MBH spin evolution are generalizations based on extrapolating simulations that were constructed explicitly to represent the disk and cluster of NGC 4258. They should be scaled and tested for a wider range of MBHs, disks, and nuclear clusters.

Limitations of the disk model. The physics underlying accretion disks are still not well understood. In particular the nature of the azimuthal and vertical viscosities, and the relation between them are quite uncertain (Ogilvie 1999; Lodato & Pringle 2007; Lodato & Price 2010). The disk is not modeled self-consistently. The thermal structure of the disk is based for simplicity on a Kramer opacity law with a free normalization constant, which is clearly non-physical, as indicated by the very high opacity that is required to produce the masing conditions. External X-ray heating, which is likely important, is neglected. MBH spin evolution due to mass accretion can be relevant on longer timescales, but is not taken into account, and the artificial BP torque suppression on small scales that is needed to stabilize the results numerically may both lead to an underestimation of the effect of disk warping on the MBH spin evolution. Finally, the disk is approximated as Newtonian and Keplerian, even near the ISCO, which is held fixed at its Schwarzschild value.

Limitations of the numerical scheme. The numerical scheme used here (Papaloizou & Pringle 1983; Pringle 1992) cannot describe azimuthal modes and as implemented is limited to moderate warp angles (although viscosities for arbitrary large angles can be calculated numerically; Ogilvie 1999). In particular, disk flipping from co-rotation to counterrotation, or disk destruction by large deformations, cannot be modeled reliably.

We analyzed and simulated numerically the evolution of a thin accretion disk around a MBH that is surrounded by a stellar cluster. We took into account the disk's internal viscous torques, the frame-dragging torques of a spinning MBH, and the stellar orbit-averaged gravitational torques. We show that the evolution of the MBH mass accretion rate, the MBH spin growth rate, and the covering fraction of the disk relative to the central ionizing continuum source, are all strongly coupled to the stochastic fluctuations of the stellar potential via the warps that the stellar torques excite in the disk. These lead to fluctuations by factors of up to a few in these quantities over a wide range of timescales, with most of the power on timescales ≳  (M_•/M_d)P(R_d). The response of the disk is stronger the lighter it is and the more centrally concentrated the stellar cusp. We demonstrated these effects by simulating the evolution of the maser disk in NGC 4258 and show that its observed O(10°) warp can be driven by the stellar torques. We also show that the frame dragging of a massive AGN disk couples the stochastic stellar torques to the MBH spin and can excite a jitter of a few degrees in its direction relative to that of the disk's outer regions.