Viscous Radiative Simulations of Thin Disks I.

# Relativistic Viscous Radiation Hydrodynamic Simulations of Geometrically Thin Disks: I. Thermal and Other Instabilities

P. Chris Fragile111Kavli Institute for Theoretical Physics, Santa Barbara, CA., Sarina M. Etheridge, Department of Physics and Astronomy, College of Charleston, Charleston, SC 29424, USA; fragilep@cofc.edu Peter Anninos, Lawerence Livermore National Laboratory, P.O. Box 808, Livermore, CA 94550, USA Bhupendra Mishra111Kavli Institute for Theoretical Physics, Santa Barbara, CA. JILA, University of Colorado and National Institute of Standards and Technology, 440 UCB, Boulder, CO 80309-0440, USA Włodek Kluźniak111Kavli Institute for Theoretical Physics, Santa Barbara, CA. Nicolaus Copernicus Astronomical Center, Bartycka 18, Warsaw, 00-716, Poland
###### Abstract

We present results from two-dimensional, general relativistic, viscous, radiation hydrodynamic numerical simulations of Shakura-Sunyaev thin disks accreting onto stellar mass Schwarzschild black holes. We consider cases on both the gas- and radiation-pressure-dominated branches of the thermal equilibrium curve, with mass accretion rates spanning the range from to . The simulations directly test the stability of this standard disk model on the different branches. We find clear evidence of thermal instability for all radiation-pressure-dominated disks, resulting universally in the vertical collapse of the disks, which in some cases then settle onto the stable, gas-pressure-dominated branch. Although these results are consistent with decades-old theoretical predictions, they appear to be in conflict with available observational data from black hole X-ray binaries. We also find evidence for a radiation-pressure-driven instability that breaks the unstable disks up into alternating rings of high and low surface density on a timescale comparable to the thermal collapse. Since radiation is included self-consistently in the simulations, we are able to calculate lightcurves and power density spectra (PDS). For the most part, we measure radiative efficiencies (ratio of luminosity to mass accretion rate) close to 6%, as expected for a non-rotating black hole. The PDS appear as broken power laws, with a break typically around 100 Hz. There is no evidence of significant excess power at any frequencies, i.e. no quasi-periodic oscillations are observed.

accretion, accretion disks — instabilities — X-rays: binaries

## 1 Introduction

From the earliest work on thin () accretion disks based on the -viscosity prescription (Shakura & Sunyaev, 1973), there have been notable problems whenever the steady-state solution calls for the disk to be radiation pressure dominated, as it will be above a certain mass accretion rate, , where g s is the Eddington mass accretion rate, and inside a certain radius, , where . Specifically, the solution is predicted to be both thermally (Shakura & Sunyaev, 1976) and viscously (Lightman & Eardley, 1974) unstable. The thermal instability arises owing to different dependencies of the disk heating, , and cooling, , on mid-plane temperature, , such that small deviations in the temperature are expected to lead to runaway heating or cooling. Similarly, the viscous instability arises from an inverse dependence of the vertically integrated stress, , on surface density, , leading to an anti-diffusion equation for such that the disk breaks up into alternating rings of high , low , and low , high .

This is particularly puzzling in light of observations of black hole X-ray binaries (BHXRBs). In these systems, the spectra look most disk-like (soft, with a prominent thermal bump around 1 keV) and stable (rms variability %) whenever (e.g. van der Klis, 2006; Done et al., 2007), corresponding to for a standard efficiency of , precisely when such disks are predicted to be unstable. So, the first puzzle one would like to solve is precisely when are these disks unstable.

The second point of interest is to find out what the non-linear manifestations of these instabilities are in real accretion disks and their consequences over periods comparable to a viscous timescale. One way to get a handle on these questions is through direct numerical simulations. At a minimum, this requires evolution of the equations of viscous, radiation hydrodynamics. Alternatively, since the true source of the “viscosity” in many accreting systems is thought to be turbulence driven by the magneto-rotational instability (MRI Balbus & Hawley, 1992, 1998), one could perform radiation magnetohydrodynamic (MHD) simulations. The codes required to perform such direct numerical simulations have only recently become available (e.g. Ohsuga et al., 2005; Krolik et al., 2007; Jiang et al., 2013; Sa̧dowski et al., 2013; McKinney et al., 2014), and we apply one such code, Cosmos++ (Anninos et al., 2005; Fragile et al., 2014), in this work.

The thermal stability of radiation-pressure-dominated accretion disks were first numerically studied in the context of limit-cycle-behavior using time-dependent, one-dimensional, height-integrated disk structure equations (Honma et al., 1991, 1992; Janiuk et al., 2002). Local, shearing box simulations that treat a small patch of the accretion disk, but include vertical and azimuthal structure, followed (Agol et al., 2001; Turner et al., 2002). Although there was an early claim of thermal stability (Hirose et al., 2009), the consensus now is that, indeed, radiation-pressure-dominated disks are thermally unstable (Jiang et al., 2013; Ross et al., 2017). We confirmed this behavior, and in addition found evidence for the viscous instability, in our own global radiation MHD simulations (Mishra et al., 2016). Global simulations are critical to capturing the viscous instability, as it requires a larger radial range than is typically treated in shearing box simulations. They also offer a better potential for comparison with observations. However, the disks simulated in the Mishra et al. (2016) paper were not true Shakura-Sunyaev disks and did not begin on the thermal equilibrium curve.

