Accretion in Radiative Equipartition (AiRE) Disks

Since the early years after Shakura & Sunyaev (1973) introduced thin disks, they have been known (Shakura & Sunyaev 1976; Piran 1978) to be fraught with thermal instabilities. The instability occurs in the inner regions of the disk, where the pressure is dominated by radiation pressure and the opacity is mostly due to Thomson scattering, resulting in a much stronger temperature dependence in the heating of the disk compared to the cooling of the disk.

Decades later, the resolution of thermal instabilities still remains one of the major outstanding problems in understanding thin and slim disks. One of the major uncertainties in the Shakura & Sunyaev (SS) α disk model is the assumption that the viscous stress is proportional to the total pressure. Early attempts to model thermal (and viscous) stability, such as the works of Sakimoto & Coroniti (1981), Stella & Rosner (1984), and Merloni (2003), explored the possibility that the viscous stress might instead be proportional only to the gas pressure. These disks were called β disks, where ${t}_{r\phi }=\beta \,{p}_{\mathrm{gas}}$. Recent numerical simulations (Brandenburg et al. 1995; Stone et al. 1996; Jiang et al. 2013; Mishra et al. 2016; Sa̧dowski 2016a, 2016b) see the presence of thermal instabilities, where the onset of thermal instability causes the disk to expand or collapse at the timescale of only a few orbits. These local simulations do not see evidence for such β disks.

Radiative MHD simulations such as those in Sa̧dowski (2016a, 2016b) find stable radiatively efficient and strong radiation-pressure-dominated disks in the presence of strong magnetic fields. This has led to the claim that strong magnetic fields could stabilize disks against thermal instabilities (Begelman & Pringle 2007; Oda et al. 2009; Sa̧dowski 2016a, 2016b). If there is not enough magnetic flux, however, instability once again sets in and a different means for stabilization must be sought. An iron opacity bump has also been suggested by Jiang et al. (2016) as a means to postpone (but not avoid) thermal instabilities. In this paper, we present accretion in radiative equipartition (AiRE) disks as an alternative solution to the thermal instability problem in thin and slim disks.

In slim (and thin) accretion disks, the main approximations made are that the disk is axisymmetric, stationary, and $h\lt r$. The background metric is assumed to be Kerr.

Let us define the following expressions involving the black hole spin:

$$\begin{eqnarray}{\rm{\Delta }} & = & {M}^{4}({r}_{* }^{2}-2{r}_{* }+{a}_{* }^{2}),\\ A & = & {M}^{4}({r}_{* }^{4}+{r}_{* }^{2}{a}_{* }^{2}+2{r}_{* }{a}_{* }^{2}),\\ { \mathcal C } & = & 1-3{r}_{* }^{-1}+2{a}_{* }{r}_{* }^{-3/2},\\ { \mathcal D } & = & 1-2{r}_{* }^{-1}+2{a}_{* }^{2}{r}_{* }^{-2},\\ { \mathcal H } & = & 1-4{a}_{* }{r}_{* }^{-3/2}+3{a}_{* }^{2}{r}_{* }^{-2},\end{eqnarray} \tag{ 1 }$$

where ${a}_{* }=a/M$ and ${r}_{* }=r/M$. When quantities have been vertically integrated, a polytropic equation of state has been assumed: $p={ \mathcal K }\ {\rho }^{1+1/N}$ where N = 3 and ${ \mathcal K }=\mathrm{const}$. In the notation that follows, ${\rm{\Sigma }},{T}_{c},$ and P are vertically integrated quantities (whereas p, for example, is not). Unless otherwise noted, we use the geometrized units, $G=c=1$. $\dot{M}$ is the mass accretion rate, and we use the convention that ${\dot{M}}_{\mathrm{Edd}}=16{L}_{\mathrm{Edd}}/{c}^{2}$.³

Following the notation of Sa̧dowski (2011) and Abramowicz & Fragile (2013), the relativistic equations describing slim disks are:

• 1.  
  Conservation of (rest) mass:

  $$\begin{eqnarray}&&\dot{M}=-2\pi {\rm{\Sigma }}\sqrt{{\rm{\Delta }}}\displaystyle \frac{v}{\sqrt{1-{v}^{2}}},\end{eqnarray} \tag{ 2 }$$

  where v is the radial velocity and Σ is the surface density.
