Local Stability and Global Instability in Iron-opaque Disks

The energy output observed in Galactic X-ray binaries (XRBs), and active galactic nuclei (AGNs), suggests that the root of emitted power in these sources must be connected with the gravitational potential energy of a compact object. In most of the former, and all of the latter, this object is a black hole. Its mass can range from a few solar masses up to a few billion of M_⊙. The material that is accreted onto a black hole and emits radiation may possess substantial angular momentum. In this case, the accretion flow forms a geometrically thin disk, which is located in the equatorial plane of an XRB, or co-aligned with the plane perpendicular to the black hole rotation axis in AGN, if the latter is powered by a Kerr black hole. The classical solution of Shakura & Sunyaev (1973) with their α-viscosity prescription describes a stationary accretion disk where the dissipated heat is balanced by thermal radiation. As studied by Pringle et al. (1973), Lightman & Eardley (1974), and Shakura & Sunyaev (1976), the viscous stress tensor scaling with a total (i.e., gas plus radiation) pressure leads to the runaway instability of the disk structure. Alternatively, the viscosity prescription may be given by a gas pressure only (in this case the thermal instability does not develop), or by some intermediate law that invokes a combination of gas and radiation with different weights (see, e.g., Szuszkiewicz 1990). For instance, a general prescription that was recently discussed by Grzȩdzielski et al. (2017) reads as

$$\begin{eqnarray}&&{\tau }_{r\phi }=\alpha {P}_{\mathrm{tot}}^{\mu }{P}_{\mathrm{gas}}^{1-\mu }.\end{eqnarray} \tag{ 1 }$$

The above model leads to the unstable disk behavior, which manifests in a limit-cycle type of oscillation of the emitted luminosity, characteristic for the periodically heated and cooled inner regions of the accretion disk. This kind of cycle is possible if only the thermal runaway is captured by some stabilizing process. This might be due to the advection of heat onto a black hole, as proposed for the so-called "slim disk" solution (Abramowicz et al. 1988). The presence of radiation-pressure instability in the action of cosmic sources has been a matter of debate (see, e.g., a review by Blaes 2014). Nevertheless, there are strong observational hints that support the limit-cycle type of behavior in at least two well-known microquasars, GRS 1915+105 and IGR J17091-324, in some of their spectral states (Belloni et al. 2000; Altamirano et al. 2011). The limit-cycle oscillations were also detected in the ultraluminous X-ray source HLX-1, which is claimed to contain an intermediate-mass black hole (Farrell et al. 2009, 2012; Lasota et al. 2011; Servillat et al. 2011; Wu et al. 2016). Furthermore, a type of nonlinear dynamics characteristic for an underlying unstable accretion flow was suggested for a number of other XRBs (Sukova et al. 2016), while the statistical studies of a large sample of sources support the "reactivation" scenario in the case of compact radio sources hosting supermassive black holes (Czerny et al. 2009). On the other hand, many of XRBs and AGN seem to be powered by a stable accretion, despite even large accretion rates. Therefore, some stabilizing mechanisms in the accretion process have been invoked, apart from the viscosity prescription itself. For instance, the propagating fluctuations in the flow (Janiuk & Misra 2012) may suppress the thermal instability, or at least delay the development of the instability (Ross et al. 2017). As recently discussed by Jiang et al. (2016), a possible stabilizing mechanism for the part of an accretion flow might be the changes in opacity that are connected with the ionization of heavy elements. Using the shearing-box simulations of MRI-driven fluid in the gravitational potential of a supermassive black hole (a particular value of black hole mass, M = 5 × 10⁸ M_⊙, was used), Jiang et al. (2016) have shown that the flow is stable against the thermal instability, if the opacity includes transitions connected with absorption and scattering on iron ions. This is because the cooling rate, which includes now not only the Thomson scattering (constant) term, but also the absorption and line emission in the Roseland mean opacity, will depend strongly on density and temperature in specific regions of the disk. In fact, as we discuss below in more detail that the dominant term from the opacity changes may completely stabilize the flow locally. The simulations of Jiang et al. (2016) showed that effect, but they did not describe the global evolution of the flow, which is the subject of our present work. In our paper we consider the disk stability for the range of black hole masses, characteristic for either XRBs or AGN. We find that the influence of the opacity changes on the global time evolution of the flow is essential, although it does not prevent the instability from developing. To illustrate this, we show examples of light curves produced by our numerical simulation. The paper is organized as follows. In Section 2 we present the analytical condition for local thermal instability and in Section 3 we present the values of the κ opacities for the typical accretion disk densities, temperatures, and range of the iron opacity bump to reveal its influence on global disk behavior in Section 4.