In the present work, we set up numerical experiments to directly test the thermal and viscous stability of Shakura-Sunyaev disks around Schwarzschild black holes. We do this by adding a full, relativistic Navier-Stokes treatment to the momentum and energy evolution equations in our Cosmos++ computational astrophysics code. This allows us to evolve the disk with the exact -viscosity treatment of Shakura & Sunyaev (1973). As expected the radiation-pressure-dominated cases that we consider show unmistakable evidence for thermal instability, plus another instability, that is possibly radiation-driven, but likely not the viscous instability. Because we include radiation in the simulations, we are able to directly measure the luminosities and efficiencies associated with our simulated disks, as well as characterize their variability. Each of these results are described, in turn, in Sections 3.23.4. Before getting to those results, though, we describe our numerical methods in Section 2. Finally, in Section 4, we present some discussions and conclusions. Unless otherwise noted, standard index notation is used for labeling spacetime coordinates: repeated indices represent summations, with Latin and Greek indices running over the spatial and 4-space dimensions, respectively.

## 2 Numerical Methods

All of the simulations presented in this paper use the Cosmos++ computational astrophysics code (Anninos et al., 2005; Fragile et al., 2012, 2014) to numerically evolve the equations of general relativistic radiative, viscous hydrodynamics, as described in the next section. Cosmos++ is a parallel, multi-dimensional, fully covariant, modern object-oriented (C++) radiation hydrodynamics and magnetohydrodynamics (MHD) code, written to support structured and unstructured adaptively-refined meshes, and for both Newtonian and general relativistic astrophysical applications. For this work, we utilize the High Resolution Shock Capturing (HRSC) scheme as described in Fragile et al. (2012), which shares many elements with the non-oscillatory central difference (NOCD) method presented in Anninos & Fragile (2003). In the next sections, we describe how we treat the viscosity and radiation in these simulations.

### 2.1 Relativistic, Viscous, Radiation Hydrodynamics

The contravariant stress energy tensor for a viscous fluid can be represented as a linear combination of the hydrodynamic and viscous contributions as

 Tαβ=ρhuαuβ+Pgasgαβhydro−(QB/c2)uαuβ−QBgαβ−QαβSviscous . (1)

Collecting terms, this can be written more conveniently as

 Tαβ=(ρh−QB/c2)uαuβ+(Pgas−QB)gαβ−QαβS , (2)

where is the fluid rest-mass density, is the specific enthalpy, is the speed of light, is the contravariant four-velocity, is the fluid pressure [for an ideal gas , where is the fluid internal energy density and is the adiabatic index], is the bulk (scalar) viscosity, is the (symmetric) shear tensor viscosity, and is the curvature metric. The viscous forces are explicitly split into bulk and shear viscosities for implementation convenience. We emphasize the and may include both molecular and artificial (for shock capturing) viscosity contributions.

For the radiation, we employ a covariant formulation of the closure scheme (Levermore, 1984; Sa̧dowski et al., 2013), which assumes that the radiation is isotropic in the radiation rest frame. If and represent the radiation energy density in the radiation rest frame and the spatial components of the radiation rest frame 4-velocity, respectively, then the radiation stress tensor becomes (Sa̧dowski et al., 2013; Fragile et al., 2014)

 Rαβ=43ERuαRuβR+13ERgαβ . (3)

The closure scheme has been shown to underestimate the flux coming out of plane-parallel atmospheres compared to semi-analytic solutions (Yan-fei Jiang, private communication). For our work, this means that we are most likely underestimating the radiative cooling of our disks. Since we find that our radiation-pressure-dominated disks are thermally unstable, usually with cooling dominating over heating, this suggests that our conclusions would be even stronger if we were using a more accurate radiation hydrodynamics scheme.

The eight equations of radiation hydrodynamics (energy plus three components of momentum for both the fluid and radiation fields) are derived from the conservation of fluid and radiation stress energy: and , where is the radiation 4-force density, which couples the fluid and radiation fields (Mihalas & Mihalas, 1984). In this work, we generalize the radiation four-force density from that presented in Fragile et al. (2014) to include a Compton scattering component. In this work

and we assume Kramers-type opacity laws. Since free-free absorption is the most relevant atomic absorption process for stellar mass black hole accretion disks, the appropriate Planck and Rosseland means (for solar metallicity and a hydrogen mass fraction of ) are and (e.g. Hirose et al., 2009), respectively, where is the ideal gas temperature of the fluid in Kelvin and is the density in g cm. In this work, we assume the flux mean, , is the same as the Rosseland mean and the J-mean, , is the same as the Planck mean. Because we neglect ionization/recombination processes and composition effects, we use a constant electron scattering opacity, . In this work, we assume that the electron-ion equilibration time is sufficiently short for the electrons to be at the same temperature as the ions (see Ressler et al., 2017; Sa̧dowski et al., 2017, for discussions of the validity and impact of this choice). In addition to the conservation of energy and momentum, we also require an equation for the conservation of mass: . Note that we do not presently include an equation for the conservation of photon number, though this choice is not expected to significantly affect our results (Sa̧dowski & Narayan, 2015).

#### 2.1.1 Conservation Equations