• 2.  
  Conservation of radial momentum:

  $$\begin{eqnarray}&&\displaystyle \frac{v}{1-{v}^{2}}\displaystyle \frac{{dv}}{{dr}}=\displaystyle \frac{{ \mathcal A }}{r}-\displaystyle \frac{1}{{\rm{\Sigma }}}\displaystyle \frac{{dP}}{{dr}},\end{eqnarray} \tag{ 3 }$$

  where

  $$\begin{eqnarray}&&{ \mathcal A }=-\displaystyle \frac{{MA}}{{r}^{3}{\rm{\Delta }}{{\rm{\Omega }}}_{k}^{+}{{\rm{\Omega }}}_{k}^{-}}\displaystyle \frac{({\rm{\Omega }}-{{\rm{\Omega }}}_{k}^{+})({\rm{\Omega }}-{{\rm{\Omega }}}_{k}^{-})}{1-{\tilde{{\rm{\Omega }}}}^{2}{\tilde{r}}^{2}},\end{eqnarray} \tag{ 4 }$$

  $\tilde{r}=A/({r}^{2}\sqrt{{\rm{\Delta }}})$, and Ω is the angular velocity with respect to the stationary observer, while $\tilde{{\rm{\Omega }}}={\rm{\Omega }}-\tfrac{2{Mar}}{A}$ is the angular velocity with respect to the inertial observer. ${{\rm{\Omega }}}_{k}^{\pm }=\pm \sqrt{M}/({r}^{3/2}\pm a\sqrt{M})$ are the angular velocities of corotating and counter-rotating Keplerian (or circular geodesic) orbits.
• 3.  
  Conservation of angular momentum:

  $$\begin{eqnarray}&&\displaystyle \frac{\dot{M}}{2\pi }({ \mathcal L }-{{ \mathcal L }}_{\mathrm{in}})=\displaystyle \frac{\sqrt{A{\rm{\Delta }}}\gamma }{r}\alpha P,\end{eqnarray} \tag{ 5 }$$

  where ${ \mathcal L }$ is the angular momentum per unit mass, ${{ \mathcal L }}_{\mathrm{in}}$ is an integration constant (approaching the angular momentum at the innermost stable circular orbit (ISCO) in the thin disk limit), and

  $$\begin{eqnarray}&&\gamma =\sqrt{\displaystyle \frac{1}{1-{v}^{2}}+\displaystyle \frac{{{ \mathcal L }}^{2}{r}^{2}}{A}},\end{eqnarray} \tag{ 6 }$$

  is the Lorentz factor. We adopt the α viscosity prescription, where the ${t}_{r\phi }$ component (in the comoving frame of the fluid) of the viscous stress tensor can be expressed as ${t}_{r\phi }=-\alpha \,p$ (Shakura & Sunyaev 1973). The relation between Ω and ${ \mathcal L }$ is

  $$\begin{eqnarray}&&{\rm{\Omega }}=\displaystyle \frac{2{aMr}}{A}+\displaystyle \frac{{r}^{3}{ \mathcal L }\sqrt{{\rm{\Delta }}}}{{A}^{3/2}\gamma }.\end{eqnarray} \tag{ 7 }$$
• 4.  
  Vertical equilibrium:

  $$\begin{eqnarray}&&{h}^{2}{{\rm{\Omega }}}_{\perp }^{2}=9\displaystyle \frac{P}{{\rm{\Sigma }}},\end{eqnarray} \tag{ 8 }$$

  where ${{\rm{\Omega }}}_{\perp }=\sqrt{\tfrac{M}{{r}^{3}}\tfrac{{ \mathcal H }}{{ \mathcal C }}}$ is the vertical epicyclic frequency and h is the half disk thickness.