The domination of radiation pressure in the accretion disk leads to the thermal instability (Pringle et al. 1973; Lightman & Eardley 1974; Shakura & Sunyaev 1976; Janiuk et al. 2002). At the instability, the heating rate Q₊ grows faster with temperature than cooling rate Q₋. The appearance of local thermal instability is given by the condition

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{log}{Q}_{-}}{d\mathrm{log}T}\lt \displaystyle \frac{d\mathrm{log}{Q}_{+}}{d\mathrm{log}T}.\end{eqnarray} \tag{ 2 }$$

The analysis performed in Grzȩdzielski et al. (2017), under the assumption that constant surface density during thermal timescales leads to the following formula on heating rate derivative:

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{log}{Q}_{+}}{d\mathrm{log}T}=1+7\mu \displaystyle \frac{1-\beta }{1+\beta },\end{eqnarray} \tag{ 3 }$$

where β = P_gas/P. For the case in this work, we assume α = 0.02 and μ = 0.56 to be typical for the IMBH and AGN disk case. We have chosen these values since the dynamics of the outbursts of the disk match those of the observed properties of the sources as shown in Grzȩdzielski et al. (2017). These values reproduce the correlation between the bolometric luminosity and the outburst duration that is known for the observed sources, especially for microquasars and intermediate-mass black holes. The radiative cooling rate depends on the disk surface density, Σ, the physical constants of Stefan–Boltzmann σ_b, and the speed of light. The radiative cooling rate is given by the following formula (Janiuk et al. 2002, 2015; Grzȩdzielski et al. 2017):

$$\begin{eqnarray}&&{Q}_{-}=\displaystyle \frac{4{\sigma }_{b}}{3c{\rm{\Sigma }}}\displaystyle \frac{{T}^{4}}{\kappa }.\end{eqnarray} \tag{ 4 }$$

For the radiation-pressure dominated disk, combined with Equations (2) and (3), from Equation (4) for μ = 0.56 we get

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{log}{Q}_{-}}{d\mathrm{log}T}\lt 4.92.\end{eqnarray} \tag{ 5 }$$

In our previous papers, only Thomson scattering was assumed, which resulted in the appearance of global radiation-pressure instability among all scales of the sub-Eddington accretion disks. The recent results of Jiang et al. (2016) are based on the stabilizing influence of the iron opacity components. To confront the results of their short-time, local 3D MHD shearing-box simulation, we propose the global model, which was previously used for the case of the intermediate-mass BH HLX-1 accretion disk (Wu et al. 2016; Grzȩdzielski et al. 2017). According to the lower temperatures in radiation-pressure dominated areas of the accretion disks, we also include in our model the atomic opacity components. Assuming κ to be a function of ρ and T, we get the following formula for log derivative of Q₋:

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{log}{Q}_{-}}{d\mathrm{log}T}=4+\displaystyle \frac{\partial \mathrm{log}\kappa }{\partial \mathrm{log}T}-\displaystyle \frac{4-3\beta }{1+\beta }\displaystyle \frac{\partial \mathrm{log}\kappa }{\partial \mathrm{log}\rho }.\end{eqnarray} \tag{ 6 }$$

Equation (6) shows the possible stabilizing influence of the negative slope of κ that is dependent on T, and the destabilizing influence of the positive slope of κ. Also, the dependence on ρ has an influence on the disk stability. The value of κ itself is not important in the local stability analysis. In the case of an unstable disk, the effect of less efficient cooling of the greater κ lowers the temperature of the unstable equilibrium solution, and enlarges the temperature of stable solution, which can modify the duty cycle quantitatively, but not qualitatively. Similarly to Grzȩdzielski et al. (2017), we can derive the necessary value for the thermal instability as

$$\begin{eqnarray}&&\beta \lt \displaystyle \frac{7\mu +3-\tfrac{\partial \mathrm{log}\kappa }{\partial \mathrm{log}T}+4\tfrac{\partial \mathrm{log}\kappa }{\partial \mathrm{log}\rho }}{7\mu -3+\tfrac{\partial \mathrm{log}\kappa }{\partial \mathrm{log}T}-4\tfrac{\partial \mathrm{log}\kappa }{\partial \mathrm{log}\rho }}.\end{eqnarray} \tag{ 7 }$$

We can also define the thermal stability parameter s as

$$\begin{eqnarray}&&s=\displaystyle \frac{d\mathrm{log}{Q}_{+}}{d\mathrm{log}T}-\displaystyle \frac{d\mathrm{log}{Q}_{-}}{d\mathrm{log}T}.\end{eqnarray} \tag{ 8 }$$