The conservation of stress energy for the fluid can be expanded and rearranged as

 ∂t(√−gT0ν)+∂i(√−gTiν)=√−gTμσ Γσμν+√−gGν . (5)

By defining the total energy density and momentum density from the full stress tensor as and , the corresponding conservation equations in this framework are

 ∂tE+∂i(−√−gTi0)=−√−gTμσ Γσμ0−√−g G0 (6)

or

 ∂tE+∂i(EVi)+∂i[√−g(P−QB−Q00)Vi+√−g Qi0]=−√−gTμσ Γσμ0−√−g G0 , (7)

for energy, and

 ∂tSj+∂i(√−gTij)=√−gTμσ Γσμj+√−g Gj (8)

or

 ∂tSj+∂i(SjVi)+∂i{√−g[(P−QB)g0j−Q0j]Vi+√−g Qij}=√−gTμσ Γσμj+√−g Gj , (9)

for momentum. Likewise, defining the radiation total energy density and radiation momentum density as and , the corresponding conservation equations are

 ∂tR+∂i(√−g Ri0)=√−g Rαβ Γβ0α−√−g G0 , (10)

and

 ∂tRj+∂i(√−g Rij)=√−g Rαβ Γβjα−√−g Gj , (11)

respectively. These equations are discretized using a second-order finite volume representation, with third-order piecewise-parabolic interpolations for the flux reconstructions, and evolved forward using a second-order, low-storage Euler time-stepping scheme.

#### 2.1.2 Primitives

At the end of each time cycle a series of coupled nonlinear equations are solved to extract the primitive fields from the evolved conserved fields, using the semi-implicit method described in Fragile et al. (2014). Essentially the method utilizes Newton-Raphson iteration to solve a linear matrix system constructed from the primitive field dependency of all of the conserved quantities. Thus within each iteration one constructs a ( for the case of radiation hydrodynamics) Jacobian matrix evaluated using initial guesses for the primitive solutions. Here is a vector list of conserved fields, while is a vector list of corresponding primitive fields. We introduce , with , as the primitive velocity in place of , where . The energy and momentum matrix elements are of particular interest since they are based on the stress energy tensor, , , and differ from ideal hydrodynamics by the addition of viscous terms. Keeping in spirit with the quasi-static approach to solving the primitive value problem, the stress tensor built from primitive fields within each Newton iteration includes viscous terms, but the row-wise Jacobian matrix elements neglect their variation as a higher order correction. Viscous terms are thus assumed constant during the primitive solve, but are updated during each hydrodynamic cycle with the same high order temporal and spatial discretization as the radiation hydrodynamics. The one exception is the treatment of the 4-velocity time derivatives, which are approximated as first order differences over each sub-cycle step.

#### 2.1.3 Viscosity Stress Tensor

Physical (molecular) viscosity takes the generic form

 Qαβ=μhσ(β∇σuα)+(μB−μ3)hαβ∇σuσ , (12)

where and are the shear and bulk coefficients, respectively, in units of mass/(lengthtime), and is the projection operator. The symmetric tensor satisfies the orthogonality condition , offering a convenient way to compute its time-time and time-space components.

Following the standard theory of thin disks, the shear viscosity coefficient is calculated as

 μ=νρ=αρcsH , (13)

where is the thermal sound speed (including both gas and radiation contributions), is the disk height, and is the Shakura-Sunyaev viscosity parameter. In this work, is assumed to be a constant, while and are calculated from the local conditions. Specifically, , where is the total pressure, and , although we also tested one case where , where is the Keplerian angular frequency. Note that this height, , is simply used to scale ; it does not restrict where viscosity is applied. In all cases, the effective height is limited to ; this was done to prevent having very large viscosity in the background gas where fluctuates considerably. We also include a bulk viscosity of the same form, , whose main effect is to help damp out spurious sound waves that come from mismatches in the initial conditions and outer boundary.

For added numerical stability and to preserve causality, we restrict the hydrodynamical timestep to be less than the minimum calculated from the local viscosity speed (). We further limit the viscous source tensor to

 Qμν=Qμνmax[1,√QαβQβα/(χQPtot)] , (14)

where is a parameter of order unity that controls the amount of limiting. This is useful for preventing the anisotropic component of viscosity from exceeding the isotropic pressure.

### 2.2 Simulation Setup

To initialize our simulations, we start from the Novikov & Thorne (1973) generalization of the Shakura-Sunyaev thin disk. As we are only considering a limited radial range, we do not require all three regions of the solution. Instead, whenever we want to initialize a gas-pressure-dominated case (appropriate for ), we only initialize the so-called “middle” region (for this , the so-called “inner” region only exists inside the innermost stable circular orbit, or ISCO). On the other hand, whenever we want to initialize a radiation-pressure-dominated case, we only initialize the “inner” region. We follow the form of the Novikov-Thorne solutions given in Abramowicz & Fragile (2013), although all we actually require are the radial dependencies of and , where is the mid-plane density. We also include the small radial velocity, , from Penna et al. (2012), associated with the slow radial drift of material through the disk.

For the vertical profile, we solve for the vertical hydrostatic equilibrium, either assuming an isothermal disk in the gas-pressure-dominated case or a polytropic EOS in the radiation-pressure-dominated case. For the former, the solution yields

 ρ(R,z)=ρ0e−z2/2H2 (15)