• 5.  
  Conservation of energy:

  $$\begin{eqnarray}&&{{ \mathcal Q }}^{\mathrm{adv}}={{ \mathcal Q }}^{\mathrm{vis}}-{{ \mathcal Q }}^{\mathrm{rad}}.\end{eqnarray} \tag{ 9 }$$

  Here, the advective cooling is defined as

  $$\begin{eqnarray}&&\,{{ \mathcal Q }}^{\mathrm{adv}}=\displaystyle \frac{1}{r}\displaystyle \frac{d}{{dr}}[{{ru}}^{r}(E+P)]-{u}^{r}\displaystyle \frac{{dP}}{{dr}}-{\int }_{-h}^{h}{u}^{z}\displaystyle \frac{\partial p}{\partial z}{dz},\end{eqnarray} \tag{ 10 }$$

  where E is the vertically integrated energy density $E={\int }_{-h}^{h}\left(\tfrac{3}{2}{p}_{\mathrm{gas}}+3{p}_{\mathrm{rad}}\right)\ {dz}$. Performing the z integral in ${Q}^{\mathrm{adv}}$, we get

  $$\begin{eqnarray}{{ \mathcal Q }}^{\mathrm{adv}} & = & \displaystyle \frac{\dot{M}}{2\pi {r}^{2}}\left({\eta }_{3}\displaystyle \frac{P}{{\rm{\Sigma }}}\displaystyle \frac{d\mathrm{ln}P}{d\mathrm{ln}r}-(1+{\eta }_{3})\displaystyle \frac{P}{{\rm{\Sigma }}}\displaystyle \frac{d\mathrm{ln}{\rm{\Sigma }}}{d\mathrm{ln}r}\right.\\ & & +\left.{\eta }_{3}\displaystyle \frac{P}{{\rm{\Sigma }}}\displaystyle \frac{d\mathrm{ln}{\eta }_{3}}{d\mathrm{ln}r}+{{\rm{\Omega }}}_{\perp }^{2}{\eta }_{4}\displaystyle \frac{d\mathrm{ln}{\eta }_{4}}{d\mathrm{ln}r}\right),\end{eqnarray} \tag{ 11 }$$

  where

  $$\begin{eqnarray}&&{\eta }_{1}=\displaystyle \frac{128}{315}\,h,\end{eqnarray} \tag{ 12 }$$

  $$\begin{eqnarray}&&{\eta }_{2}=\displaystyle \frac{8}{9},\end{eqnarray} \tag{ 13 }$$

  $$\begin{eqnarray}&&\,{\eta }_{3}=\displaystyle \frac{1}{P}\left(\displaystyle \frac{1}{5/3-1}\displaystyle \frac{k}{\mu {m}_{p}}\displaystyle \frac{8}{9}{\rm{\Sigma }}\,{T}_{c}+\displaystyle \frac{256}{315}{a}_{\mathrm{rad}}\,{T}_{c}^{4}\,h\right),\end{eqnarray} \tag{ 14 }$$

  $$\begin{eqnarray}&&{\eta }_{4}=\displaystyle \frac{1}{18}\,{h}^{2},\end{eqnarray} \tag{ 15 }$$

  and where k is the Boltzmann constant, m_p is the mass of the proton, and $\mu =0.62$ is the solar abundance.The viscous heating ${{ \mathcal Q }}^{\mathrm{vis}}$ is fixed by the α-prescription

  $$\begin{eqnarray}&&{{ \mathcal Q }}^{\mathrm{vis}}=-\alpha P\displaystyle \frac{A{\gamma }^{2}}{{r}^{3}}\displaystyle \frac{d{\rm{\Omega }}}{{dr}},\end{eqnarray} \tag{ 16 }$$

  while the radiative cooling ${{ \mathcal Q }}^{\mathrm{rad}}$ is

  $$\begin{eqnarray}&&{{ \mathcal Q }}^{\mathrm{rad}}=\displaystyle \frac{64\sigma {T}_{c}^{4}}{3{\rm{\Sigma }}\kappa },\end{eqnarray} \tag{ 17 }$$

  where T_c is the temperature at the equatorial plane, and the opacity coefficient κ is given by Kramer's formula (Sa̧dowski 2011)

  $$\begin{eqnarray}&&\kappa ={\kappa }_{\mathrm{es}}+{\kappa }_{{ff}}=0.34+6.4\times {10}^{22}\rho {T}^{-7/2}\end{eqnarray} \tag{ 18 }$$

  in cgs units, and using the solar abundance $\rho \sim \tfrac{{\rm{\Sigma }}}{2h}$ and $T\sim {T}_{c}$. The factor $\tfrac{64}{3}$ in front of the radiation term in (17) is somewhat arbitrary, as it depends on the details of the assumed vertical structure and radiative transfer, and thus its exact value should be taken with a grain of salt.⁴
• 6.  
  The vertically integrated equation of state:

  $$\begin{eqnarray}&&P={\eta }_{2}\displaystyle \frac{k}{\mu {m}_{p}}{\rm{\Sigma }}\,{T}_{c}+\displaystyle \frac{2}{3}\,{\eta }_{1}\,{a}_{\mathrm{rad}}\,{T}_{c}^{4},\end{eqnarray} \tag{ 19 }$$

  where the first term is the gas pressure ${P}_{\mathrm{gas}}={\eta }_{2}\tfrac{k}{\mu {m}_{p}}{\rm{\Sigma }}\,{T}_{c}$, and the second term is the radiation pressure ${P}_{\mathrm{rad}}\,=\tfrac{2}{3}\,{\eta }_{1}\,{a}_{\mathrm{rad}}\,{T}_{c}^{4}$.

In the thin disk (SS) limit of the above equations, Keplerian orbits are assumed (${\rm{\Omega }}\to {{\rm{\Omega }}}_{k}$), while the disk becomes radiatively efficient (${{ \mathcal Q }}^{\mathrm{adv}}\to 0$). The system then simplifies to an algebraic one where closed-form solutions may be found in three regimes: (1) an outer region where ${P}_{\mathrm{gas}}$ and ${\kappa }_{{ff}}$ dominate, (2) a middle region where ${P}_{\mathrm{gas}}$ and ${\kappa }_{\mathrm{es}}$ dominate, and (3) an inner region where ${P}_{\mathrm{rad}}$ and ${\kappa }_{\mathrm{es}}$ dominate.