Using the Equations (3) and (6), we can write Equation (8) as

$$\begin{eqnarray}&&s=-3+7\mu \displaystyle \frac{1-\beta }{1+\beta }-\displaystyle \frac{\partial \mathrm{log}\kappa }{\partial \mathrm{log}T}+\displaystyle \frac{4-3\beta }{1+\beta }\displaystyle \frac{\partial \mathrm{log}\kappa }{\partial \mathrm{log}\rho }.\end{eqnarray} \tag{ 9 }$$

The s parameter is connected with the Lyapunov exponent for the system described by the energy equation in the accretion disk with the stress tensor given by Equation (1). The value s > 0 means that the disk is locally thermally unstable, and the values ≤ 0 means the the disk is locally thermally stable. Assuming μ = 0.56 and β ≪ 1 for ρ = 10⁻⁸ g cm⁻³ and T = 10⁵ K, we get s = −11 which corresponds to the local thermal stability. Nevertheless, for this density the parameter s gains positive values for T > 3 × 10⁵ K. Below this temperature, for T > 1.75 × 10⁵ K, the disk is locally thermally stable because of the negative slope of the bump. For temperatures in the range of (1.1–1.75) × 10⁵ K, the disk is locally thermally unstable.

The opacity κ is the local function describing the interaction of photons with matter from accretion disks. Under the assumptions of the local thermal equilibrium, the radiation and gas contribution to the total pressure, and the local vertical hydrostatic equilibrium, both the heating and cooling rates can be described as a function of radius r, local density ρ, and local temperature T. Although for the radius r, affecting the angular momentum transport is important for the heating rate in the α-disk model, it affects the stability analysis only indirectly via the parameters of the stationary solutions. In Figure 1 we present the profiles of the total opacity for solar metallicity, computed as a function of density and temperature (Alexander et al. 1983; Seaton et al. 1994; Różańska et al. 1999). The combined conditions (3) and (4) to the opacity values result in the significant local stabilization for the temperatures of (1–4) × 10⁵ K and densities of about 10⁻⁸ g cm⁻³ that are typical for the AGN accretion disks. We fitted the κ function with the following formulae:

$$\begin{eqnarray}\kappa & = & {\kappa }_{\mathrm{Th}}+{\kappa }_{{pl}}+{\kappa }_{\mathrm{bump}},{\kappa }_{\mathrm{Th}}=0.34\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1},\\ {\kappa }_{{pl}} & = & 4.6\times {10}^{23}\rho {T}^{-3.5},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray*}{\kappa }_{\mathrm{bump}} & = & 39.8{\rho }^{0.2}\left(0.8\exp \,-\,{\left(\displaystyle \frac{T-1.75\times {10}^{5}{\rm{K}}}{8.2\times {10}^{4}{\rm{K}}}\right)}^{2}\right.\\ & & +\left.6.3\exp \,-\,{\left(\displaystyle \frac{T-4\times {10}^{4}{\rm{K}}}{3\times {10}^{4}{\rm{K}}}\right)}^{2}\right).\end{eqnarray*}$$

Zoom In Zoom Out Reset image size

The negative stabilizing slope of the iron opacity bump is visible in Figure 1. In Figure 2, the upper panel shows the opacities directly from the tables (Alexander et al. 1983; Seaton et al. 1994; Różańska et al. 1999), and in the lower panel we present the analytical approximation of the opacity function from Equation (3). The detailed results of the dynamical model are presented in Section 4.

Zoom In Zoom Out Reset image size

4.1. Values of $\rho$ and T

In Figure 3 we present typical values of ρ and T for the wide range of accretion disks, computed via the GLADIS code. For the values of ρ and T in the upper left corner, the power-law term of κ dominates, but matter with this parameter is too dense and too cold for the central areas of sub-Eddington accretion disks. The oblique belt below presents typical values of ρ and T for the accretion disks with different masses. This area, for $5.1\lt \mathrm{log}T\lt 5.4$, is covered by the iron opacity bump with a stabilizing negative slope (see Section 2). According to these results, the stabilizing effect of the negative slope of the iron opacity bump can be visible only for the AGN accretion disks with M ≈ 10⁸–10⁹ M_⊙.

Zoom In Zoom Out Reset image size

4.2. Results for the Full Model with Bump

