The Dynamics of Truncated Black Hole Accretion Disks I: Viscous, Hydrodynamic Case

TL;DR

This paper presents the first detailed hydrodynamic simulation of a truncated black hole accretion disk, revealing how thermal instability and convective motions influence disk dynamics and outflows.

Contribution

It introduces a well-resolved viscous hydrodynamic model of truncated disks, highlighting the role of thermal instability and convective cells in their evolution.

Findings

Thermal instability causes quasi-periodic transitions in disk cooling states.

Outflows are launched from hot/cold gas interfaces, affecting disk structure.

Convective cells driven by thermal instability influence angular momentum transport.

Abstract

Truncated accretion disks are commonly invoked to explain the spectro-temporal variability from accreting black holes in both small systems, i.e. state transitions in galactic black hole binaries (GBHBs), and large systems, i.e. low-luminosity active galactic nuclei (LLAGNs). In the canonical truncated disk model of moderately low accretion rate systems, gas in the inner region of the accretion disk occupies a hot, radiatively inefficient phase, which leads to a geometrically thick disk, while the gas in the outer region occupies a cooler, radiatively efficient phase that resides in the standard geometrically thin disk. Observationally, there is strong empirical evidence to support this phenomenological model, but a detailed understanding of the dynamics of truncated disks is lacking. We present a well-resolved viscous, hydrodynamic simulation that uses an ad hoc cooling prescription to drive a thermal instability and, hence, produce the first sustained truncated accretion disk. With this simulation, we perform a study of the dynamics, angular momentum transport, and energetics of a truncated disk. We find that time variability introduced by the quasi-periodic transition of gas from efficient cooling to inefficient cooling impacts the evolution of the simulated disk. A consequence of the thermal instability is that an outflow is launched from the hot/cold gas interface which drives large, sub-Keplerian convective cells in the disk atmosphere. The convective cells introduce a viscous $\theta-\phi$ stress that is less than the generic $r-\phi$ viscous stress component, but greatly influences the evolution of the disk. In the truncated disk, we find that the bulk of the accreted gas is in the hot phase.

Figures

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{v})$$

$$\frac{\partial }{\partial t}(\rho \mathbf{v})+\nabla \cdot (\rho \mathbf{v}\mathbf{v}+P\mathcal{I})$$

$$\frac{\partial }{\partial t}(E+\rho \Phi )+\nabla \cdot \left[(E+P+\rho \Phi )\mathbf{v}\right]$$

$$E=u+\frac{1}{2}\rho |\mathbf{v}{|}^2,$$

$$\Phi =-\frac{GM}{r}.$$

$$\rho(R,\theta)=\rho_{0}R^{-3/2}\exp\Bigg{(}-\frac{z^{2}}{2c_{s}^{2}R^{3}}\Bigg{)}\Bigg{[}1-\Bigg{(}\frac{R_{*}}{R}\Bigg{)}^{\frac{1}{2}}\Bigg{]}$$

$$\nu =\frac{2}{3}\alpha P{\Omega }_K^{-1}=\frac{2}{3}\alpha \rho T{\Omega }_K^{-1},$$

$$\mathcal{H}=\frac{3}{2}\alpha \mathcal{R}\rho T{\Omega }_K,$$

$$\Lambda ={\Lambda }_0(T)\rho T^{2}f(r),$$

$$\frac{h}{r}=\sqrt{\frac{3\alpha {\mathcal{R}}^2}{2{\Lambda }_0(T)r^2{\Omega }_{K}f(r)}}.$$

$$\frac{h}{r}=\sqrt{\frac{3\alpha }{2{\Lambda }_0(T)}},$$