The condition for thermal instability can be written as (Pringle 1976)

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial }{\partial h}({{ \mathcal Q }}^{\mathrm{vis}}-{{ \mathcal Q }}^{\mathrm{rad}})\right|}_{{\rm{\Sigma }}}\gt 0.\end{eqnarray} \tag{ 20 }$$

For stellar mass black holes, the slim disks in the previous section are thermally unstable above $\dot{M}/{\dot{M}}_{\mathrm{Edd}}\sim 0.06$ in their radiation-pressure-dominated regimes. For supermassive black holes, this limit is lower (see Figure 1). The thermal instability sets in because in the innermost regions of the disk we have $p\sim {p}_{\mathrm{rad}}$ and $\kappa \sim$ const., leading to

$$\begin{eqnarray}&&{{ \mathcal Q }}^{\mathrm{vis}}=-\displaystyle \frac{3}{2}{\rm{\Omega }}\,{t}_{r\phi }h\propto p\,h\propto {h}^{2},\end{eqnarray} \tag{ 21 }$$

where we used the vertical equilibrium of $p\propto {{\rm{\Omega }}}^{2}\,{\rm{\Sigma }}\,h$ and assumed that Σ and α are constant, while

$$\begin{eqnarray}&&{{ \mathcal Q }}^{\mathrm{rad}}=16\displaystyle \frac{c\,{p}_{\mathrm{rad}}}{\kappa {\rm{\Sigma }}}\propto h.\end{eqnarray} \tag{ 22 }$$

Zoom In Zoom Out Reset image size

Since $h\propto p\propto {T}_{c}^{4}$, a temperature increase would change the rate of heating to be much faster than that of cooling, making the disk thermally unstable.

There is also a relation between the shape of the ${T}_{c}({\rm{\Sigma }})$ curve and the stability of the disk. The equilibrium states of the disk can be described by the ${T}_{c}({\rm{\Sigma }})$ solutions, at a fixed radius, to ${{ \mathcal Q }}^{\mathrm{adv}}=0$. This ${T}_{c}({\rm{\Sigma }})$ relation is sometimes referred to as an S-curve (Abramowicz et al. 1988; Lasota 2016) because of the shape it takes in the models and the temperature ranges that are often considered. Points on this curve with a positive slope correspond to stable equilibria, while points on this curve with a negative slope correspond to unstable equilibria (Lasota 2016). In the former case small perturbations to the temperature bring the state back to equilibrium, while in the latter case these perturbations lead to runaway heating or cooling.

Figure 1 shows the radii below which radiation pressure, ${P}_{\mathrm{rad}}$, starts to dominate over gas pressure, ${P}_{\mathrm{gas}}$, for black holes of different masses and α's. We will then assume that at this point in the accretion flow, and inward, the onset of thermal instability creates an inhomogeneous disk structure, where photons can effectively escape faster than they would through diffusion in a smooth disk. In other words, we hypothesize that the effective optical depth (or opacity) of the disk drops due to the instability, or we have more efficient cooling. This could be, e.g., due to rising photon bubbles, low density funnels, or other inhomogeneous structures that could invalidate (17).

We assume that this efficient cooling can cool the disk down to pressure equipartition. However, the cooling would not drive the disk to become gas pressure dominated again, as the thermal instability responsible for faster cooling ceases in this regime and the disk heats back up within a viscous time. Therefore, the state of marginal thermal instability, or pressure equipartition, is a stable fixed point (see Section 4.1 for details). Indeed, since disks at high Eddington ratios are observed in a steady state, we conjecture that the cooling is efficient enough to keep the disk at equipartition, i.e., ${P}_{\mathrm{gas}}\approx {P}_{\mathrm{rad}}$.⁵

We will call these flows AiRE disks. This equipartition tames the temperature dependence of the disk in the inner region. Within the radii of Figure 1, AiRE disks have different properties from slim and thin disks. In this model, the equipartition condition replaces the radiative energy loss condition (17), which is only valid for atmospheres in vertical equilibrium with planar symmetry. More generally, many MHD/fluid instabilities occur when there are different fluid components that dominate inertia and (isotropic or anisotropic) stress. Some examples are magnetorotational instability (MRI), convective instability, Rayleigh–Taylor instability, Parker instability, and potentially whatever sets the maximum mass of main-sequence stars. In all of these, the instability saturates at equipartition, where stress and inertia are (somewhat) equally distributed in dominant components. For example, in gas-dominated disks, MRI saturates when ${P}_{\mathrm{mag}}={ \mathcal O }(0.1)\times {P}_{\mathrm{gas}}$ (Hawley et al. 1995).

4.1. Efficient Cooling and Pressure Equipartition

Let us now consider how efficient cooling can lead to pressure equipartition.