and

 Ptot(R,z)=GMBHH2R3ρ(R,z) . (16)

In the latter case,

 ρ(R,z)=ρ0[1−z22H2]1/(ΓNT−1) (17)

and

 Ptot(R,z)=κρΓNT , (18)

where

 κ=GMBHH2ΓNT(ΓNT−1)ρΓNT−10R3 . (19)

Assuming the gas and radiation are in local thermodynamic equilibrium for the initial, analytic solution, we partition the pressure as

where and is the radiation constant. We can now solve for . The initial azimuthal velocity is taken to be Keplerian, . Note that we neglect the additional corrections to the Novikov-Thorne solution suggested by Penna et al. (2012), but since we are just using these conditions to initialize our simulations, this should not matter too much. For the background, we initialize a cold (), low density (), static () fluid.

The temperature found above is also used to set the radiation field. In the frame of the fluid, the radiation energy density is

while the flux, , is initially set equal to the gradient of this quantity. To get the radiation density in the radiation rest frame, , and the radiation rest-frame four-velocity, , we follow the transformation procedure outlined in Sa̧dowski et al. (2013).

For simulation S01E, the disk initially extends from to the outer domain boundary, . For all other simulations, the disk is initially truncated at . We chose to do this because the analytic solution becomes extremely thin close to the ISCO and leads to very large gradients that cause numerical difficulties at the start of the simulations. Instead, we let these simulations fill in this region as they evolve. For the inner boundary, we do not go all the way to the black hole event horizon, , but rather truncate our domain at some intermediate radius, (between and ). This is done strictly to reduce computational cost. We utilize logarithmic spacing of the form

 x1=1+ln(rrBH) (22)

in the radial direction. We also only consider a small wedge in , with resolution concentrated toward the mid-plane by a function of the form

 θ=x2+12[1−q]sin(2x2) , (23)

with values of provided in Table 1. All simulations that we report in this paper use a resolution of . Although this may seem like a modest resolution, the combination of a limited range in and the concentrated latitude coordinate results in an effective resolution in the vertical direction as small as 24 m in some simulations. Lower resolution test simulations exhibit the same basic behaviors, suggesting that our conclusions are robust. The inner, top, and bottom boundaries use outflow conditions, where matter and radiation are only allowed to flow off the grid. At the outer boundary we hold the ghost zone values at those prescribed by the Novikov-Thorne solution, though material (and radiation) can still leave the grid through this boundary.

In this work, we consider five cases, which can be illustrated nicely by considering their positions in the - plane, compared to an appropriate thermal equilibrium () curve (Figure 1). All five simulations assume a Schwarzschild black hole with . The viscosity parameter is . The value of is chosen to make connection with earlier numerical work (e.g. Hirose et al., 2009; Jiang et al., 2013; Mishra et al., 2016), while the value of is motivated by values typically seen in global MHD simulations of MRI turbulence with weak, local fields (Hawley et al., 2011, 2013, and references therein). As can be seen in Figure 1, simulation S01E, with , lies on the putative stable, gas-pressure-dominated branch, while simulation S1E, with , lies near the transition between the stable and unstable branches, and S3E, S3Ep, and S10E, with , 3, and 10, respectively, lie on or near the unstable, radiation-pressure-dominated branch. Simulations S3E and S3Ep differ only in that the temperature is perturbed by a factor of 1.5 above its initial, equilibrium value in simulation S3Ep. Table 1 gives the full parameters for each simulation.

## 3 Results

The Shakura-Sunyaev -viscosity model should, in general, only be studied on timescales longer than the viscous time, ), where the effects of turbulence can, in principle, be averaged over. For thin disks with , s at , i.e., roughly of the same order as our simulations. For reference, the orbital period at the ISCO in these simulations is s, meaning all these simulations ran for hundreds of ISCO orbits. This is especially important when considering issues of stability, as simulations would need to run for at least this long to be deemed stable. In our case, many of our simulations showed clear signs of instability on considerably shorter timescales. For instance, of the simulations shown in Figure 2, only simulation S01E appears to have remained stable (i.e., close to its initial configuration) for its duration. The other simulations, S3E and S10E in this case, by contrast, show clear evidence of vertical collapse, indicative of the thermal instability. This same conclusion is even more apparent in Figure 1, where we follow the evolutionary tracks of each simulation in the - plane. While simulation S01E remains close to its initial location, most of the others drift a considerable distance in , again as expected for a thermal instability. In the next few sections, we will explore each of these cases more thoroughly.

### 3.1 Stable, Gas-Pressure-Dominated Disk

Before presenting further evidence that some of our simulations exhibit true instabilities, it may be helpful for us to demonstrate that our code can successfully evolve a stable, thin accretion disk, which was the main purpose of simulation S01E. This was the only gas-pressure-dominated (i.e., middle region) case that we considered. Figure 3 shows time histories of two of the key disk variables – and – for this simulation. The takeaway point is that, other than a brief initial adjustment, there is very little evolution in these variables, indicating stable evolution at close to the prescribed configuration.