We perform the simulations of the global disk behavior using the time-dependent global code GLADIS (Grzȩdzielski et al. 2017). The GLADIS code is a time-dependent code, which solves hydrodynamic equations, describing the long-time behavior of accretion disks under the assumption of the vertical hydrostatic equilibrium. The code models the time evolution of the flow in thermal and viscous timescales. Moreover, the code assumes axial symmetry and Keplerian angular velocity. In this paper we set α = 0.02, μ = 0.56, and M = 5 × 10⁸ M_⊙. The major change in comparison to Grzȩdzielski et al. (2017) is the replacement of the constant Thomson κ with Equation (3). Similarly to Jiang et al. (2016), we assumed the Eddington rate to be $\dot{m}=0.03$. The results of the time-dependent model are presented in Figure 4. The local shearing-box simulation resulted in the significant stabilization of the disk (Jiang et al. 2016). However, the global model does not confirm these results. The stabilization of the disk, which appears according to local prediction, is not found in global models considering a large range of radii. Figure 4 presents the stability parameter profile s as defined in Equation (8). For the inner area of the disk, the Thomson component of opacity dominates and the temperature is too large to expect any form of stabilization. The outer area of the disk characterizes the larger value of total opacity temperature to be about (1.5–2) × 10⁵ K and the negative values of s (the bump temperatures). The significant gradient of s is correlated with the gradient of temperature and gradient of κ in the opposite direction. The stability parameter, as presented in Figure 4, can reach values between −3 and −3 + 7 μ (0.92 for our choice of μ). In Figure 4 the typical profile of the s parameter is presented. As the bump is an approximately Gaussian function centered at 1.75 × 10⁵ K with a standard deviation of σ = 0.82 × 10⁵ K, it is expected that the strongest effect of the stabilizing slope would be visible for such temperatures. The combined outcome of the stabilizing influence of the negative slope of the bump and destabilizing influence of the positive slope of the bump for the dynamical model is presented in Figure 5. In contrast to the results of Jiang et al. (2016), the iron bump does not stabilize the disk. However, it complicates the light-curve pattern (many small short flares preceding main outbursts instead of one simple flare) due to the complexity of the photon absorption process, but the inner regions of the disk remain hot enough to perform the limit-cycle oscillations. In effect, the bump partially stabilizes the disk—the amplitude L_max/L_min decreases from 156 for the model with constant κ (lower panel of Figure 5) to 16.1 for the model with the bump (upper panel of Figure 5). The detailed analysis of the light-curve shape is presented in Table 1. Since the timescales presented in Figure 5 are much longer than that of the duration of the observation (lasting up to several decades), it is impossible to find such a light curve using the direct method. For such a black hole mass, during the phase of fastest growth of the luminosity, the luminosity change can reach 1% per year. The shape of the light curve can be reflected in the Eddington rate statistic where similar objects, being in the different phase of the limit-cycle presented in Figure 4, can emit the radiation with a different luminosity and spectra. However, in cases of much smaller black hole masses, new timescales are perhaps accessible to the observations. For example, digitalization of Harvard plates (Grindlay et al. 2012)¹ will bring light curves on the order of a hundred years and perhaps the outburst of AGN disks can be discovered.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Flare Parameters for the Model from Figure 5 with κ Described by Formula (3) and Thomson κ

Model	Amplitude A	Period P (years)	Width Δ
κ with bump	16.1	12635	0.0023
Thomson κ	156	70197	0.0044

Note. Those parameters were described more precisely in Grzȩdzielski et al. (2017).

Download table as:  ASCIITypeset image

In our previous paper (Grzȩdzielski et al. 2017) we computed large grids of models confirming the universality of radiation-pressure instability across the BH mass-scale. In this paper we changed the opacity prescription to examine the influence of heavy atoms on the accretion disk instability. The comparison between the two models presented in Table 1 leads to the conclusion that heavy atoms stabilize the disk partially, but do not imply that the variability vanishes. This stabilizing effect manifests itself rather well in a significant period and amplitude decrease, without a relative broadening of the outbursts with respect to their separation. Additionally, some mild precursors, being an outcome of a non-monotonic profile of the s parameter distribution, are also visible. That partial stabilization, being an important effect for the AGN, has only weak influence on the radiation pressure among all of the BH mass-scales in accretion disks under the assumption of solar metallicity. In the case of sources with different metallicity, this effect can change its extent. Finally, we conclude that the radiation-pressure driven limit-cycle oscillations, suffering some disturbances from the Iron Opacity Bump in the case of the AGN disks, are also expected, at least for moderately large supermassive black holes (M = 5 × 10⁸ M_⊙).

This work was supported in part by the grants DEC-2012/05/E/ST9/03914 and 2015/18/M/ST9/00541 from the Polish National Science Center.

Software: GLADIS, (Janiuk et al. 2002), (Janiuk et al. 2015).