Since the development of inhomogeneities that (could) lead to more efficient cooling is driven by thermal instability, we can write a phenomenological toy model for cooling rate as

$$\begin{eqnarray}{\tilde{{ \mathcal Q }}}^{\mathrm{rad}} & \equiv & {{ \mathcal Q }}^{\mathrm{rad}}{\left(1+\mathrm{Re}[\sqrt{2-5\beta }]\displaystyle \frac{{t}_{\mathrm{vis}}}{{t}_{\mathrm{th}}}\right)}^{n}\\ & = & \displaystyle \frac{64\sigma {T}_{c}^{4}}{3{\rm{\Sigma }}\kappa }{\left(1+\mathrm{Re}[\sqrt{2-5\beta }]\displaystyle \frac{{t}_{\mathrm{vis}}}{{t}_{\mathrm{th}}}\right)}^{n},\end{eqnarray} \tag{ 23 }$$

where $\beta =\tfrac{{P}_{\mathrm{gas}}}{{P}_{\mathrm{gas}}+{P}_{\mathrm{rad}}}$, t_th, and t_vis are the thermal and viscous timescales, respectively, and $n\gt 0$ is a free parameter. n = 0 corresponds to the standard disk model from Equation (17), for which ${\tilde{{ \mathcal Q }}}^{\mathrm{rad}}\equiv {{ \mathcal Q }}^{\mathrm{rad}}$. As the rate of thermal instability is ${\lambda }_{\mathrm{th}}\sim \mathrm{Re}[\sqrt{2-5\beta }]{t}_{\mathrm{th}}^{-1}$, which operates over an accretion/viscous time of $\sim {t}_{\mathrm{vis}}$, we expect the cooling to be roughly enhanced by ${({\lambda }_{\mathrm{th}}{t}_{\mathrm{vis}})}^{n}$ with $n={ \mathcal O }(1)$. Therefore, Equation (23) gives a reasonable toy model for the instability-driven cooling.

Figure 2 shows the thermal equilibria solutions to ${\tilde{{ \mathcal Q }}}^{\mathrm{adv}}\equiv {{ \mathcal Q }}^{\mathrm{vis}}-{\tilde{{ \mathcal Q }}}^{\mathrm{rad}}=0$ at different fixed radii ($r=15\,M,20\,M,30\,M,40\,M$ from top left to bottom right) for n = 0, 0.05, 0.1, and 1. The dashed curve, L, represents solutions with conserved angular momentum. The intersection of this L curve with each equilibrium curve denotes an equilibrium point that we are interested in. The intersection with the standard disk model is marked with a circle and the intersection with the AiRE disk model is marked with a star. At larger radii these intersections are closer to one other, but as we move inward in the disk, the AiRE disk gives us different equilibria compared to that of the standard model. The slope at the circle is negative, corresponding to an unstable equilibrium, while the slope at the star is positive, corresponding to a stable equilibrium. This confirms our expectation that efficient cooling (23), driven by thermal instability, leads to radiative pressure equipartition, ${P}_{\mathrm{gas}}\approx {P}_{\mathrm{rad}}$, which is thermally stable.

Zoom In Zoom Out Reset image size

We think that local shearing box simulations such as Jiang et al. (2013) have not yet seen this pressure equipartition realized, either due to insufficient running time or more likely due to the limited box size. Most of the runs in these simulations show either expanding or collapsing vertical height. However, in a big enough box, both should happen in different places, which could lead to the inhomogeneities we have just described.

Mishra et al. (2016) also see thermal instability in simulations of initially radiation-pressure-dominated thin disks. They argue that a comparison of their evolution with the relevant thin disk thermal equilibrium curve suggests that their disk may be headed for the thermally stable, gas-pressure-dominated branch, meaning that it is moving toward equipartition. This supports our conjecture. Their simulations, however, had to be terminated before equipartition could occur because they reached a point where they could no longer resolve the disk.

4.2. AiRE Disk Equations

Using ${P}_{\mathrm{gas}}\approx {P}_{\mathrm{rad}}$, which for AiRE disks effectively replaces Equation (17) in the slim disk model, the equation of mass conservation (2) and vertical equilibrium (8), T_c and v, can be related:

$$\begin{eqnarray}&&\dot{M}\approx -2\pi \sqrt{{\rm{\Delta }}}v\displaystyle \frac{\sqrt{18k/\mu {m}_{p}}\,\xi \,{{T}_{c}}^{7/2}}{{{\rm{\Omega }}}_{\perp }},\end{eqnarray} \tag{ 24 }$$

where $\xi \equiv \tfrac{2}{3}\times \tfrac{128}{315}\times \tfrac{{a}_{\mathrm{rad}}}{{\eta }_{2}k/\mu {m}_{p}}$.