Although such thermal stability is the expected behavior for gas-pressure-dominated disks based upon the Shakura-Sunyaev model, there is observational evidence that this particular disk configuration may not manifest itself in nature. At the low accretion rate associated with this model (), BHXRBs uniformly appear in the so-called “Hard” state (Maccarone, 2003; Remillard & McClintock, 2006), which is usually interpreted as the inner part of the accretion flow being hot, thick, and advection-dominated (Esin et al., 1997; Done et al., 2007), very unlike the geometrically thin, Shakura-Sunyaev solution studied here. We can only speculate that additional physics that is not represented in our simulation (e.g. thermal conduction) may prevent this particular solution from being realized.

### 3.2 Thermal Instability

The rest of the simulations we consider all lie on the radiation-pressure-dominated branch and are thus subject to various instabilities. We start by considering the thermal instability. To understand this instability, we note that the local cooling rate per unit area of the disk surface is proportional to the radiation energy density at the mid-plane, . Since the opacity in these disks is dominated by electron scattering, the cooling should scale as . The heating rate per unit area, on the other hand, is given by the vertically integrated stress times the rate of strain, or . In the case where the disk is supported vertically predominantly by radiation pressure, the disk half thickness is proportional to the surface radiation flux, . For a standard alpha disk, where , this means that the heating rate scales as . Figure 4 demonstrates that our simulated radiation-pressure-dominated disks indeed follow the and scalings. The result is that perturbative increases (decreases) in the midplane temperature should result in runaway heating (cooling). Furthermore, the resulting instability should grow on the thermal timescale, .

Figure 5, which shows the density-squared-weighted height,

 =  ⎷∫θmaxθmin√−gρ2(θ−π/2)2dθ∫θmaxθmin√−gρ2dθ , (24)

for simulations S3E and S10E, clearly points to a thermal collapse (follow a vertical line at any radius and notice the height decreasing with time until it stabilizes at a new value). From the figure, it appears that the collapse happens slightly faster than the local thermal time (indicated by the solid, white line), which was also seen in Mishra et al. (2016). Simulation S3Ep, the perturbed case, initially shows a modest growth in height, but still ultimately collapses on roughly the thermal timescale (left panel of Figure 6). The collapse of simulation S1E is not as dramatic as these others, as it represents an intermediate case.

Figure 7 confirms that the collapse seen in Figure 5 happens because cooling (blue shades) rather quickly exceeds heating (red shades) in the disk. The process seems to accelerate as the disks collapse. At later times, simulation S3E appears to reach a new thermal equilibrium (green shades), at least in part of the disk (). This is consistent with Figure 1, where it appears that this simulation collapses down to the lower stable (gas-pressure-dominated) branch and settles into a new solution there. Similar transitions between the unstable to stable branches were suggested in some earlier shearing box studies (Agol et al., 2001; Turner et al., 2002). Simulation S10E, on the other hand, never regains its thermal equilibrium, not quite reaching the lower branch in Figure 1. We suspect that this is simply the result of the simulation not having sufficient resolution to follow the collapse further. We expect that with additional resolution this simulation would continue to collapse until it reached the stable branch and established a new thermal equilibrium.

It is important to note the thermal collapse (i.e., runaway cooling) of all three simulations (S3E, S3Ep, and S10E), especially S3Ep. For that simulation, we purposefully perturbed the initial disk solution to a higher temperature to put it above the thermal equilibrium curve (see Figure 1). According to our scaling arguments above, this should have resulted in runaway heating and thermal expansion of the disk. Figure 6 shows that this was the case for a very brief period of time. Ultimately, however, this simulation followed an evolutionary track very similar to the unperturbed simulation S3E, although on a somewhat delayed timescale. This has been a common occurrence in numerical simulations of thermally-unstable, radiation-pressure-dominated disks (Jiang et al., 2013; Mishra et al., 2016), i.e., they tend to exhibit thermal collapse preferentially over thermal expansion with only a few exceptions (e.g. Fig. 1 of Jiang et al., 2013). In our case, one thing to consider is that our disks do not actually start in hydrostatic equilibrium, even though we initiate them as close as possible to the Shakura-Sunyaev solution. The initial setup may, in some way, favor cooling over heating to such a degree that all of the simulations end in collapse. However, the reason for this apparent preference, and even whether it is physical or numerical, remains unclear at this time.

The longer term evolution of these disks would of course be of interest. As the new solutions have lower than what is being fed in from larger radii, mass must be accumulating somewhere in the disk, as is indeed happening near the outer boundary of each simulation. This could potentially lead to some sort of limit-cycle behavior, where the disk switches between low- and high- solutions (Janiuk et al., 2002; Ohsuga, 2006). This possibility will be explored in future work.

In this paper we are only considering radiation hydrodynamic simulations, i.e., the role of magnetic fields is neglected, other than that they are presumed to be the underlying source of the stresses that we model as a viscosity. However, if strong magnetic fields are present, they may, in principle, stabilize radiation-pressure-dominated disks against the thermal instability seen in our work (Begelman & Pringle, 2007; Oda et al., 2009; Sa̧dowski, 2016). Even with weaker magnetic fields, the cooling of the disk may be altered by magnetic buoyancy.

### 3.3 Viscous Instability

In the radiation-pressure-dominated regime, the vertically integrated stress, , in a standard Shakura-Sunyaev disk can be shown to be inversely proportional to the surface density, . This means that the evolution equation for , which amounts to a nonlinear diffusion equation,

 ∂Σ∂t∝∂∂r[r2WRϕ] , (25)

