S-curve slope for generalized opacities

March 18, 2020, Revision: 1.37

1 Introduction

In the standard Shakura & Sunyaev (1973) accretion disk solution, the mass accretion rate Ṁ increases with surface density Σ. Since the early 1980s, it is well known that a monotonous increase of Ṁ with Σ implies viscous stability of the disk, while a negative dependence where the slope

s ≡ d ln Ṁ/d ln Σ  (1)

implies instability (Bath & Pringle, 1981, 1982; Meyer & Meyer-Hofmeister, 1981, 1982; Cannizzo et al., 1982). Shakura & Sunyaev (1973) found s = 10/7 for moderately hot disks in which the opacity is governed by bound-free transitions, which is well approximated by the Kramers opacity. For hotter disks – when the dominant opacity is given by electron scattering, the slope s is still positive, but the value is now ??. The slope also depends on the equation of state. For moderately hot disks, the perfect gas equation is appropriate, but in hotter disks radiation pressure becomes dominant.

Particularly interesting are cases where the dependence Ṁ(Σ) can be my multivalued, i.e. different values of Ṁ are possible for the same value of Σ. This is relevant for explaining so-called high and low accretion states. In such a case the dependence of Ṁ on Σ is S shaped. The different branches are believed to be associated with the onset of convection (ref?).

A major uncertainty in modelling accretion discs is the nature of turbulence. Since the 1990s, the source of turbulence is generally believed to be the magnetorotational instability (Balbus & Hawley, 1991, 1992; Hawley & Balbus, 1991, 1992). Here, a moderately strong magnetic field leads to the development of turbulence (Hawley et al., 1995). The magnetic field can either be external to the disk, e.g., the field from the central object, or it could be generated by an internal dynamo within the disk. In the latter case, the dynamo is driven entirely from MRI-driven turbulence. This requires a doublypositive feedback between the magneto-rotational and dynamo instabilities. Local numerical simulations using the shearing box approximation have shown that this mechanism can lead to sustained turbulence with a finite accretions stress (sum of Reynolds and Maxwell stresses) of the appropriate sign for being able to tap energy from the differential rotation (Brandenburg et al., 1995; Hawley et al., 1996; Stone et al., 1996). This energy feeds a turbulent cascade toward progressively smaller length scales where eventually microphysical viscous and ohmic dissipation convert this energy into heat. The resulting heating is usually described in terms of a turbulent viscosity and would thus be proportional to the square of the shear rate. The turbulent viscosity is assumed to be proportional to a certain fraction of the local sound speed and the local pressure scale height. This assumption is plausible, but details are not borne out by numerical simulations, which instead suggest significant corrections to the assumed quadratic dependence of the shear rate (Abramowicz et al., 1996; Hawley et al., 1999; Pessah et al., 2007). Also the vertical dependence of the heating rate follows a more complicated profile than what is expected from the simple parameterization in terms of a turbulent viscosity. For these reason it is crucial to perform numerical simulations in suitable settings with radiation transport included.

In a local model of an accretion disk, owing to the use of shearing-periodic boundaries, no net mass enters or leaves the computational domain. Thus, the mean density and therefore also Σ is on average constant. In such a simulation, we can therefore not expect an evolution of Σ. The shearingperiodic boundary conditions also do not permit an accretion flow, because left and right sides of the box are equivalent, and there is no way of telling whether the accreting body would lie to the left or the right of the shearing box. However, there can still be a turbulent stress and there could also be corresponding disk solutions that could in principle jump between different states.

The goal of the present work is to explore the nature of high and low accretion state solutions in a parameter regime that is relevant to cataclysmic variables. Given the importance of the opacity for radiation transport and to be able to identify the dominant physical processes, it is useful to keep certain properties of the opacity adjustable. We do this by using a combination of Kramers-like opacities with adjustable exponents. This has been useful in previous studies of stellar atmospheres, where it turned out that, at least in the gray approximation, one naturally obtains polytropic layers when the exponents of of the Kramers opacity are held had fixed. To keep the model sufficiently simple, we also ignore here variable ionization.