Thus, with the assumption of equipartition of gas pressure with radiation pressure, the equations of the previous section can be combined into one ordinary differential equation, namely the equation of radial momentum conservation (3). This equation can be rewritten in the form

$$\begin{eqnarray}&&\displaystyle \frac{{dv}}{{dr}}=\displaystyle \frac{{ \mathcal N }(r,v)}{{ \mathcal D }(r,v)},\end{eqnarray} \tag{ 25 }$$

where

$$\begin{eqnarray}&&{ \mathcal N }(r,v)=\displaystyle \frac{{ \mathcal A }}{r}\end{eqnarray} \tag{ 26 }$$

and

$$\begin{eqnarray}&&{ \mathcal D }(r,v)\approx v-1.3\displaystyle \frac{k}{v\,\mu {m}_{p}}{\left(-\displaystyle \frac{{\eta }_{2}\dot{M}{{\rm{\Omega }}}_{\perp }\sqrt{\tfrac{k}{\mu {m}_{p}}}}{{a}_{\mathrm{rad}}v\sqrt{{\rm{\Delta }}}}\right)}^{2/7}.\end{eqnarray} \tag{ 27 }$$

Figures 3 and 4 show the (dimensionless) Toomre parameter and central temperature for AiRE disks in comparison with SS disks, using the Keplerian approximation (${\rm{\Omega }}\to {{\rm{\Omega }}}_{k}$) for different nonspinning black hole masses with $\dot{M}/{\dot{M}}_{\mathrm{Edd}}=0.1$ and $\alpha =0.01$. Work is in progress to find the spectrum of the AiRE disks and to compare it to the SS spectrum (Y. K. Yazdi & N. Afshordi 2017, in preparation). It would be interesting to see if the AiRE disk spectrum may resolve some of the discrepancies that come from using the SS spectrum, such as in the determination of black hole spins (Abramowicz & Fragile 2013).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The thin disk solution in the outer region is used to set the boundary condition far outside the disk. In the innermost region, the flow must continuously make a transition from subsonic to supersonic flow. This occurs at a radius called the sonic point. The denominator, ${ \mathcal D }$, in Equation (25) vanishes at the sonic point, making the equation singular, unless the numerator, ${ \mathcal N }$, also vanishes at this point. Thus, at the sonic point we must have the inner boundary condition of

$$\begin{eqnarray}&&{ \mathcal N }{| {}_{r={r}_{\mathrm{sonic}}}={ \mathcal D }| }_{r={r}_{\mathrm{sonic}}}=0.\end{eqnarray} \tag{ 28 }$$

The value of the angular momentum constant ${{ \mathcal L }}_{\mathrm{in}}$ must be chosen carefully to achieve this. If we call the value of ${{ \mathcal L }}_{\mathrm{in}}$, which achieves this ${{ \mathcal L }}_{\mathrm{in}}^{t}$, then for smaller values of this constant ${{ \mathcal L }}_{\mathrm{in}}\lt {{ \mathcal L }}_{\mathrm{in}}^{t}$, ${ \mathcal N }$ will vanish but  ${ \mathcal D }$ will not (there will not be a transition from subsonic to supersonic flow), while for greater values of this constant ${{ \mathcal L }}_{\mathrm{in}}\gt {{ \mathcal L }}_{\mathrm{in}}^{t}$, ${ \mathcal D }$ will vanish but  ${ \mathcal N }$ will not (there will be a singularity).

The radial momentum in Equation (25) may be solved by the shooting method, the relaxation method, or a combination of the two (Press et al. 2007). The topology of v differs for solutions with ${{ \mathcal L }}_{\mathrm{in}}\gt {{ \mathcal L }}_{\mathrm{in}}^{t}$ and ${{ \mathcal L }}_{\mathrm{in}}\lt {{ \mathcal L }}_{\mathrm{in}}^{t}$, and this topology change is key to finding the solution via the shooting method. Figure 5 shows radial velocity profiles for different values of ${{ \mathcal L }}_{\mathrm{in}}$ for $M=10\,{M}_{\odot }$, $\dot{M}/{\dot{M}}_{\mathrm{Edd}}=0.1$, a = 0, and $\alpha =0.01$. Solutions with a minimum have ${{ \mathcal L }}_{\mathrm{in}}\lt {{ \mathcal L }}_{\mathrm{in}}^{t}$, whereas the solutions that diverge have ${{ \mathcal L }}_{\mathrm{in}}\gt {{ \mathcal L }}_{\mathrm{in}}^{t}$. The solution with the desired inner boundary condition (28) lies at the transition between these two branches of solutions.

Zoom In Zoom Out Reset image size

An additional challenge to this problem is that the location of the sonic point is not known in advance. This can be handled by treating the sonic point as a free boundary (Press et al. 2007) when using the relaxation method.