has a negative effective diffusion coefficient (Lightman & Eardley, 1974). This implies that the disk should break up into high and low surface density rings with on a timescale , which turns out to be roughly the same as the thermal timescale, . We indeed see a breaking of our disks into rings, as exhibited by the long, vertical streaks of surface density in the spacetime diagrams of Figure 8.

Although the disks in both simulations 3E and 10E break up into distinct radial rings with on roughly the local thermal timescale (solid line), as predicted by Lightman & Eardley (1974), we do not believe this is the same instability identified in that paper. Our main reason for arguing this is that, as shown in Figure 9, the vertically integrated stress, , is actually high in regions where is high and low where is low, exactly opposite what is predicted by the viscous instability.

The first clue as to what may, in fact, be driving our simulated disks to break into rings is to note that it only happens in the radiation-pressure-dominated cases (compare for simulation S01E in Figure 3 to for simulations S3E and S10E in Figure 8). In fact, the low-density gaps in the disk appear to form where there are local peaks in the ratio of radiation-to-gas pressure, as shown in Figure 10. Also consistent with this idea is the fact that the opacity is notably lower in the low-density gaps, as shown in Figure 11, providing the radiation an easier path to escape from the disk. This is expected behavior for a radiation-pressure-dominated medium; once low density channels form, the radiation will naturally follow these channels and thus reinforce the perturbations. The only question is, what allows the initial density inhomogeneities to form? In this case, it appears it may be convection. We find that our radiation-pressure-dominated disks are convectively unstable according to the Schwarzschild condition

 dTdz<(1−1Γ)TPdPdz , (26)

whereas our gas-pressure-dominated case (S01E) is stable by the same criterion.

### 3.4 Luminosity & Variability

Because these are global, radiation hydrodynamic simulations, we can directly extract light curves and compare the measured luminosity with predictions for thin accretion disks. First, let us look at the measured accretion rates of the simulations, as it is mass accretion that ultimately powers the luminosity. Not surprisingly, given the instabilities we have already identified, the measured mass accretion rates often do not match the input value for that simulation. This is because, as the thermally unstable solutions collapse and seek new, stable equilibria, the mass accretion rate must necessarily drop. This follows from the fact that, for a given , the stable, gas-pressure-dominated branch has a lower than the unstable, radiation-pressure-dominated one. Figure 12 shows the values measured near the inner boundaries of each simulation domain, where

 ˙M=−2π∫√−gρurdθ . (27)

For clarity, we only show the mass accretion histories of three of our simulations: S01E, S1E, and S10E. Note that simulation S01E accretes at exactly its targeted accretion rate (), as expected for a stable, gas-pressure-dominated case, whereas simulations S1E and S10E accrete at significantly below their target rates, consistent with their thermal collapse. Simulations S3E and S3Ep (not shown) also accrete at significantly below their target levels. The thermally unstable simulations also exhibit far greater variability in their mass accretion rates than the stable S01E case.

Next we consider how efficiently this matter accretion is converted into radiation. To calculate the luminosity, we integrate

 L=−∫S√−gRθtdAθ , (28)

where is the radiative flux in the -direction, normal to an area element . For this work, we take the surface to be at constant near the top and bottom boundaries of each grid, out to a radius of . We are obviously underestimating the total luminosity as we are ignoring any radiation originating outside , as well as any radiation passing through the inner and outer radial boundaries of our defined volume. However, both contributions should be relatively minor in terms of their contribution to the luminosity that would be measured at infinity. Indeed, we find that the radiated luminosity (Figure 13, left panel) closely matches the theoretical expectations corresponding to the observed mass accretion rates in Figure 12. For a Schwarzschild (non-rotating) black hole, the radiative efficiency, , should be approximately 6%. In the right panel of Figure 13, we plot the measured radiative efficiencies from some of our simulations as a function of time. Clearly, for simulations S1E and S10E, the efficiency varies considerably, but for the simulations that do find a stable solution, the value of is generally within a factor of a couple of its predicted value. Keep in mind, too, that because is a measure of the instantaneous mass accretion rate at the inner boundary while the luminosity is measured over a much larger region, the two are not expected to vary on the same timescale, so should be expected to deviate significantly on short timescales.

With light curves in hand, we can now look at their behavior in frequency space to explore the nature of any variability. As we dump data three times per ISCO orbital period (i.e., every 1 ms), we have a Nyquist frequency of 500 Hz. Figure 14 gives the Fourier power of the full light curves for simulations S01E and S10E (not the smoothed light curves used for display purposes in Figure 13). The full power spectra for simulations S01E (thick line, left panel) and S10E (right panel) both appear to have an underlying broken power law dependence on , although of differing slopes before and after the break in each case. Interestingly, the break occurs around 100 Hz in both cases. While intriguing, it is not obvious why this frequency, which corresponds to about 3 orbital periods at the ISCO, should be special. None of the power spectra show significant evidence of excess power near certain discrete frequencies, i.e. there is no evidence for quasi-periodic oscillations (QPOs) (Remillard & McClintock, 2006).