We begin by constructing one-dimensional hydrostatic equilibrium solutions. These can be obtained by simple integration from the surface inwards. However, it is at least equally insightful to obtain such solutions using the same timestepping method as the one used for obtaining MRI solutions in the shearing box approximation.

2 Governing equations

We solve the hydrodynamic equations for logarithmic density ln ρ, velocity u, and specific entropy s, in the form

\frac{D ln ρ}{Dt} = −∇ · u,
ρ\frac{Du}{Dt} = −∇P − 2Ω × ρu + ρg_{eff} + ∇ · (2νρS),
ρT \frac{DS}{Dt} = −∇ · F_{rad} + 2νρS^{2} + H,

where P is the gas pressure, g is the gravitational acceleration, ν is the viscosity, S = \frac{1}{2}[∇u + (∇u)^{T}] − \frac{1}{3}I∇ · u is the traceless rate-of-strain tensor, I is the unit tensor, T is the temperature, and F_{rad} is the radiative flux. For the equation of state, we assume a perfect gas with p = (R/μ)T ρ, where R is the universal gas constant and μ is the mean molecular weight. The pressure is related to ρ and S via P = ρ^{γ} exp(S/c_{v}), where the adiabatic index γ = c_{p}/c_{v} is the ratio of specific heats at constant pressure and constant volume, respectively, and c_{p} − c_{v} = R/μ. To obtain the radiative flux, we adopt the gray approximation, ignore scattering, and assume that the source function S (not to be confused with the rate-of-strain tensor S) is given by the frequency-integrated Planck function, so S = (σ_{SB}/π)T^{4}, where σ_{SB} is the Stefan– Boltzmann constant. The negative divergence of the radiative flux is then given by

−∇ · F_{rad} = κρ ∮_{4π}(I − S) dΩ,  (2)

where κ is the opacity per unit mass (assumed independent of frequency) and I(x, t, n̂) is the frequency-integrated specific intensity corresponding to the energy that is carried by radiation per unit area and per unit time in the direction n̂ through a solid angle dΩ. We obtain I(x, t, n̂) by solving the radiative transfer equation,

−n̂ · ∇I = κρ (I − S),  (3)

along a set of rays in different directions n̂ using the method of long characteristics.

To explore different solutions, it is useful to be able to control the value and functional form of the opacity. We therefore choose a Kramers-type opacity given by

κ = κ_{0}ρ^{a}T^{b},  (4)

where a and b are considered as free parameters that characterize the relevant radiative processes. We compute the radiative conductivity K(ρ, T), which is given by

K(ρ, T) = \frac{16σ_{SB}T^{3}}{3κρ} = \frac{16σ_{SB}T^{3−b}}{3κ_{0}ρ^{a+1}} .  (5)

We note that, in a plane-parallel polytropic atmosphere, T (z) varies linearly with height z and in the stationary state, K(ρ, T) is constant in the optically thick part. This implies that ρ is proportional to T^{n}, where

n = \frac{3 − b}{1 + a}  (6)

is the polytropic index (not to be confused with the direction of a ray n̂).

3 Comparison with vertically integrated disks

Use

\frac{4σ_{SB}T_{c}^{4−b}}{3κ_{0}Σ^{1+a}} H = \frac{1}{2}ν_{t}Σ(\frac{3}{2}Ω)^{2}  (7)

or

\frac{4σ_{SB}(R/μ)^{−(4−b)}c_{s}^{2(4−b)}}{3κ_{0}Σ^{1+a}} H = \frac{1}{2}ν_{t}Σ(\frac{3}{2}Ω)^{2} (8)

which leads to

c_{s}^{2(4−b)}H = Kν_{t}Σ^{2+a}Ω^{2}  (9)

where

K = \frac{27κ_{0}}{32σ_{SB}} (\frac{R}{μ})^{4−b}  (10)

Eliminate sound speed and ν_{t}, so

Ω^{2(4−b)}H^{2(4−b)}H = Kα_{SS}ΩH^{2}Σ^{2+a}Ω^{2} (11)