$${\Lambda }_0(T)=\{\begin{array}{l}{\Lambda }_c(T<T_T)\\ {\Lambda }_h(T\ge T_T)\end{array}$$

$${\left(\frac{h}{r}\right)}_{\mathrm{hot}}=\sqrt{\frac{3\alpha }{2{\Lambda }_h}}\mathrm{for}r<r_{tr,i}:=\frac{3\alpha GM}{{\Lambda }_c{\mathcal{R}}^2{T}_T}\text{vspace}0.01cm$$

$${\left(\frac{h}{r}\right)}_{\mathrm{cold}}=\sqrt{\frac{3\alpha }{2{\Lambda }_c}}\mathrm{for}r>r_{tr,o}:=\frac{3\alpha GM}{{\Lambda }_h{\mathcal{R}}^2{T}_T}$$

$${\left(\frac{h}{r}\right)}_{\mathrm{cold}}\mathrm{or}{\left(\frac{h}{r}\right)}_{\mathrm{hot}}\mathrm{for}{r}_{tr,i}<r<r_{tr,o}.$$

$$\frac{r_{tr,o}}{r_{tr,i}}={\left[\frac{(h/r{)}_{\mathrm{hot}}}{(h/r{)}_{\mathrm{cold}}}\right]}^2$$

$$\left(\frac{h}{r}\right)$$

$$r_{tr,i}$$

$$\frac{h_{g}(r)}{r}=\Bigg{\langle}\sqrt{\frac{\int(\theta(r)-\bar{\theta}(r))^{2}\rho d\Omega}{\int\rho d\Omega}}\Bigg{\rangle},$$

$$\bar{\theta }(r)=\frac{\int \theta (r)\rho d\Omega }{\int \rho d\Omega }$$

$$\frac{h_{T}(r)}{r}=\Bigg{\langle}\frac{\int\sqrt{Tr}\rho d\Omega}{\int\rho d\Omega}\Bigg{\rangle},$$

$$\Sigma =\frac{1}{2\pi r}\int \rho r^{2}d\Omega .$$

$$v_K=\frac{1}{\sqrt{r}}.$$

$$\dot{M}(R)=\int -\rho v_r{r}^{2}d\Omega .$$

$$\displaystyle\frac{\partial}{\partial t}\bigg{(}\rho rv_{\phi}\bigg{)}=\nabla\cdot\bigg{(}\eta r^{2}\frac{\partial}{\partial r}\bigg{(}\frac{v_{\phi}}{r\sin\theta}\bigg{)}\boldsymbol{\hat{\textbf{e}}}_{r}$$

$$\displaystyle+\eta\sin\theta\frac{\partial}{\partial\theta}\bigg{(}\frac{v_{\phi}}{\sin\theta}\bigg{)}\boldsymbol{\hat{\textbf{e}}}_{\theta}-\rho rv_{\phi}\hbox{${\bf v}$}\bigg{)},$$

$$\frac{\partial \mathcal{E}}{\partial t}+\mathcal{C}+{\mathcal{A}}_{\mathcal{K}}+{\mathcal{A}}_{\mathcal{T}}+\mathcal{G}-\mathcal{W}=0,$$

$$\displaystyle\frac{\partial\mathcal{E}}{\partial t}=\frac{\partial}{\partial t}\bigg{[}\int d\Omega\bigg{(}\frac{1}{2}\rho v^{2}+\frac{P}{\gamma-1}+\rho\Phi\bigg{)}\bigg{]},$$

$$\mathcal{C}=\int d\Omega \Lambda ,$$

$$\displaystyle\mathcal{A_{K}}=\frac{1}{r^{2}}\frac{\partial}{\partial r}\bigg{[}r^{2}\int d\Omega\bigg{(}\frac{1}{2}\rho v^{2}\bigg{)}v_{r}\bigg{]},$$

$$\displaystyle\mathcal{A_{T}}=\frac{1}{r^{2}}\frac{\partial}{\partial r}\bigg{[}r^{2}\int d\Omega\bigg{(}\frac{P}{\gamma-1}\bigg{)}v_{r}\bigg{]},$$

$$\displaystyle\mathcal{G}=\frac{1}{r^{2}}\frac{\partial}{\partial r}\bigg{[}r^{2}\int d\Omega\bigg{(}\rho\Phi\bigg{)}v_{r}\bigg{]},$$

$$\displaystyle\mathcal{W}=\frac{1}{r^{2}}\frac{\partial}{\partial r}\bigg{\{}r^{2}\int d\Omega\eta v_{\phi}\bigg{[}r\frac{\partial}{\partial r}\bigg{(}\frac{v_{\phi}}{\sin\theta}\bigg{)}$$

$$\displaystyle+\sin\theta\frac{\partial}{\partial\theta}\bigg{(}\frac{v_{\phi}}{\sin\theta}\bigg{)}\bigg{]}\bigg{\}}.$$

$$Be=\frac{1}{2}{\mathbf{v}}^2+\Phi +W,$$

Full text

1. Introduction

The radiation produced by the accretion of gas onto a compact object is highly variable, both spectrally and temporally. One compelling example of this variable behavior is the spectral state changes from outbursting black hole binaries (BHBs) and the associated hysteresis. BHBs are relatively rare, but the several dozen that have been identified are found to have common evolutionary properties during their state transitions (e.g. Cygnus X-1, 4U 1543-47, H1743-332, GRO J1655-40, GX 339-4).

Much uncertainty surrounds this behavior, but spectral timing studies find that, in a broad sense, there are two primary BHB accretion states, low and high, with distinct properties and additional transitionary intermediate states with composite properties (Remillard & McClintock, 2006). Of the two main states, the low state is less luminous and the BHB X-ray spectrum is dominated by a power law with index $\sim 1.7$. A thermal spectral component is present, but typically contributes around $20\%$ of the $2-20$ keV flux. Radio emission originating from a jet is usually observed and, in some objects with high luminosities, an outflowing wind can be found (Homan et al., 2016), too. In the high state, a BHB is more luminous. Thermal, blackbody emission dominates the X-ray spectrum and accounts for $>75\%$ of the unobscured $2-20$ keV flux. The peak of this thermal emission is around $1$ keV, which is taken as evidence the geometrically thin disk extends close to the inner-most stable circular orbit (ISCO). In this state, a jet is absent, but signatures of outflowing winds are usually found in the X-ray spectrum. Additionally, quasiperiodic oscillations (QPOs) are observed during accretion state transitions, with high-frequency QPOs observed in certain intermediate states and low-frequency QPOs found in some low-states.

While these common behaviors exist, every outbursting BHB is unique and has its own peculiarities. Nevertheless, the gross differences between the low and high state properties seem to naturally be explained with a model where the accretion occurs through either a hot or cold mode. These different modes are associated with different accretion flow structures and varying degrees of inner disk truncation (for review see Zdziarski & Gierliński, 2004; McClintock & Remillard, 2006; Done et al., 2007). The preferred mode depends on the cooling efficiency of the inner region of the accretion flow as the accretion rate evolves; specifically, the flow is in a radiatively inefficient mode in the low state and a radiatively efficient mode in the high state. In the low state, the inefficient cooling in the inner region effectively truncates the cold, thin disk since the thermal energy of the gas cannot be radiated away. The inner disk is replaced by a hot, diffuse flow, hence, the accretion flow should inflate and be geometrically thick and optically thin. In the high state, radiation can efficiently remove thermal energy, allowing the disk to take on an optically thick and geometrically thin configuration like the $\alpha$-disk of Shakura & Sunyaev (1973). Since the disk can cool efficiently, there is no disk truncation and the thin disk extends down to the ISCO.

The bifurcation of accretion states in active galactic nuclei (AGNs) is less clear, but there is evidence that accretion in low-luminosity AGN (LLAGN), including Sgr A* in our galactic center, is dominated by a hot mode (Rees et al., 1982; Narayan & Yi, 1995a). In AGNs whose accretion rate is sufficiently low, less than $1\%$ of their Eddington limit, hot mode accretion (Ho, 2008; Yuan & Narayan, 2014) and inner disk truncation (Trump et al., 2011; Yu et al., 2011; Nemmen et al., 2014) seems to be common. Qualitatively, the properties of LLAGN resemble the BHB hard state. The SED of a standard, unobscured, luminous Seyfert ($1-30\%$ Eddington), on the other hand, typically has a large thermal component peaking in the far-UV, analogous to the soft state of BHBs. The prevalence with which thick-to-thin disk transitions are invoked to explain the phenomenology of both BHBs and AGNs warrants a rigorous investigation of the behavior of the accretion flow at this interface to solidify the theoretical underpinnings of the truncation model.

For gas to accrete onto the central black hole, it must lose angular momentum. The seminal work of Shakura & Sunyaev (1973) parameterized the strength of the internal stress relative to the gas pressure through the so-called “$\alpha$-parameter,” further described in Section 2.2. The source of this anomalous viscosity was assumed to be turbulence, but the driving mechanism was not elucidated until the magnetorotational instability (MRI) was proposed as a viable instability (Balbus & Hawley, 1991). The connection between MRI-driven turbulence and the $\alpha$-parameter has since been well-studied in global MHD accretion disk simulations (e.g. Armitage, 1998; Hawley, 2000; Beckwith et al., 2011; Flock et al., 2011) and thin disks have been shown to behave viscously.

The behavior, energetics, and dynamics of the hot phase accretion, on the other hand, are not well constrained. Many types of hot flows have been suggested to fill this role including advection dominated accretion flows (ADAFs; Narayan & Yi, 1994, 1995a, 1995b) and convection dominated accretion flows (CDAFs; Quataert & Gruzinov, 2000). Simulations of non radiative, thick disks with varying levels of sophistication have been used to shed light on the behavior of hot accretion flows. Both HD (Stone & Balbus, 1996; Igumenshchev & Abramowicz, 1999; Igumenshchev et al., 2000; Proga & Begelman, 2003; Yuan & Bu, 2010) and MHD (Stone & Pringle, 2001; Hawley et al., 2001; Hawley & Balbus, 2002; Igumenshchev et al., 2003) numerical models ubiquitously develop outflows. In the HD cases, and some early MHD models, circulation often developed within the outflow and is interpreted as convection. In more recent MHD models, though, the convection is found to be stabilized (e.g. Yuan et al., 2012; Narayan et al., 2012; Yuan et al., 2015).

Despite the thoroughness with which numerical models have investigated the thin and thick disk configurations separately, the details of how the hot and cold flow coexist are not understood. Several global accretion models have been proposed and either explored analytically or with simple numerical models to determine how a truncated disk might develop. For instance, a truncated disk could exist because the inner disk evaporates due to the cool, Keplerian disk being heated by a hot corona (Meyer & Meyer-Hofmeister, 1994; Narayan et al., 1996; Liu et al., 1999; RóżaŃska & Czerny, 2000). Alternatively, the diffusion of energy from the interior of the disk could also heat the gas beyond the viral temperature, thereby driving the inner disk to a hot accretion solution and truncating the disk (Honma, 1996; Manmoto et al., 2000; Manmoto & Kato, 2000). Local thermal instability can also develop if the $\alpha$-parameter is proportional to the magnetic Prandtl number (Potter & Balbus, 2014, 2017). However, there is still much uncertainty underlying the trigger mechanism.

Still, headway can be made in characterizing the flow dynamics of a truncated disk. Of the outstanding issues, there is a particular need to understand the structure of the transition from the assumed thin disk to the thick disk and the influence the truncation has on the disk evolution. For instance, how does the flow conserve mass, angular momentum, and energy when it passes through an abrupt transition zone? Can a steady-state configuration develop, or is time-variability a necessary ingredient? On these fronts, numerical models provide a valuable tool and way to glean insight into the properties of truncated disks.

Radiation HD (RHD; Das & Sharma, 2013; Wu et al., 2016) simulations have tried to directly address the viability of the state transitions of BHBs. They find that when realistic radiation prescriptions which account for bremsstrahlung, synchrotron, and synchrotron self-Comptionization are included in the models, disks that are qualitatively similar to the BHB state transition model readily develop. These models naturally form truncated, thin, cold disks embedded in a hot, diffuse medium. Recently, Takahashi et al. (2016) used a simulation which accounted for general relativity, radiation and MHD (GRRMHD) and found a truncated disk developed where gas resided in a two-phase medium configuration. In the inner region of their simulations, the gas is radiatively inefficient and the viscous time is shorter than the cooling time. While it was not explicitly stated, this presumably sets up an ADAF-like flow. Evidence for thermal, and possibly viscous instability, has also been found in the GRRMHD simulations in Mishra et al. (2016).

The truncated disks found in the RHD and GRRMHD simulations suggest these transitions can naturally happen, but relying on computationally intensive algorithms comes at the sacrifice of resolution and simulation time. A pervasive issue in these works is that the cold disk cools to the scale of the grid and the internal dynamics are unresolved. This poses a problem because the defining evolution of the accretion flow can occur on small scales and on long times, which makes a proper simulation of a truncated accretion disk very difficult to do with current computational capabilities. Thus, to conduct a detailed study of the accretion flow itself and get a handle on how a truncated accretion disk will evolve, we must make the problem tractable by forgoing some of the more computationally expensive physics, like a full treatment of the radiation physics and the general relativistic effects which may impact the evolution of the disk in the very inner region close to the black hole. Here, we assume the truncation found in the simulations of Das & Sharma (2013); Wu et al. (2016); Takahashi et al. (2016); Mishra et al. (2016) can naturally occur and seek to fully characterize the accretion flow of a truncated $\alpha$-type viscous disk and understand the observational consequences.

In this paper we produce the first sustained truncated accretion disk, i.e. a truncation that can be studied for thousands of dynamical times. Our simulation is well-resolved and allows for a detailed study of the dynamics, energetics, and angular momentum transport in the disk. We do not use a full prescription for the radiation, but rather use an ad hoc, bistable cooling function to drive a thermal instability interior to some radius, thereby creating a transition between a cold and hot disk that can be characterized. We emphasize that the aim of this work is not to reproduce the state transitions observed in BH systems, instead it is to study the putative truncated disk found in hard state BHBs and LLAGNs. To isolate the behavior of the transition zone, we study the most ideal scenario possible, an axisymmetric viscous HD accretion disk in a Newtonian potential. The work presented here (Paper I) is a precursor to a study using a 3D high-resolution MHD simulation with the same bi-stable cooling function (Paper II). Together, these investigations will significantly improve our understanding of truncated disk dynamics.

The viscous formalism used here gives the freedom to set the viscosity and balance the heating rate to the cooling rate according to an $\alpha$-parameter. To better match the accompanying MHD truncated disk simulation, we set $\alpha=0.065$. This value was determined through a calibration MHD simulation, described in Paper II. In setting the HD $\alpha$-parameter to that of an effective $\alpha$-parameter from an MHD model, we can provide a more direct association of the disk behaviors between the HD and MHD regime.

By first establishing a benchmark with a simple $\alpha$-disk, we hope to gain a better understanding of the more complex disk model set in MHD, where heating and angular momentum transport occur solely thorough MRI-driven turbulence, and more honestly assess the viability of truncated disk models in the literature. Much of the previous theoretical work is developed in HD, so the simplified model considered here will help bridge the gap between expectations of models developed in the literature and the more MHD model.

This paper is organized as follows. In Section 2 we discuss our numerical model, simulation setup, and bistable cooling function. In Section 3 we discuss our analysis of dynamics, angular momentum transport, and energetics of the truncated disk. In Section 4 we discuss the implications of our results and the caveats associated with the simulation. In Section 5 we provide a summary and closing remarks.

2. Numerical Model

The goal of our work is to emulate a truncated accretion disk and perform a detailed study of its dynamics, angular momentum transport, and energetics. We therefore idealize our test problem so that we can direct our computational resources to the physics necessary to capture the hot-cold flow boundaries, run the model at high resolution, and run for many evolutionary times. In this work, we consider a two-dimensional, axisymmetric, viscous disk. The simulation is divided into two phases: an initialization period before the switch to the bistable cooling function and the latter portion used in our analysis where the bistable cooling function drives thermal instability and creates a sustained truncation. For numerical reasons that become especially important in the MHD case, we initially evolve the simulation as a thin disk with constant aspect ratio $h/r=0.1$ where the disk thickness is maintained by a cooling function with a single target disk thickness. Once a quasi-steady state develops, the cooling is switched to a bistable cooling function where the gas can cool efficiently with a target scale height $h/r=0.1$ or inefficiently with a target scale height of $h/r=0.4$, depending on its temperature. If the gas temperature is above the transition temperature, $T_{T}$, it cools according to the inefficient prescription; if it is below it cools according the the efficient prescription. The analysis focuses on the truncated disk.

The additional computational expense of accounting for full GR in this hydrodynamic simulation is manageable; however, it would be very restrictive for our follow-up simulation in MHD and hinder our ability to achieve the goal of performing a detailed study of the flow behavior. Therefore, we opt to consider the simplified Newtonian case here to aid in comparing between the HD and MHD scenarios. Additionally, these transitions are expected to occur at $30-100$ $r_{g}$ (Nemmen et al., 2014), which is far enough from the black hole that relativistic effects are minimal, so there is observational motivation for neglecting GR as well.

2.1. Simulation Code

For this study we use the second-order accurate PLUTO v4.2 code (Mignone et al., 2007) to model the gas dynamics. PLUTO employs a Godunov scheme that evolves the equations of hydrodynamics,

[TABLE]

where $\rho$ is the gas density, ${\bf v}$ is the fluid velocity, $P$ is the gas pressure, ${\cal I}$ is the unit rank-two tensor, $\sigma$ is the viscous stress tensor, $E$ is the total energy density of the fluid,

[TABLE]

and $\Lambda$ accounts for radiative losses through cooling. PLUTO was used in the dimensionally-unsplit mode with the hllc Riemann solver (a three-wave, four state solver) with piecewise linear spatial interpolation and second-order Runge-Kutta time integration. The Godunov scheme is a finite-volume method which conserves the total energy of simulation to machine accuracy. The only energy losses come from outflow across the domain boundary and the cooling function of the simulation.

2.2. Simulation Setup

In our simulation we adopt an axisymmetric, two-dimensional, spherical coordinate system ($r$, $\theta$, $\phi$) for our computational grid. The domain spans $r\in[10r_{g},1000r_{g}]$ and $\theta\in[\pi/2-1.0,\pi/2+1.0]$ with $N_{R}\times N_{\theta}=512\times 512=2.62\times 10^{5}$ zones. The radial spacing of the grid increases logarithmically so that $\Delta r/r$ is constant. In the $\theta$-direction, the simulation uses 320 uniformly grid zones in the midplane region from $\theta=\pi/2\pm 0.5$ and 96 geometrically stretched grid zones in each of the coronal upper regions ($\pi/2+0.5\leq\lvert\theta\rvert\leq\pi/2+1.0$).

Outflowing boundary conditions are used for the inner and outer $r$-boundaries and the $\pm\theta$ boundaries, but were modified to prevent the inflow of material. Periodic boundary conditions are used for the axisymmetric $\phi$-direction. A density and pressure floor is implemented to protect against values becoming artificially small or negative. All fluid variables (e.g. $\rho$, $p$, $T$) are reported in normalized units from the simulation.

We consider an ideal gas with a $\gamma=5/3$ adiabatic equation of state in a Newtonian gravitational potential where,

[TABLE]

Initially, the simulated disk is thin ($h/r=0.1$) and in vertical hydrostatic equilibrium. The steady state $\alpha$-disk solution is used for initialization:

[TABLE]

where $\rho_{0}$ is a normalization constant to set the maximum disk density to unity, $R$ is the cylindrical radius ($R=r\cos\theta$), $z$ is the vertical disk height ($z=r\sin\theta$), and $c_{s}$ is the sound speed ($c_{s}=\sqrt{P/\rho}=(h/r)v_{K}$). $R_{*}$ is taken to be just inside the inner simulation domain ($R_{*}=11\>GM/c^{2}$) to account for the torque-free inner boundary condition.

In astrophysical accretion disks, the effective viscosity originates from MRI driven turbulence. In our hydrodynamic simulation, however, we approximate this by treating the angular momentum transport with an $\alpha$-type viscosity according to Shakura & Sunyaev (1973). We set the shear viscosity term to be,

[TABLE]

where $\Omega_{K}$ is the Keplerian angular velocity. A dimensionless $\alpha$-parameter is set to $\alpha=0.065$ to match the test MHD run discussed in Paper II, as we described in Section 1.

2.3. Cooling Function

Our goal is to construct a model accretion disk that undergoes a radial transition from a geometrically-thin radiatively-efficient state to a geometrically-thick inefficient state as the matter flows inwards. Once such a state has been constructed, we can then characterize the dynamics of this transition region. Here, we describe the toy optically-thin cooling function that we use to achieve this goal.

The volumetric heating due to the $r-\phi$ shear stresses is

[TABLE]

where $\Omega_{K}$ is the Keplerian angular velocity, and ${\cal R}$ is the gas constant. In a steady-state, radiatively efficient disk we will have ${\cal H}=\Lambda$. We seek to construct a cooling function $\Lambda(\rho,T,r)$ that (i) is positive-definite everywhere (in contrast to many of the “target thickness” cooling functions employed in previous work; Noble et al. 2009), (ii) yields a disk model that is viscously stable, (iii) produces a hot, geometrically-thick solution at small radii, (iv) produces a cool, geometrically-thin solution at large radii. We adopt a cooling function of the form,

[TABLE]

where $f(r)$ shall be specified later to maximize convenience. The linear dependence on density ensures that the thermal balance is independent of density, removing the possibility of viscous instability. Noting that $c_{s}^{2}={\cal R}T=(h/r)^{2}\Omega_{K}^{2}$, we see that the condition ${\cal H}=\Lambda$ gives

[TABLE]

We choose $f(r)={\cal R}^{2}r^{-2}\Omega_{K}^{-1}\propto r^{-1/2}$ so that

[TABLE]

and, by construction, the disk has a constant aspect ratio (i.e., fixed $h/r$) if $\Lambda_{0}$ is a constant.

In this work, we induce thermal instability and hence a radial transition into the flow by introducing an abrupt decrease in $\Lambda_{0}(T)$ at a prescribed transition temperature $T_{T}$,

[TABLE]

where $\Lambda_{c}$ and $\Lambda_{h}$ are constants, and we will assume here that $\Lambda_{c}>\Lambda_{h}$. With $\Lambda_{0}$ so defined, the thermal balance condition ${\cal H}=\Lambda$ delineates three radial regions within the disk. Inside of a radius $r_{tr,i}=3\alpha GM/2\Lambda_{c}{\cal R}^{2}T_{T}$ (the inner transition radius), local thermal balance requires that the flow is hot with geometric thickness

[TABLE]

Outside of a radius $r_{tr,o}=3\alpha GM/2\Lambda_{h}{\cal R}^{2}T_{T}$ (the outer transition radius), thermal balance requires that the flow is cold with geometric thickness

[TABLE]

Between these two transition radii, the flow is bistable and thermal balance can be achieved either in the cold or the hot branch with geometric thickness

[TABLE]

Note that, due to the form of our cooling law, the ratio of the cold/hot geometric thickness is related to the ratio of the inner/outer transition radii,

[TABLE]

In this work, we set

[TABLE]

which implies $r_{tr,o}=1600r_{g}$. As we shall see, this toy cooling function results in a successful transition from a cold/thin to a hot/thicker disk close to the inner transition radius. Interestingly, the flow in the bistable region remains predominantly in the cold/thin state.

3. Results

We will now discuss our analysis of the hydrodynamic disk simulation. We show that a truncated disk does indeed develop and investigate the dynamics, energetics, and angular momentum transport associated with the accretion flow in the truncated state. Interestingly, we find that a sustained stream of sub-Keplarian gas flows inwards from the thick-to-thin disk transition and undergoes delayed heating. As the gas heats to $T_{T}$ and beyond, pockets of hot gas form in the dense disk which buoyantly rise and cool, inducing a circular flow in the coronal regions of the disk and CDAF-like behavior. As the disk refills at the transition radius, the cycle repeats quasi-periodically.

3.1. The Truncated Disk

We begin by discussing the initial development of the truncated disk in the final phase of the initialization and the general properties of the truncated disk. As discussed in Section 2.3, the simulation is first evolved with constant $h/r=0.1$, from $t=0$ $GM/c^{3}$ to $t=1.98\times 10^{6}$ $GM/c^{3}$ (the first $10,000$ inner orbits) before the cooling is switched to include both a hot and a cold branch. Immediately after the switch, there is a large, transient heating event within $r\approx 200\>r_{g}$ that lasts roughly $\Delta t=1.4\times 10^{4}$ $GM/c^{3}$ and then cools over $\Delta t=8.6\times 10^{4}\>GM/c^{3}$. At $100$ $r_{g}$, $t_{dyn}\sim 6.3\times 10^{3}GM/c^{3}$ so the heating event lasts approximately 2.5 dynamical times and the cooling event occurrs on roughly the thermal time of the inner disk. Figure 1 shows four snapshots of $\rho$ during the heating event.

The gas in the inner region of the disk responds to the change in cooling by inflating since it is beyond $T_{T}$ and has an excess of thermal energy which it cannot rapidly dissipate on the inefficient cooling branch. Additionally, viscous stresses continue to do work on the gas, further increasing the temperature. A thin, dense stream of material plunges into the inner region between $t=2.01\times 10^{6}\>GM/c^{3}$ and $t=2.04\times 10^{6}\>GM/c^{3}$. Since this gas is near $T_{T}$, it ablates away, ultimately leaving a truncated disk. The final remnants of the transient can be seen in a very thin stream of dense gas in the inner region of the disk at $t=2.01\times 10^{6}\>GM/c^{3}$ and large, buoyant bubbles rising away from the transition radius.

The abrupt change in the cooling law is, obviously, unphysical and the heating event is the reaction of the gas to the sudden switch. Following this, however, the disk quickly settles into a new, sustained configuration with a “transition zone” from $r=80-100\>r_{g}$. In this transition zone, the region is cyclically filled with material through accretion, but then evacuated as the gas heats to $T_{T}$. Unless explicitly stated, the analysis we detail is conducted following the implementation of the bistable cooling function.

Figure 2 shows two representative snapshots of the gas density and scale height temperature during initialization and during the analysis. Before the change in cooling function, the disk has a uniform scale height ratio and very little hot gas is present in the simulation. In the density and temperature maps of the truncated disk, wispy regions of cold, dense gas are seen surrounding diffuse, hot bubbles. The cold features are very dynamic and short lived, as they are close to the transition temperature and can spontaneously transition to the hot phase.

Figure 2 also shows the time averaged density of the disk with two measures of the characteristic disk scale height overlaid. The first measure is that of the geometric scale height ratio,

[TABLE]

where $d\Omega=\sin\theta d\theta d\phi$ is the solid angle element in spherical coordinates and ,

[TABLE]

is the average polar angle of the gas. The second measure is the density weighted, temperature scale height ratio,

[TABLE]

where $T$ is the gas temperature, $T=P/\rho$. The values of $h_{g}/r$ and $h_{T}/r$ were calculated for every data dump and time averaged, denoted by the angled brackets. Throughout this paper we take gas in the “hot phase” to be that with $T>T_{T}$ and gas in the “cold phase” to have $T<T_{T}$.

Spatially, the time-averaged gas density shows a distinct truncation at $r=100$ $r_{g}$, which also corresponds to location where the disk becomes thicker, as shown by $h_{g}/r$ and $h_{T}/r$. In the inner region of the disk $\rho$ is lower than the disk body by approximately an order of magnitude. The averaged scale height ratio of the inner disk is only twice as thick as the rest of the disk and does not reach the “thick” scale height target of $h_{h}/r=0.4$ due to a significant amount of advective cooling.

The signature of the thick-to-thin transition is also present in the thermal profile of the disk. Figure 3 shows the time averaged radial profile of the fraction of gas in the hot and cold phases. The fractions do vary in time, but, nevertheless, from this profile, we can see hot phase gas typically resides in the inner region of the disk and cold phase gas resides beyond the density transition. It is also worth noting that the transition from a predominately hot to cold phase, i.e. where the fraction gas in the hot phase equals the fraction of gas in the cold phase, does not occur at $r=100\>r_{g}$ as we would expect from the cooling function, but rather just inwards of this radius at $r=80\>r_{g}$. In brief, this is due to delayed heating as the gas is close to, but not yet at, $T_{T}$ as it crosses $r=100\>r_{g}$. At the hot/cold flow interface it moves inwards and then makes the transition within the expected transition radius.

The truncation is also present in the time averaged surface density,

[TABLE]

The hot and cold gas phases as a function of radius is also shown in Figure 3. Compared to the surface density at the end of the initialization, the inner disk surface density is lower by a factor of five. Further out, beyond the transition radius where the gas is predominantly in the cold phase, the overall surface density is very near the intialization level. Thus, through the launching of the convective cells by the transitioning gas, there is essentially an evacuation of the inner disk relative to the untruncated disk.

3.2. Disk Dynamics

The disk dynamics in a truncated disk are of interest since hot and cold accretion flows should behave quite differently, but the truncated disk model used to explain astrophysical sources expects the two flows to exist in close proximity. In the canonical thin accretion disk, the disk is rotationally supported and thermal motions are a small component of the overall velocities, i.e. $h/r=c_{s}/v_{K}$ where $v_{K}$ is the Keplarian velocity (Shakura & Sunyaev, 1973).

In Figure 4 the time averaged ratios of $\langle v_{\phi}\rangle/v_{K}$ and $\langle v_{r}\rangle/v_{K}$ are shown as a function of position in the simulation domain. In the normalized simulation units (i.e. $G=M=c=1$), the Keplerian velocity is:

[TABLE]

The disk establishes four dynamically distinct regions. First is the outer disk body where the gas follows Keplerian orbits with very small radial velocity, as in a standard thin accretion disk. Second, are the coronal regions, where the gas is sub-Keplerian with an average azimuthal velocity of $\approx 15\%$ below $v_{K}$. For comparison, previous hydrodynamic studies of hydrodynamic CDAF-like flows (e.g. Stone et al., 1999), typically find that the azimuthal velocities are nearly-Keplerian and greater than we find here. Additionally, as can be seen in the plot of $\langle v_{r}\rangle/v_{K}$, there is a substantial outflow with $v_{r}\approx 10\%$ of the local $v_{K}$ immediately above and below the disk, and a comparable fallback velocity at higher altitude. This is due to a single convective cell that occupies each of the coronal regions. Third is the inner disk body found within $r=100\>r_{g}$ and one scale height above and below the disk. In Figure 2 we saw that this region was less dense, but contained a fair amount of cold gas. In Figure 4 we can see that this cooler region is modestly sub-Keplerian and inflowing at $1\%$ of $v_{K}$. Finally, there is a thin layer of gas enveloping the disk at one scale height where the $\langle v_{r}\rangle/v_{K}$ map reveals that there is a prominent inflow directly on the disk body. It is in this region that the bulk of the mass accretion from the outer disk occurs.

The launching mechanism is interesting, especially considering this is a simple hydrodynamic simulation so a magnetocentrifugally driven wind (Blandford & Payne, 1982) or a radiation driven wind (Proga et al., 2000; Proga & Kallman, 2004) cannot form. Here, we find the outflow in the simulation is launched by pockets of hot gas that spontaneously undergo a transition from the cold phase to the hot phase in the cool stream. As the gas expands, it becomes less dense and, thus, buoyantly unstable according to the the Solberg-Høiland criterion describing the stability of a rotating fluid. When the gas is launched, there is already an angular momentum deficit, so it has little rotational support at a large radius. Unlike a pure CDAF, the gas accretes when it falls back because, while the cooling of the hot gas is inefficient, it is not zero and Bernoulli parameter is not conserved.

The launching processes is shown in a sequence of images in Figure 5. The cold, sub-Keplerian stream can be seen in the dense gas. As the gas spontaneously transitions to the inefficient cooling branch, it heats and expands. In doing so, these pockets of newly heated gas are surrounded by denser, cooler material so they buoyantly rise which disrupts the cold stream and disturbs the flow. Once the gas reaches a low density region, at approximately one density scale height, it moves radially against the temperature-gravity gradient.

Figure 6 shows a time trace of the fraction of the gas in the hot phase, $f_{h}$, and the fraction of gas in the cold phase, $f_{c}$, at $r=100\>r_{g}$. The transient heating event is seen at early times in time traces of $f_{h}$ and $f_{c}$. Once the disk relaxes into its truncated configuration, $f_{h}$ and $f_{c}$ oscillate as heating events quasi-periodically occur in the cold streaming material and hot gas bubbles rise. These fluctuations are relatively small, but rapid with the fluctuations on average taking $t=2.65\times 10^{4}\>GM/c^{3}$, or four thermal times. At $r=100\>r_{g}$ the average gas fractions in each phase are $f_{c}=0.84$ and $f_{h}=0.15$. The oscillations in the gas fraction are asymmetric and have a standard deviation of $0.03$. As the oscillations occur, $f_{h}$ rises rapidly and then experiences a slower decay. Consequently, $f_{c}$ decreases sharply and gradually increases before the cycle repeats.

3.3. Mass Accretion and Angular Momentum Transport

Previous numerical studies of convective-like, hydrodynamic accretion flows have found that very little gas actually accretes onto the black hole (Igumenshchev & Abramowicz, 1999; Stone et al., 1999; Igumenshchev & Abramowicz, 2000; Narayan et al., 2000). Here, we report conflicting results, as we find that the hot gas does indeed accrete in our simulation. One possibility for this divergence from previous studies of hot flows is that the convective cells in our simulation are continually fed at the truncation of the thin disk so the gas originates from deeper in the the gravitational potential, rather than being a hot flow throughout the entire domain and overall being less gravitationally bound.

In our 2D, axisymmetric disk, the instantaneous mass accretion rate in the simulation is given by,

[TABLE]

Using our definition of the gas phases from Section 3.1, we show the time trace of the accretion rate for the hot ($\dot{M}_{H}$) and cold phase ($\dot{M}_{C}$) gas at the inner boundary of our simulation in Figure 7. The pre-analysis disk mass accretion rate of cold gas, $\langle\dot{M}_{pre}\rangle=6.65\times 10^{-3}$, is also shown for reference.

We find, somewhat surprisingly, the mass flux across the inner boundary is slightly higher in the truncated disk than in the initialization phase. Additionally, there is a striking disparity between the $\dot{M}_{h}$ and $\dot{M}_{c}$. Overwhelmingly, by approximately three orders of magnitude, hot gas is accreted rather than cold gas. Initially, the overall mass accretion rate is higher during the transient phase of the heating event, but it settles into a sustained accretion rate of $\langle\dot{M}\rangle=1.0\times 10^{-2}$ for the remaining length of the simulation.

Since an appreciable amount of gas is accreted, there must be efficient angular momentum transport in the simulation. In our disk, the angular momentum can be transported three ways: through advection, stresses arising from the radial shear of the disk, and stresses from shear in the atmosphere of the disk as it rotates within the sub-Keplerian convective outflow. The conservation of angular momentum equation takes the following form,

[TABLE]

where the three components appearing in the divergence term are, respectively, the advection term, the $r-\phi$ viscous stress term, and the $\theta-\phi$ stress term. The time-averaged maps of these terms shown in Figure 8.

The importance of the transport terms at different locations in the disk is quite clear. The $r-\phi$ viscous stress is negative everywhere and largest within the disk body, closer in towards the inner domain. This structure is expected since we consider an $\alpha$-type viscosity where the viscous stress is proportional to the gas pressure and $\phi$-shear and the pressure is largest in the inner regions of the disk. Above and below the disk, at approximately one scale height, there are localized layers of $\theta-\phi$ stress where the Keplerian disk interacts with the sub-Keplerian outflow. Also, within the very inner regions, the sign of the stress in each respective hemisphere reverses where the convective cell closes, reflecting the opposite signs of the vertical $v_{\phi}/\sin(\theta)$ gradient. The largest $\theta-\phi$ stresses occur on the upper and lower edges of the truncated disk where the cold, dense disk is buffeted by the outflowing, hot, buoyant gas. The stress gradually decreases at larger radii as the strength of the interaction between the disk and outflow decreases.

The time averaged map of the angular momentum advection terms is also stratified, reflecting the structure seen in Figure 4 of the radial velocities. It is worth noting that within the disk body meridional circulation establishes a very small, but positive, net radial velocity at the disk midplane. This behavior was originally predicted by Kippenhahn & Thomas (1982) and Urpin (1984) and has since been reproduced numerically for viscous accretion disks (Kley & Lin, 1992; Rozyczka et al., 1994). Moving away from the disk midplane to higher latitude, the circulation reverses the flow and gas moves inwards.

The radial profiles of the time averaged angular momentum fluxes resulting from the terms in Equation 26 are shown in Figure 9. On average, the viscous $r-\phi$ stress is balanced by advection, at least within the region that has established a quasi-steady state and viscous equilibrium (out to $\approx 500\>r_{g}$). In a purely convective flow, angular momentum is advected inwards, rather than transported outward (Ryu & Goodman, 1992; Stone & Balbus, 1996). However, in this model, the viscous stresses continually act on the gas, allowing for the loss of angular momentum. Consistent with the removal of angular momentum by viscous torques, the specific angular momentum profile, also shown in Figure 9, is Keplerian, $l\propto r^{1/2}$, when averaged over the duration of the simulation.

The angular momentum transport at the hot-cold phase interface is complex and depends on the time variability of $\Sigma$, $T$, the $r-\phi$ stress component, and the advective stress component. Figure 10 shows these quantities at $r=80r_{g}$, in the heart of the transition zone, over several of the bubble launching cycles. In Figure 10, $\Sigma$, $T$, the $r-\phi$ stress, and advective stress have been averaged over one thin-disk scale height ($h=\pi/2\pm 0.1$).

At the truncation, a cyclical building and expulsion of material occurs with a period of $t\approx 2.65\times 10^{4}\>GM/c^{3}$. In this cycle, peaks in $T$ occur when the density is lowest and vice versa. The angular momentum terms vary with $\Sigma$ and $T$, but there is a slight delay, $\Delta t\approx 2.05\times 10^{3}\>GM/c^{3}$, in the maxima and minima. For the majority of the cycle, approximately $80\%$, the surface density increases. Eventually, a critical threshold is reached when the heating rate exceeds the cooling rate, causing the temperature to increase and the cooling to become inefficient as pockets of gas reach $T_{c}$. This causes a runaway heating event, seen in the rapid increase of $T$, and the outflow to be launched which can be seen as turn over in the advection term to positive values.

A key feature throughout this cycle is that the gas pressure remains relatively constant as the gas expands to maintain a pressure equilibrium, a facet which is discussed to greater depth in Section 3.4. This, in turn, means that the specific angular momentum of the lower density, transitioning gas at the truncation is more rapidly lost to the neighboring dense gas of the disk body through the $r-\phi$ viscous stress since the amplitude of the viscosity is proportional to the pressure. Additionally, as seen in Figure 5, the hot transitioning gas is in contact with the dense gas of the cooler stream which cannot lose angular momentum as well. As the specific angular momenta of these two regions changes due to the differences in the viscous efficiency, a vertical shear is introduced. This shear allows the $\theta-\phi$ viscous stress to transfer angular momentum from the dense gas in the stream to the diffuse, outflowing bubble where it is carried away. This changes the shear rate of the dense gas, resulting in a spike in the $r-\phi$ stress component. This rapid exchange of angular momentum drives the angular momentum transport in the stream which then flows to smaller radii, as evidenced by the sharp negative dip in the advective angular momentum transport term corresponding to the peak in the $r-\phi$ stress term.

After this involved transitory event, the disk refills and the cycle repeats.

3.4. Disk Energetics

3.4.1 Energy Terms

In our hydrodynamic model, the conservation of energy equation is,

[TABLE]

where $\mathcal{E}$ is the total disk energy, $\mathcal{C}$ is the “radiative” cooling done by the cooling function, $\mathcal{A_{K}}$ is the advection of kinetic energy, $\mathcal{A_{T}}$ is the advection of thermal energy, $\mathcal{G}$ is the advection of gravitational potential energy, and $\mathcal{W}$ is the work done by viscous torques. We will consider these terms in a spherical shell integrated form, such that they are

[TABLE]

The time averaged, radial profiles of these terms are shown in Figure 11. The presence of a truncated disk influences every energy term, with the exception of $\mathcal{G}$. Given that there is sustained accretion and the average specific angular momentum profile is Keplerian, this is not necessarily surprising in the time average profile. Essentially, the smooth, time averaged $\mathcal{G}$ indicates that, globally, material is accreting and there is no mass build-up anywhere in the disk. Even though material may temporarily be caught in convective cells, it is eventually accreted so gravitational potential energy is liberated.

In the outer disk the $\mathcal{A_{T}}$ term, which is positive, is suppressed and accounts for very little of the overall energy transport of the gas. This changes within $r<R_{T}$, however, as $\mathcal{A_{T}}$ and $\mathcal{A_{K}}$ become major contributors in the truncated disk region. Advection plays a dominate role in the energy transport within the transition radius, overtaking the cooling. The dynamical influence of the convective cells acts to advect energy across the inner boundary as material gets swept inwards.

There is a slight deficit in the $\mathcal{A_{K}}$ term and a significant deficit in the $\mathcal{W}$ term inside of $r<R_{t}$ because of the disk truncation. The impact on the $\mathcal{A_{K}}$ term is twofold. In the inner disk the density is lower and the flow is sub-Keplerian, thus overall advection of kinetic energy in this region is slightly less what it might otherwise be. Regardless, the sign does not change as the bulk of the kinetic energy is swept inwards in the dense gas that streams inwards through the viscous interaction of the hot outflow. In a general, time averaged sense, the energy dissipation of the $\mathcal{W}$ term is suppressed because the density is lower and the fluid motions are more chaotic in this region. However, as seen in the variability of the viscous stress components, there are brief periods of enhanced viscous stress and, therefore, heating. Nevertheless, the disruption of the flow by the hot gas bubbles makes it more difficult to establish the sustained radial and vertical gradients in $v_{\phi}$ that generate the viscous stresses and heating in the outer disk, thereby diminishing the contribution of the $\mathcal{W}$ term in the inner region of the disk.

The transition has the strongest affect on the $\mathcal{A_{T}}$ term, as its sign actually flips at the disk transition. Within the inner disk, material is generally swept inwards, either in the streaming material or through the return of the convective cells. The sign of $\mathcal{A_{T}}$ flips at the transition, however, because the majority of the thermal energy resides in the buoyant, outflowing gas. The apparent suppression of this term is somewhat artificial because it accounts for two competing components; the net inflow of thermal energy through the accretion in the dense, thin disk and the outflow of the diffuse, hot gas.

Further insight into the periodic launching of the buoyant bubbles can be gained from examining the individual components of $\mathcal{E}$, i.e. the kinetic energy (\mathcal{E}_{K}=1/2\rho\hbox{{\bf v}}), the internal energy ($\mathcal{E}_{I}=P/(\gamma-1)$), and the gravitational potential energy ($\mathcal{E}_{G}=\rho\Phi$) of the gas. The variability of these energy components at $r=80\>r_{g}$ is shown in Figure 12. Similar to Figure 10, these terms have been averaged over $h=\pi/2\pm 0.1$.

The cyclical behavior observed in $\Sigma$, $T$, the $r-\phi$ stress, and the angular momentum advection is also seen in the gas energy. Unlike the involved angular momentum transport interaction observed in this region, the temporal variability of the gas energy terms is simply due to the removal of gas in the buoyant bubbles, which can be seen in the corresponding dips in $\Sigma$ in Figure 10 and in the variability of $\mathcal{E}_{G}$ since $\Phi$ is a constant at $r=80\>r_{g}$. The removal of the gas is tied to an increase in the gas energy as it heats (i.e. increased $\mathcal{E}_{I}$). Throughout the launching of the hot bubbles, the internal energy and gas pressure have very little variability ($\approx 2\%$). When gas transitions in the flow, the bubble expands adiabatically as it buoyantly rises vertically along the density gradient.

3.4.2 Bound Outflow

Throughout the simulation, the Bernoulli parameter,

[TABLE]

where $W$ is the specific enthalpy

[TABLE]

of the gas, is negative. When $Be$ is negative, the kinetic and thermal energy of the gas is less than the gravitational potential energy, so the gas is bound and cannot escape.

Shown in Figure 13 is the distribution of outflowing mass fraction ($\dot{M}_{out}$), i.e. gas with $v_{r}>0$, at $r=100\>r_{g}$ as a function of $Be$. On the surface, it looks like $\dot{M}_{out}$ is nearly independent of $Be$. When we dig deeper into the distribution, we find that, in fact, there are two components to the outflow: one component that is diffuse gas with high velocities and other that is dense gas with low radial velocity. Also shown in Figure 13 are the average density of the outflow and velocity as a function of $Be$. When a bubble of hot gas is launched from the inner region of the disk, the gas has a large, positive radial velocity, but low density. As it rises, higher density gas gets entrained in the wake of the bubble and pulled along to larger radius, although at a lower velocity.

3.4.3 Synthetic Light Curves

To finish our analysis, we use the cooling in our simulation to generate a synthetic light curve. Since PLUTO conserves total energy in the simulation, the “radiative losses” our cooling function imitates can be used as a proxy for the broadband emission. There are obvious caveats to this exercise, namely that we use an ad hoc cooling function to approximate the real, detailed radiation physics of the disk. Regardless, it will provide a sense of what type of photometric variability could originate from this disk.

Observationally, in astrophysical black hole systems, the higher energy emission is attributed to the hot, radiatively inefficient gas. The lower energy, blackbody emission is attributed to the cooler, radiatively efficient, dense gas. Following this paradigm, we use our previously defined temperature criteria to generate “hard,” $L_{h}$, and “soft,” $L_{s}$, band light curves, shown with the volume integrated cooling as a function of time of the hot and cold gas, respectively, in Figure 14.

During the initialization, the large heating event causes an immediate spike in $L_{h}$ that decays down to a level of $L_{h}=3.0\times 10^{-3}$ which it maintains for the rest of the simulation. This is higher than the total pre-analysis emission level, $L_{p}$, of $2.8\times 10^{-3}$. $L_{s}$ is approximately a quarter less than the $L_{h}$, due in part to the fact that the hot gas fills a large portion of the simulation volume. Interestingly, there is very little sign of the oscillatory behavior seen Figure 3 of the gas fraction at $r_{trans}$. Instead, the light curves of each gas phase qualitatively resemble the the time trace of the respective mass accretion rates, shown in Figure 7.

4. Discussion

By introducing a bistable cooling law to a simple hydrodynamic viscous disk, we find that gas transitioning from the hot phase to the cold phase develops significant time variability which affects the evolution of the disk. The transition between efficient and inefficient gas cooling leads to disk truncation and causes large, convective cells to develop in the inner region of the disk. The presence of sub-Keplerian, outflowing gas in close proximity to the accretion disk generates a large $\theta-\phi$ viscous stress which has previously been neglected, but is influential here. The inner, truncated region of the disk is occupied by hot gas that cools inefficiently. Beyond the transition, the bulk of the gas resides in cool phase in a configuration similar to a standard Shakura $\&$ Sunyaev thin disk. In general, the truncated disk in our simulation is qualitatively similar to the phenomenological model invoked to explain the spectral properties of BHBs in the low/hard and intermediate states like Cygnus X-1 (Done & Zycki, 1999), XTE J1650-500 (Done & Gierliński, 2006), H1743-322, (Sriram et al., 2009), MAXI J1543-564 (Rapisarda et al., 2014), and GX 339-5 (Plant et al., 2015), which was expected by construction. There are several notable differences between the simulation and the standard truncated disk model, which could manifest as observable signatures in BHB and LLAGN systems.

In the canonical HD truncated disk model, the inner accretion flow is assumed to be an ADAF and many of the observational results are interpreted with this assumption (e.g. Gammie et al., 1999; Shields et al., 2000; Ho et al., 2000; Ulvestad & Ho, 2001; McClintock et al., 2003; Yuan et al., 2005). Hence, it is interesting that a flow which closely resembles a CDAF naturally forms. Even though there are qualitative similarities, the convective flow is ultimately distinct from the standard CDAF. The primary way in which it differs is that the gas does not stay in the convective flow indefinitely as it is constantly cooling. Once it is launched, the gas rises to larger radii, but then stalls and falls back. The accretion rate is higher following the transition because of the interaction of the outflowing gas and the thin disk body. It is, in fact, primarily the gas in the hot phase that accretes with very little cold gas making it to the inner boundary of the simulation.

One might hope the oscillatory gas fraction could be connected to observable variability, like a QPO, but that is not currently supported by our simulation. QPOs are present in the intermediate states between the low-hard state and the high-soft states when the disk is assumed to be truncated, similar to what we modeled. However, neither the mass accretion rate nor our simple cooling prescription show appreciable variability. If the buoyant bubbles we capture in our simulation can be physically realized and related to QPO behavior, their signature must ultimately arise from more complex radiative processes than we have included.

Despite the obvious issues with connecting this HD accretion disk model to real black hole disks, it still plays an important role in providing a needed benchmark against which future MHD models of truncated disks can be compared, as we originally motivated. In particular, the surprising complexity in the disk structure, angular momentum transport, and energetics demonstrate how poorly the behavior of disks in this state is understood, even in what should be a simple case. The results from this simulation help to provide a bridge between the simple disk model currently used to understand state transitions in BHBs and LLAGNs and the more complex MHD physics that actually govern their evolution.

Of the many caveats and pitfalls associated with using an HD model to describe a real accretion disk, we will briefly discuss two. The first concern is that the behavior is sensitive to the strength of the anomalous viscosity. As the viscosity of the disk changes, the ability of the fluid to develop small scale structure and buoyantly rise also changes, which has been shown to affect the size and behavior of hydrodynamic, convective accretion flows (Igumenshchev & Abramowicz, 1999). To test the sensitivity of this model we reran the full simulation with $\alpha=0.04$. The development of this low $\alpha$-parameter was similar. Granted, the amplitudes of the stress, $\dot{M}$, etc. changed because of the lower viscosity, but the distinguishing flow properties like the rapid inflow of gas along the skin of the disk, a single large convective eddy in the atmosphere, and the quasi-periodic development of buoyant bubbles at the truncation were all reproduced. In all, we can be confident that the behavior we present in this paper is robust for reasonable $\alpha$-parameter values.

The biggest concern is the role MHD will play in altering the behavior of the disk. One of the most conspicuous features of our simulation is the presence of a “cold” stream of gas extending from the transition radius down to near the inner boundary. Observationally, there is evidence that cool gas may be present near the black hole, even when the disk is thought to be truncated (Miller et al., 2002a, b). While it is tempting to make a connection between what develops in our simulation and what is empirically observed, this must be done cautiously. In an astrophysical accretion disk the gas motions are dominated by vigorous, MHD turbulence with correlation lengths of order the disk scale height (Guan et al., 2009; Beckwith et al., 2011; Hogg & Reynolds, 2016). The stream of material in this simulation is smaller than that nominal scale, hence, it is reasonable to expect that even if this stream could form, it would quickly be disrupted by the turbulence and unable to exist in the way it does in our simulation. Additionally, if this cold stream is destroyed by MRI driven turbulence, it should follow that the convective cells would disappear since they are launched from bubbles of gas trapped in the cold stream. The absence of outflowing material against the disk would dramatically change the character from what we capture in this simulation since the vertical shear in the velocity field leads to the interesting behavior we observe.

Additionally, the properties of HD and MHD stresses are markedly different. In this HD model the stresses, and subsequent angular momentum transport, originate from a true viscosity. In an MHD disk, however, correlated radial and azimuthal fluctuations in both magnetic field and fluid velocity in anisotropic turbulence govern the angular momentum transport (Balbus et al., 1994). This casts doubt on whether the vertical velocity shear plays an important role in the disk evolution and angular momentum transport. Furthermore, the viscosity in this simulation does not vary. In MHD, the stresses fluctuate on a dynamical time and are modulated by dynamo behavior on tens of dynamical times. This natural variability in the effective viscosity could interact with the thermal instability and further manipulate the launching of hot bubbles of gas in the inner region of the disk. An examination of this model within a MHD framework in Paper II will address these concerns and provide a rigorous comparison between the HD and MHD regimes.

5. Conclusions

Using a viscous, hydrodynamic simulation of an accretion disk with a bistable cooling function, we conduct a detailed analysis of the dynamics, angular momentum transport and energetics. We find that a truncated configuration develops and that convective cells are launched from the evacuated region as cold gas flows inwards and undergoes delayed heating. Pockets of hot gas form within the cold gas stream and buoyantly rise along the surface of the cooler disk which introduces a $\theta-\phi$ viscous stress component through a vertical shear between the coronal atmosphere and the disk itself. This viscous stress component has been neglected in previous hydrodynamic treatments of accretion disks, but it plays in important role in the evolution of the disk by transferring additional angular momentum from the dense disk to the outflow and enhancing accretion in a thin layer long the surface of the disk. In this truncated disk, the bulk of the accreted material is in the hot phase, rather than the cold phase, and the accretion rate is slightly elevated compared to that of the thin disk during the initialization. The spontaneous transition of gas to the hot phase in the cold stream of material causes a quasi-periodic oscillation in the hot/cold gas fraction at the transition radius, but this does not translate in any appreciable fluctuations in the mass accretion rate, nor the volume integrated cooling ratio.

With this baseline hydrodynamic model of a truncated disk, we can begin to understand more complex future MHD models. Ultimately, the stresses we approximate as a hydrodynamic viscosity arise through the MRI, which behaves differently. Additionally, the interesting behavior displayed by this simulation develops through purely HD processes and on small scales. MHD turbulence should disrupt the flow enough that the global properties will be different; however, until a rigorous MHD test of a truncated accretion disk is completed, it is unknown how the disk evolution will be changed.

The authors thank the anonymous referee for useful comments that helped strengthen and clarify the paper. CSR thanks support from a Sackler Fellowship (hosted by the Institute of Astronomy, Cambridge) and NASA under grant NNX15AC40G. JDH thanks support by NASA under the NASA Earth and Space Science Fellowship program. The authors acknowledge the University of Maryland supercomputing resources (http://www.it.umd.edu/hpcc) made available for conducting the research reported in this paper.