Overall, the rms variability in these power spectra are notably high. Remember that in the soft state, BHXRBs often exhibit rms variability % (van der Klis, 2006; Done et al., 2007). For simulation S10E, the high variability is obvious in the lightcurve in Figure 13, and is intimately tied to the instability of that disk. However, even the stable, gas-pressure-dominated simulation, S01E, has a surprisingly high fractional rms variability. It seems, though, that most of the power in the S01E spectrum actually comes from trying to fit the luminosity spikes that occur around , 38700, 51900, 73400, and in Figure 13. If, instead, we only analyze the rms variability during a “quiet” interval of the S01E lightcurve, say from 10000 to , then we find a much weaker power density spectrum, as shown by the thin, red, dashed line in the left panel of Figure 14. Those spikes in luminosity come from spikes in the mass accretion rate at the inner edge of the disk (note the corresponding spikes in Figure 12). During these spikes, it seems the disk temporarily fills in the region inside the ISCO. Figure 15 shows an image of the disk gas density during a normal, quiescent phase (left panel), where the disk is truncated near , and during one of the accretion spikes (right panel), where the disk extends more or less to the inner edge of the simulation domain ( in this case).

The underlying source of true variability in all these simulations seems to be a spectrum of sound waves and diskoseismic modes that are causing fluctuations in the surface density. These modes will be the subject of a future paper.

## 4 Discussion & Conclusions

In this paper we presented results from a set of numerical simulations of thin-disk accretion onto a Schwarzschild black hole that were designed to be directly comparable to the Shakura-Sunyaev disk model, with particular focus on the unstable, radiation-pressure-dominated branch. As predicted by theory and shown in earlier shearing box simulations, all disks that started on the unstable branch exhibited runaway behavior on timescales comparable to the local thermal time.

Somewhat surprisingly, though, all of the thermally unstable disks underwent runaway cooling and collapse. We did not see any examples of sustained runaway heating and expansion, even in the case of our simulation S3Ep, for which the disk temperature was intentionally perturbed upward by a factor of 1.5 to try to trigger runaway heating. This tendency of thermally unstable disks to collapse, rather than expand, was also noted earlier in shearing box simulations (Jiang et al., 2013). This could indicate that there is something we are still not understanding about this instability.

We also found that the thermally unstable disks broke up into rings of high and low surface density. Although this behavior is consistent with what is predicted for the viscous instability (Lightman & Eardley, 1974), we argued that it is more likely a radiation pressure instability. Whereas the Lightman-Eardley instability predicts an inverse correlation between surface density and vertically integrated stress, we found a positive correlation. Instead, we argued that the low density “gaps,” triggered initially by a convective instability, provide preferred channels for radiation to escape, which reinforces the evacuation of these regions. We speculate that in three dimensions, this instability may break the disk up into inhomogeneous clumps of high and low surface density.

Because radiation was included self-consistently in our simulations, we were able to calculate lightcurves and power density spectra corresponding to each simulation. The lightcurves, when coupled with the measured mass accretion rates, revealed radiative efficiencies close to the expected value of 6% for a Schwarzschild black hole. The rms variability, on the other hand, was seen to be significantly higher in our numerical simulations (up to order unity) than in observed systems (). In most cases, this variability was directly related to the instability of the disks. For simulation S01E, the variability was limited to very brief episodes where the disk experienced sharp spikes in mass accretion rate and luminosity. These spikes will be explored further in future work.

Although our results are broadly consistent with theoretical predictions and previous numerical work, they appear to be at odds with observations of BHXRBs. Our unstable simulations – S1E, S3E, S3Ep, and S10E – span the exact luminosity range where disks in nature appear to be most stable (e.g. van der Klis, 2006; Done et al., 2007). Furthermore, our one stable simulation – S01E – is at a mass accretion rate where BHXRBs appear to be in a spectrally hard state, inconsistent with a simple blackbody disk spectrum (Maccarone, 2003; Remillard & McClintock, 2006). In both cases, therefore, it seems there is still important physics that is not being properly captured by standard theory or numerical simulations. Strong magnetic fields have been proposed as a potential solution to the first problem, as such fields may help stabilize radiation-pressure-dominated disks (Begelman & Pringle, 2007; Oda et al., 2009; Sa̧dowski, 2016). Thermal conduction has been suggested as a possible solution to the second problem, by providing a mechanism by which a cold, thin disk may evaporate into a hot, thick flow at low accretion rates (Mayer & Pringle, 2007). Each of these will be considered in future work.

The authors would like to thank Omer Blaes and Chris Done for useful discussions and feedback regarding this work. This project was supported by National Science Foundation grants AST-1211230, AST-1616185, and PHY-1125915. It used resources from the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575 and PROMETHEUS supercomputer in the PL-Grid infrastructure in Poland. Additionally, SME acknowledges support from the College of Charleston Undergraduate Research and Creative Activities Board, through SURF grants 2015-10 and 2016-9. Work by PA was performed in part under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. BM and WK acknowledge support from Polish NCN grants 2013/08/A/ST9/00795 and 2014/15/N/ST9/04633. \softwareCosmos++ (Anninos et al., 2005; Fragile et al., 2012, 2014)

## References