Once we have found the solution with the correct inner boundary condition, we can study the nature of the sonic point by considering the Jacobian matrix ${ \mathcal J }=\left(c\begin{array}{cc}\tfrac{\partial { \mathcal D }}{\partial r} & \tfrac{\partial { \mathcal D }}{\partial v}\\ \tfrac{\partial { \mathcal N }}{\partial r} & \tfrac{\partial { \mathcal N }}{\partial v}\end{array}\right)$. The relative signs of the eigenvalues of this matrix, at the sonic point, characterize the sonic point. If the eigenvalues have the same sign, the sonic point is a nodal point (Figure 6), if they have opposite signs, it is a saddle point (Figure 7), and if they are complex, the sonic point is a spiral (unphysical) point (Figure 7). Given that perturbations can only propagate downstream beyond the sonic point, we conjecture that saddle-type sonic points correspond to stable disk configurations, while nodal points correspond to unstable ones. To see this, note that small perturbations inside nodal⁶ sonic points grow as they get dragged deeper into the black hole, while small perturbations inside saddle sonic points shrink (see Figures 8–10). Other authors have made similar arguments about the instability of nodal sonic points (e.g., Kato et al. 1988, 1993).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For fixed values of $\dot{M}/{\dot{M}}_{\mathrm{Edd}},\,M$, and a, there is a transition from the saddle to the nodal type of sonic point as we increase α. Figure 11 shows some transition lines for different masses and spins of black holes. The transition occurs at

$$\begin{eqnarray}{\alpha }_{\mathrm{saddle}} & \lesssim & 0.77{\left(-\displaystyle \frac{{v}_{r}}{{v}_{\phi }}\right)}^{1/3}\\ & \simeq & 0.097{(1-a)}^{2/15}{(3+a)}^{1/3}{\left(\displaystyle \frac{\dot{M}/{\dot{M}}_{\mathrm{Edd}}}{M/{M}_{\odot }}\right)}^{1/24},\end{eqnarray} \tag{ 29 }$$

where ${v}_{r}=\tfrac{1}{\gamma }\tfrac{v}{\sqrt{1-{v}^{2}}}$ and ${v}_{\phi }=\tfrac{2{Ma}}{r\sqrt{{\rm{\Delta }}}}+\tfrac{1}{\gamma }\tfrac{{ \mathcal L }r}{\sqrt{A}}$. Based on our conjecture, we think that values of α greater than $\sim 0.77{\left(-\tfrac{{v}_{r}}{{v}_{\phi }}\right)}^{1/3}$ correspond to nodal sonic points and are thus unstable, while AiRE disks with smaller α's have healthy saddle sonic points. A similar expression for this transition was found by Afshordi & Paczynski (2003) in pseudo-Newtonian isothermal disks, but with $\tfrac{h}{r}$ instead of $\tfrac{{v}_{r}}{{v}_{\phi }}$. Our expression is more general in that it holds for arbitrary spin in full general relativity.

Zoom In Zoom Out Reset image size

Assuming that the value of α is fixed by the MHD physics in thin accretion disks, Equation (29) suggests a physical origin for the minimum Eddington ratio observed for the soft states of X-ray binaries (e.g., see Kalemci et al. 2013), which is around 1%–3%. Combining these values with the current measurements of black hole masses and spins (e.g., McClintock et al. 2014) gives a range of $\alpha =0.11\mbox{--}0.13$ for the viscosity parameter.

We should note that our proposed mechanism for the hard to soft transition would only be consistent with observations if α lies in this narrow range. Current variations of α in simulations are driven by their numerical and/or physical limitations (e.g., thermal instability). However, a better understanding of the physics involved could sharpen these predictions. In other words, we have a strong prediction for how α should behave in more realistic accretion simulations, but a very weak prediction for the transition Eddington ratios (in the absence of a precise value of α).

Self-gravity becomes important in disks outside of some radius, r_sg. When self-gravity becomes important, matter clumps together and will no longer accrete onto the compact object. The more massive the compact object, the closer r_sg gets to the ISCO. When ${r}_{\mathrm{sg}}={r}_{\mathrm{ISCO}}$, there will no longer be any accretion of luminous matter. Since accretion is the main mechanism by which supermassive black holes grow, this can be used to set an upper mass limit for supermassive black holes, as suggested in King (2016).

The radius at which self-gravity takes over can be determined by the Toomre parameter (Toomre 1964), $Q=\tfrac{{c}_{s}\,{{\rm{\Omega }}}_{r}}{\pi \,{\rm{\Sigma }}}$, where ${{\rm{\Omega }}}_{r}=\sqrt{\tfrac{-\tfrac{3{a}^{2}}{{r}^{2}}+8a\sqrt{M}{r}^{-3/2}-\tfrac{6M}{r}+1}{{r}^{3}\left(2a\sqrt{M}{r}^{-3/2}-\tfrac{3M}{r}+1\right)}}$ is the radial epicyclic frequency for circular, equatorial Kerr geodesics (Gammie 2004). The condition for stability is that $Q\gt 1$. Figure 3 shows how Q decreases with mass and approaches 1 for masses of $\sim {10}^{11}\,{M}_{\odot }$.

We find that for the disk to remain Toomre stable ($Q\gt 1$ everywhere outside the ISCO), we must have

$$\begin{eqnarray}&&M/{M}_{\odot }\leqslant f(a,\alpha )\sqrt{{\dot{M}}_{\mathrm{Edd}}/\dot{M}},\end{eqnarray} \tag{ 30 }$$

where f(0, 0.1) ≈ 6.3 × 10²¹, f(0.999, 0.1) ≈ 1.1 × 10²², and f(0, 0.01) ≈ 6.3 × 10²⁰. Restricting to sub-Eddington mass accretion rates $\dot{M}/{\dot{M}}_{\mathrm{Edd}}\leqslant 1$, together with the constraint from Toomre stability (30), we get

$$\begin{eqnarray}\displaystyle \frac{{dM}}{{dt}} & \leqslant & \min \left[{\dot{M}}_{\mathrm{Edd}},\displaystyle \frac{{\dot{M}}_{\mathrm{Edd}}({M}_{\odot })f(a,\alpha )}{M}\right]\\ & \approx & {\dot{M}}_{\mathrm{Edd}}({M}_{\odot }){\left(\displaystyle \frac{{m}^{2}}{f{(a,\alpha )}^{2}}+\displaystyle \frac{1}{{m}^{2}}\right)}^{-1/2}.\end{eqnarray} \tag{ 31 }$$

Integrating (31) we have

$$\begin{eqnarray}&&t({M}_{f})\gtrsim {\int }_{{M}_{i}}^{{M}_{f}}{dm}\displaystyle \frac{1}{{\dot{M}}_{\mathrm{Edd}}({M}_{\odot })}\sqrt{\displaystyle \frac{{m}^{2}}{f{(a,\alpha )}^{2}}+\displaystyle \frac{1}{{m}^{2}}},\end{eqnarray} \tag{ 32 }$$

where masses and t are in units of M_⊙. Inverting this relation to get M_f(t), assuming that the first massive black holes were born at $z\sim 30$, and converting t to z to get M(z), we arrive at the upper bounds shown in Figure 12 for nonspinning and nearly maximally spinning supermassive black holes. In the range of redshift shown in Figure 12, the mass bounds are not sensitive to the starting mass of the seed black hole (M_i). Also included in Figure 12 are the upper bounds (for typical values of parameters) given in King (2016), as well as a Swift satellite observation of S50014+813 at z = 3.366, where its mass was found to be $4\times {10}^{10}\,{M}_{\odot }$ (Ghisellini et al. 2009).

Zoom In Zoom Out Reset image size

To summarize our main results, we have introduced AiRE disks as a solution to the thermal instability in thin disks. The key feature of AiRE disks is that the radiation pressure is in equipartition with the gas pressure in the inner region. We have presented some features of these flows such as their central temperature and Toomre parameter profiles. We have derived upper limits for the mass of supermassive black holes due to the gravitational Toomre instability in AiRE disks. We have also found a transition from the saddle to the nodal type of the sonic points in AiRE disks and used nodal point instability to place a lower limit on the mass accretion rate as a function of viscosity parameter α and black hole spin. We conjecture that this transition might be responsible for the observed lower limits on the Eddington ratio of the soft state in X-ray binaries.

While we introduced AiRE disks to provide a thermally stable description of thin accretion flows, they may also significantly refine our understanding of other disk properties. With new observations from the advanced Laser Interferometer Gravitational-wave Observatory and the Event Horizon Telescope, we are at the advent of a new era of black hole physics. Disk theory may play a major role in explaining some of their future findings.

In upcoming work (Y. K. Yazdi & N. Afshordi 2017, in preparation), we study the spectrum of AiRE disks and its properties. Furthermore, the onset of Toomre instability in the inner regions of AiRE disks around active galactic nuclei can lead to formation and merger of binary black hole systems, such as the ones recently detected by LIGO (Abbott et al. 2016), and lead to characteristic detectable electromagnetic signatures (e.g., Bartos et al. 2017).

Another important future direction is the study of the AiRE disk regime in full MHD radiative transfer simulations and whether enhanced cooling leading to pressure equipartition can indeed arise in a realistic setting.

We would like to thank Shane Davis, Jeremy Goodman, Ramesh Narayan, Olek Sa̧dowski, and Jim Stone for useful discussions. This research was supported in part by Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development Canada, and by the Province of Ontario through the Ministry of Research, Innovation, and Science.