collect

Ω^{5−2b}H^{7−2b} = Kα_{SS}Σ^{2+a}  (12)

Combine with angular momentum conservation ν_{t}Σ = ṁ, so

α_{SS}ΩH^{2}Σ = ṁ ≡ Ṁ/3π  (13)

to eliminate Σ, so we find

Ω^{5−2b}H^{7−2b} = Kα_{SS}(ṁ/α_{SS}ΩH^{2})^{2+a}  (14)

or

Ω^{2+a}H^{2(2+a)}Ω^{5−2b}H^{7−2b} = Kṁ^{2+a}α_{SS}^{−(1+a)} (15)

so

Ω^{7+a−2b}H^{11+2a−2b} = Kṁ^{2+a}α_{SS}^{−(1+a)}  (16)

Using Ω(r) = r^{3/2}/\sqrt{GM}, we have H = H(r), given by

H^{11+2a−2b} = Kṁ^{2+a}α_{SS}^{−(1+a)}Ω^{−(7+a−2b)} (17)

Alternatively, using Eq. (13), we can determine Σ = Σ(r) as

(ṁ/α_{SS}ΩΣ)^{(11+2a−2b)/2} = Kṁ^{2+a}α_{SS}^{−(1+a)}Ω^{−(7+a−2b)}

(18)

so

ṁ^{7/2−b} = Kα_{SS}^{9/2−b}Ω^{−3/2+b}Σ^{11/2+a−b}  (19)

or

Σ^{11/2+a−b} = ṁ^{7/2−b}K^{−1}α_{SS}^{−9/2+b}Ω^{3/2−b} (20)

The important relation is ν_{t} versus Σ, so now we find

ṁ ∝ Σ^{(11+2a−2b)/(7−2b)}  (21)

For Kramers’ opacity, a = 1 and b = −7/2, the exponent is (11 + 2a − 2b)/(7 − 2b) = 20/14 = 10/7, which is the familiar result. The exponent becomes negative when 11 + 2a − 2b < 0, i.e.,

b > 11/2 + a = 13/2 = 6.5 (for a = 1) (22)

which is just a bit less than the value for H^{−} opacity, or

b > 6 (for a = 1/2),  (23)

which is more appropriate for H^{−} opacity.

Figure 1: Slopes (11+2a−2b)/(7−2b) as a function of a and b (upper panel) and polytropic index n = (3 − b)/(1 + a) (lower panel).

Table 1:

variable code units  cgs units
Σ g cm^{−3} Mm 10^{8} g cm^{−2}
Ṁ  g cm^{−3} Mm^{2} km s^{−1} 10^{21} g s^{−1}
≈ 1.5^{−5}M_{⊙} yr^{−1}

4 Simulation strategy

In a local model of an accretion disk, the domain is (shearing-) periodic in the horizontal directions and the boundaries in the vertical direction are closed, so the mass in the domain is conserved and equal to the value given by the initial condition.

relative hydro stress:

\overline{ρu_{x}u_{y}}/\sqrt{\overline{ρu_{x}^{2}} \overline{ρu_{y}^{2}}}  (24)

is around −0.34 Alpha:

α_{SS}^{K} = \overline{ρu_{x}u_{y}}/\frac{3}{2} \overline{ρ}Ω^{2}H^{2}  (25)
α_{SS}^{M} = −\overline{B_{x}B_{y}/μ_{0}}/\frac{3}{2} \overline{ρ}Ω^{2}H^{2}  (26)

or

α_{SS}^{K} = \overline{u_{x}u_{y}}/\frac{3}{2}Ω^{2}H^{2}  (27)
α_{SS}^{M} = −\overline{B_{x}B_{y}/ρμ_{0}}/\frac{3}{2} \overline{ρ}Ω^{2}H^{2}  (28)

Figure 2: Analytic radiative equilibrium solutions

5 Hydrostatic, non-turbulent disk structure

\frac{dF}{dz} = H(z)  (29)
\frac{dT}{dz} = −F/K  (30)
\frac{dP}{dz} = −ρg  (31)

where g = Ω^{2}z is the local gravitational acceleration.

Assume H = const for z ≤ L_{heat} and 0 otherwise, so F = H min(z, L_{heat}). Then

\frac{dT}{dz} = −H min(z, L_{heat})/K  (32)

\frac{dP}{dz} = −ρΩ^{2}z  (33)

Divide by each other

\frac{dT}{dP} = \frac{H}{KρΩ^{2}} min(L_{heat}/z, 1)  (34)

\frac{d ln T}{d ln P} = \frac{HP}{KρT Ω^{2}} min(L_{heat}/z, 1)  (35)

Use P = RT ρ/μ, and define

∇_{rad} = \frac{HR/μ}{KΩ^{2}} min(L_{heat}/z, 1).  (36)

So we would need to solve d ln T/d ln P = ∇_{rad}, but with convection we have mixing, so the actual temperature gradient becomes equal to the adiabatic one and therefore

\frac{d ln T}{d ln P} = min(∇_{rad}, ∇_{ad})  (37)

and ∇_{ad} = 1 − 1/γ.

In Fig. 14 we show the magnetic field at the periphery of the computational domain. Note that the magnetic field is strongest some distance away from the midplane where the density is lower. While this results agrees with earlier findings (?) the magnetic field is here unrealistically weak near the midplane, This could be an artifact having used moderately low resolution and correspondingly larger overall magnetic diffusion.

6 Disk parameters

H = c_{s}^{max}/Ω  (38)
Ṁ = 3π \frac{2H_{0}L_{H}}{(\frac{3}{2}Ω)^{2}}  (39)

factor 2?

α_{SS} = Ṁ/3π ΣΩH^{2}  (40)
σ_{SB}T_{eff}^{4} = (Ṁ/3π) (\frac{3}{2}Ω)^{2}  (41)

7 Comparison with time integration

In Figs. 6 and 7 we verify solutions with Pencil 1D. ΔP = Δz^{2}/2H_{P}^{2} and H_{P} = \sqrt{RT/μ}/Ω.

8 Problems with vertical disk integration

−\frac{16σ_{SB}T^{3−b}}{3κ_{0}ρ^{1+a}} \frac{dT}{dz} = ∫ H dz  (42)

\frac{4σ_{SB}T^{4−b}}{3κ_{0}Σ^{1+a}} H^{a} = ν_{t}Σ(\frac{3}{2}Ω)^{2}  (43)

Figure 5: Comparison of solutions for Σ = 3.6 × 10^{−7}. κ_{0} = 2 × 10^{12}, with Kramers/H^{−} opacity using a_{H^{−}} = 1, b_{H^{−}} = 4, and T_{0} = 2 × 10^{4} K, ρ_{0} = 10^{−5} g cm^{−3}. Note that the models on the upper and middle branches are convectively unstable (topheavy and negative vertical entropy gradient).

9 Mixing-length treatment

A relatively simple treatment of convection is to set ∇ equal to ∇_{ad} when ∇_{rad} > ∇_{ad} and equal to ∇_{rad} otherwise. A shortcoming of this approach is that the energy flux balance is ignored: the radiative flux is now too low and now expression for Figure 6: Comparison between analytic solution (dashed, black) and Pencil Code solution (solid, red) for the H state solution with H_{0} = 1.4 × 10^{−5}.

Figure 7: Same as Fig. 6, but for the M solution with H_{0} = 6 × 10^{−6}. Note that the parameters Σ and α_{SS} are the same in both cases.

the convective energy flux is adopted. In standard mixing length theory (MLT), the convective energy flux is proportional to the superadiabatic gradient, ∇−∇_{ad}, which would obviously be zero in the primitive approach of setting ∇ = ∇_{ad} in the convection zone. Instead, we now set the total flux equal to the sum of radiative and convective energy fluxes.

Figure 8: Comparison of different ∇ as a function of z for H_{0} = 10^{−3}, Ṁ = 5.8 × 10^{−5}, and Σ = 6.3 × 10^{−6}.

Figure 9: Ratio ∇/∇_{rad}, which corresponds to the ratio of radiative to total flux when ∇/∇_{rad} < 0. (back) and ratio of convective to total flux (red) in the improved mixing length model.

The kinetic energy flux is usually neglected in standard MLT, but simulations of solar convection have shown that it is about 30% of the total flux in the bulk of the convection zone. Since both the enthalpy flux and the kinetic energy flux scale in the same way, it is possible to subsume both into a common expression for the convective energy flux such that

F_{tot} = F_{rad} + F_{conv}.  (44)

The radiative flux is still given by Eq. (??) and ∇ is still the actual logarithmic temperature gradient, but now it is no longer equal to ∇_{rad}. Instead, ...

• what is the KH time scale for given κ_{0}? How Figure 10: Growth of Alfvén speed in a turbulent MRI simulation. (Different solutions are patched together.)

Figure 11: Energy dissipation (magnetic red, and kinetic blue) as a function of t. No external heating has been applied.

does it depend on ρ?

• How is T_{mid} related to Σ?

• Discuss strategy of fixing H.

• Plot T vs τ.

References

Abramowicz, M. A., Brandenburg, A., & Lasota, J.-P. 1996, MNRAS, 281, L21

Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214

Balbus, S. A. & Hawley, J. F. 1992, ApJ, 400, 610 Figure 12: Energy dissipation (magnetic red, and kinetic blue) as a function of z.

Figure 13: Velocity and entropy fluctuation in a vertical plane. The τ = 1 surface is shown as a dashed line.

Bath, G. T., & Pringle, J. E. 1981, MNRAS, 194, 967

Bath, G. T., & Pringle, J. E. 1982, MNRAS, 199, 267

Brandenburg, A., Nordlund, Å., Stein, R. F., & Torkelsson, U. 1995, ApJ, 446, 741

Cannizzo, J. K., Ghosh, P., & Wheeler, J. C. 1982, ApJ, 260, L83

Cannizzo, J. K. 1992, ApJ, 385, 94 Figure 14: Visualization of B_{z} on the periphery of the computational domain.

Figure 15: Mean temperature versus t and z.

Hawley, J. F. & Balbus, S. A. 1991, ApJ, 376, 223

Hawley, J. F., Balbus, S. A. 1992, ApJ, 400, 595

Hawley, J. F., Balbus, S. A., & Winters, W. F. 1999, ApJ, 518, 394

Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742

Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1996, ApJ, 464, 690 Figure 16: Total, magnetic (red), and kinetic (blue) stresses.

Hirose, S., Blaes, O., Krolik, J. H., Coleman, M. S. B., & Sano, T. 2014, ApJ, 787, 1

Meyer, F. & Meyer-Hofmeister, E. 1981, A&A, 104, L10

Meyer, F. & Meyer-Hofmeister, E. 1982, A&A, 106, 34

Pessah, M. E., Chan, C.-k., & Psaltis, D. 2007, ApJ, 668, L51

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Shaviv, G., & Wehrse, R. 1991, A&A, 251, 117

Stone, J. M., Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1996, ApJ, 463, 656

$Header: /var/cvs/brandenb/tex/helena/SSdisk/notes.tex,v 1.37 2020/03/18 18:01:56 brandenb Exp $

Figure 3: S curve for α_{SS} = 5 × 10^{−3} (red), 2 × 10^{−3} (black), and 10^{−1} (blue). The slope for the unstable branch is close to −5, while those for the stable branches are close to 10/7. The red vertical line denotes Σ = 3 × 10^{−7} g cm^{−3} Mm = 0.3 g cm^{−2}.

9 Figure 4: S curve for α_{SS} = 5 × 10^{−3} (red), 2 × 10^{−3} (black), and with convection and α_{SS} = 5 × 10^{−3} (blue), The dash-dotted lines indicate the slopes 10/7 for the two stable branches without convection, 5/2 with convection, ad −5 for the unstable branch.