• Abramowicz & Fragile (2013) Abramowicz, M. A., & Fragile, P. C. 2013, Living Reviews in Relativity, 16, 1
• Agol et al. (2001) Agol, E., Krolik, J., Turner, N. J., & Stone, J. M. 2001, ApJ, 558, 543
• Anninos & Fragile (2003) Anninos, P., & Fragile, P. C. 2003, ApJS, 144, 243
• Anninos et al. (2005) Anninos, P., Fragile, P. C., & Salmonson, J. D. 2005, ApJ, 635, 723
• Balbus & Hawley (1992) Balbus, S. A., & Hawley, J. F. 1992, ApJ, 400, 610
• Balbus & Hawley (1998) —. 1998, Reviews of Modern Physics, 70, 1
• Begelman & Pringle (2007) Begelman, M. C., & Pringle, J. E. 2007, MNRAS, 375, 1070
• Done et al. (2007) Done, C., Gierliński, M., & Kubota, A. 2007, A&A Rev., 15, 1
• Esin et al. (1997) Esin, A. A., McClintock, J. E., & Narayan, R. 1997, ApJ, 489, 865
• Fragile et al. (2012) Fragile, P. C., Gillespie, A., Monahan, T., Rodriguez, M., & Anninos, P. 2012, ApJS, 201, 9
• Fragile et al. (2014) Fragile, P. C., Olejar, A., & Anninos, P. 2014, ApJ, 796, 22
• Hawley et al. (2011) Hawley, J. F., Guan, X., & Krolik, J. H. 2011, ApJ, 738, 84
• Hawley et al. (2013) Hawley, J. F., Richers, S. A., Guan, X., & Krolik, J. H. 2013, ApJ, 772, 102
• Hirose et al. (2009) Hirose, S., Krolik, J. H., & Blaes, O. 2009, ApJ, 691, 16
• Honma et al. (1991) Honma, F., Kato, S., & Matsumoto, R. 1991, PASJ, 43, 147
• Honma et al. (1992) Honma, F., Matsumoto, R., & Kato, S. 1992, PASJ, 44, 529
• Janiuk et al. (2002) Janiuk, A., Czerny, B., & Siemiginowska, A. 2002, ApJ, 576, 908
• Jiang et al. (2013) Jiang, Y.-F., Stone, J. M., & Davis, S. W. 2013, ApJ, 778, 65
• Krolik et al. (2007) Krolik, J. H., Hirose, S., & Blaes, O. 2007, ApJ, 664, 1045
• Levermore (1984) Levermore, C. D. 1984, JQSRT, 31, 149
• Lightman & Eardley (1974) Lightman, A. P., & Eardley, D. M. 1974, ApJ, 187, L1
• Maccarone (2003) Maccarone, T. J. 2003, A&A, 409, 697
• Mayer & Pringle (2007) Mayer, M., & Pringle, J. E. 2007, MNRAS, 376, 435
• McKinney et al. (2014) McKinney, J. C., Tchekhovskoy, A., Sadowski, A., & Narayan, R. 2014, MNRAS, 441, 3177
• Mihalas & Mihalas (1984) Mihalas, D., & Mihalas, B. W. 1984, Foundations of radiation hydrodynamics
• Mishra et al. (2016) Mishra, B., Fragile, P. C., Johnson, L. C., & Kluźniak, W. 2016, MNRAS, 463, 3437
• Novikov & Thorne (1973) Novikov, I. D., & Thorne, K. S. 1973, in Black Holes (Les Astres Occlus), ed. C. Dewitt & B. S. Dewitt, 343–450
• Oda et al. (2009) Oda, H., Machida, M., Nakamura, K. E., & Matsumoto, R. 2009, ApJ, 697, 16
• Ohsuga (2006) Ohsuga, K. 2006, ApJ, 640, 923
• Ohsuga et al. (2005) Ohsuga, K., Mori, M., Nakamoto, T., & Mineshige, S. 2005, ApJ, 628, 368
• Penna et al. (2012) Penna, R. F., Sa̧owski, A., & McKinney, J. C. 2012, MNRAS, 420, 684
• Remillard & McClintock (2006) Remillard, R. A., & McClintock, J. E. 2006, ARA&A, 44, 49
• Ressler et al. (2017) Ressler, S. M., Tchekhovskoy, A., Quataert, E., & Gammie, C. F. 2017, MNRAS, 467, 3604
• Ross et al. (2017) Ross, J., Latter, H., & Tehranchi, M. 2017, ArXiv e-prints, arXiv:1703.00211
• Sa̧dowski (2016) Sa̧dowski, A. 2016, MNRAS, 459, 4397
• Sa̧dowski & Narayan (2015) Sa̧dowski, A., & Narayan, R. 2015, MNRAS, 454, 2372
• Sa̧dowski et al. (2013) Sa̧dowski, A., Narayan, R., Tchekhovskoy, A., & Zhu, Y. 2013, MNRAS, 429, 3533
• Sa̧dowski et al. (2017) Sa̧dowski, A., Wielgus, M., Narayan, R., et al. 2017, MNRAS, 466, 705
• Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
• Shakura & Sunyaev (1976) —. 1976, MNRAS, 175, 613
• Shaviv (2001) Shaviv, N. J. 2001, ApJ, 549, 1093
• Turner et al. (2002) Turner, N. J., Stone, J. M., & Sano, T. 2002, ApJ, 566, 148
• van der Klis (2006) van der Klis, M. 2006, Rapid X-ray variability, Compact stellar X-ray sources, edited by Lewin, W.H.G. and van der Klis, M., pages 39-98, Cambridge, UK: Cambridge University Press
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters