% $Id: notes.tex,v 1.37 2020/03/18 18:01:56 brandenb Exp $
%\documentclass{article}
\documentclass[twocolumn]{article}
\setlength{\textwidth}{160mm}
\setlength{\oddsidemargin}{-0mm}
\setlength{\textheight}{220mm}
\setlength{\topmargin}{-8mm}
\usepackage{graphicx,natbib}
\usepackage{bm,color}
\graphicspath{{./fig/}{./png/}}
\input macros
\title{S-curve slope for generalized opacities}
\author{}
\date{\today,~ $ $Revision: 1.37 $ $}
\begin{document}
\maketitle

\section{Introduction}

In the standard \cite{SS73} accretion disk solution, the mass
accretion rate $\dot{M}$ {\em increases} with surface density $\Sigma$.
Since the early 1980s, it is well known that a monotonous increase
of $\dot{M}$ with $\Sigma$ implies viscous stability of the disk,
while a negative dependence where the slope
\EQ
s\equiv\dd\ln\dot{M}/\dd\ln\Sigma
\EN
implies instability \citep{BP81,BP82,MMH81,MMH82,Cannizzo82}.
\cite{SS73} 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 $\dot{M}(\Sigma)$
can be my multivalued, i.e. different values of $\dot{M}$ are possible
for the same value of $\Sigma$.
This is relevant for explaining so-called high and low accretion states.
In such a case the dependence of $\dot{M}$ on $\Sigma$ 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 \citep{BH91,BH92,HB91,HB92}.
Here, a moderately strong magnetic field leads to the development of
turbulence \citep{HGB95}.
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 doubly-positive 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
\citep{BNST95,HGB96,SHGB96}.
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 \citep{ABL96,HBW96,PCP07}.
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 $\Sigma$ is on average constant.
In such a simulation, we can therefore not expect an evolution of $\Sigma$.
The shearing-periodic 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.

\section{Governing equations}

We solve the hydrodynamic equations for logarithmic density $\ln\rho$,
velocity $\uu$, and specific entropy $s$, in the form
\EQA
{\DD \ln \rho \over \DD t}&=&-\nab\cdot\uu, \nonumber \\
 \rho{\DD \uu\over \DD t}&=&-\nab P-2\OO\times\rho\uu
+\rho\grav_{\rm eff} + \nab\cdot(2\nu\rho\SSSS),\nonumber  \\
 \rho T {\DD S \over \DD
  t}&=&-\nab\cdot\FF_{\rm rad}+2\nu\rho\SSSS^2+{\cal H},
\nonumber
\label{sRT}
\ENA
where $P$ is the gas pressure, $\grav$ is the gravitational acceleration,
$\nu$ is the viscosity,
$\SSSS=\half[\nab\uu+(\nab\uu)^T]-\onethird\IIII\nab\cdot\uu$
is the traceless rate-of-strain tensor, $\IIII$ is the unit tensor,
$T$ is the temperature,
and $\FF_{\rm rad}$ is the radiative flux.
For the equation of state, we assume a perfect gas with
$p=({\cal R}/\mu)T\rho$, where ${\cal R}$ is the universal gas constant
and $\mu$ is the mean molecular weight.
The pressure is related to $\rho$ and $S$ via $P=\rho^\gamma\exp(S/\cv)$,
where the adiabatic index $\gamma=\cp/\cv$ is the ratio of specific
heats at constant pressure and constant volume, respectively, and
$\cp-\cv={\cal R}/\mu$.
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 $\SSSS$) is given
by the frequency-integrated Planck function, so $S=(\sigmaSB/\pi)T^4$,
where $\sigmaSB$ is the Stefan--Boltzmann constant.
The negative divergence of the radiative flux is then given by
\EQ
-\nab\cdot\FF_{\rm rad}=\kappa\rho \oint_{4\pi}(I-S)\,\dd\Omega,
\label{fff}
\EN
where $\kappa$ is the opacity per unit mass
(assumed independent of frequency) and $I(\xx,t,\nnn)$
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 $\nnn$ through a solid angle $\dd\Omega$.
We obtain $I(\xx,t,\nnn)$ by solving the radiative transfer equation,
\EQ
-\nnn\cdot\nab I=\kappa\rho\, (I-S),
\label{RT-eq}
\EN
along a set of rays in different directions $\nnn$ 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
\EQ
\kappa=\kappa_0\rho^a T^b,
\label{kappa}
\EN
where $a$ and $b$ are considered as free parameters that characterize
the relevant radiative processes.
We compute the radiative conductivity $K(\rho,T)$,
which is given by
\EQ
K(\rho,T)={16\sigmaSB T^3\over 3\kappa\rho}
={16\sigmaSB T^{3-b}\over 3\kappa_0\rho^{a+1}}.
\label{K-model}
\EN
We note that, in a plane-parallel polytropic atmosphere,
$T(z)$ varies linearly with height $z$ and in the stationary state,
$K(\rho,T)$ is constant in the optically thick part.
%AB: is that actually the case?
This implies that $\rho$ is proportional to $T^n$, where
\EQ
n={3-b\over 1+a}
\label{n_from_ab}
\EN
is the polytropic index
(not to be confused with the direction of a ray $\nnn$).

\section{Comparison with vertically integrated disks}

Use
\begin{equation}
{4\sigmaSB T_c^{4-b}\over3\kappa_0\Sigma^{1+a}}H
=\half\nut\Sigma(\threehalf\Omega)^{2}
\end{equation}
or
\begin{equation}
{4\sigmaSB ({\cal R}/\mu)^{-(4-b)}\cs^{2(4-b)}\over3\kappa_0\Sigma^{1+a}}H
=\half\nut\Sigma(\threehalf\Omega)^{2}
\end{equation}
which leads to
\begin{equation}
\cs^{2(4-b)}H={\cal K}\nut\Sigma^{2+a}\Omega^2
\end{equation}
where
\begin{equation}
{\cal K}={27\kappa_0\over32\sigmaSB}\left({{\cal R}\over\mu}\right)^{4-b}
\end{equation}
Eliminate sound speed and $\nut$, so
\begin{equation}
\Omega^{2(4-b)} H^{2(4-b)}H={\cal K}\alphaSS\Omega H^2\Sigma^{2+a}\Omega^2
\end{equation}
collect
\begin{equation}
\Omega^{5-2b} H^{7-2b} ={\cal K}\alphaSS \Sigma^{2+a}
\end{equation}
Combine with angular momentum conservation
$\nut\Sigma=\dot{m}$, so
\begin{equation}
\alphaSS\Omega H^2\Sigma=\dot{m}\equiv\dot{M}/3\pi
\label{mdoteqn}
\end{equation}
to eliminate $\Sigma$, so we find
\begin{equation}
\Omega^{5-2b} H^{7-2b} ={\cal K}\alphaSS (\dot{m}/\alphaSS\Omega H^2)^{2+a}
\end{equation}
or
\begin{equation}
\Omega^{2+a} H^{2(2+a)}
\Omega^{5-2b} H^{7-2b} ={\cal K} \dot{m}^{2+a}\alphaSS^{-(1+a)}
\end{equation}
so
\begin{equation}
\Omega^{7+a-2b} H^{11+2a-2b}={\cal K} \dot{m}^{2+a}\alphaSS^{-(1+a)}
\end{equation}
Using $\Omega(r)=r^{3/2}/\sqrt{GM}$, we have $H=H(r)$, given by
\begin{equation}
H^{11+2a-2b}={\cal K} \dot{m}^{2+a}\alphaSS^{-(1+a)} \Omega^{-(7+a-2b)} 
\end{equation}
Alternatively, using \Eq{mdoteqn},
we can determine $\Sigma=\Sigma(r)$ as
\begin{equation}
(\dot{m}/\alphaSS\Omega\Sigma)^{(11+2a-2b)/2}={\cal K} \dot{m}^{2+a}\alphaSS^{-(1+a)} \Omega^{-(7+a-2b)} 
\end{equation}
so
\begin{equation}
\dot{m}^{7/2-b}={\cal K} \alphaSS^{9/2-b} \Omega^{-3/2+b} 
\Sigma^{11/2+a-b}
\end{equation}
or
\begin{equation}
\Sigma^{11/2+a-b}=\dot{m}^{7/2-b}{\cal K}^{-1} \alphaSS^{-9/2+b} \Omega^{3/2-b}
\end{equation}
The important relation is $\nut$ versus $\Sigma$, so now we find
\begin{equation}
\dot{m}\propto \Sigma^{(11+2a-2b)/(7-2b)}
\end{equation}
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.,
\begin{equation}
b>11/2+a = 13/2 = 6.5 \quad \mbox{(for $a=1$)}
\end{equation}
which is just a bit less than the value for $\Hminus$ opacity, or
\begin{equation}
b> 6 \quad \mbox{(for $a=1/2$)},
\end{equation}
which is more appropriate for $\Hminus$ opacity.

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{all_slopes}
\end{center}\caption[]{
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).
}\label{all_slopes}\end{figure}

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{pSummary_analyt}
\end{center}\caption[]{
Analytic radiative equilibrium solutions
}\label{pSummary_analyt}\end{figure}

\begin{table}[htb]\caption{
}\vspace{12pt}\centerline{\begin{tabular}{cll}
variable & code units & cgs units \\
\hline
$\Sigma$  & $\g\cm^{-3}\Mm$  &  $10^8\g\cm^{-2}$ \\
$\dot{M}$ & $\g\cm^{-3}\Mm^2\km\s^{-1}$  &  $10^{21}\g\s^{-1}$ \\
          &                              & $\approx1.5^{-5}\MSun\yr^{-1}$ \\
\label{Ttimescale}\end{tabular}}\end{table}

\section{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.

%\begin{figure}[t!]\begin{center}
%\includegraphics[width=\columnwidth]{pSummary}
%\end{center}\caption[]{
%Radiative equilibrium solutions
%}\label{pSummary}\end{figure}

relative hydro stress:
\begin{equation}
\overline{\rho u_x u_y}/
\sqrt{\overline{\rho u_x^2}\;\overline{\rho u_y^2}}
\end{equation}
is around $-0.34$

Alpha:
\begin{equation}
\alpha_{\rm SS}^{\rm K}=\overline{\rho u_x u_y}/
\threehalf\meanrho\Omega^2H^2
\end{equation}
\begin{equation}
\alpha_{\rm SS}^{\rm M}=-\overline{B_x B_y/\mu_0}/
\threehalf\meanrho\Omega^2H^2
\end{equation}
or
\begin{equation}
\alpha_{\rm SS}^{\rm K}=\overline{u_x u_y}/
\threehalf\Omega^2H^2
\end{equation}
\begin{equation}
\alpha_{\rm SS}^{\rm M}=-\overline{B_x B_y/\rho\mu_0}/
\threehalf\meanrho\Omega^2H^2
\end{equation}
\section{Hydrostatic, non-turbulent disk structure}

\begin{equation}
{\dd F\over\dd z}={\cal H}(z)
\end{equation}
\begin{equation}
{\dd T\over\dd z}=-F/K
\end{equation}
\begin{equation}
{\dd P\over\dd z}=-\rho g
\end{equation}
where $g=\Omega^2z$ is the local gravitational acceleration.

Assume ${\cal H}=\const$ for $z\leq L_{\rm heat}$ and 0 otherwise,
so $F={\cal H}\min(z,L_{\rm heat})$.
Then
\begin{equation}
{\dd T\over\dd z}=-{\cal H}\min(z,L_{\rm heat})/K
\end{equation}
\begin{equation}
{\dd P\over\dd z}=-\rho \Omega^2 z
\end{equation}
Divide by each other
\begin{equation}
{\dd T\over\dd P}={{\cal H}\over K\rho \Omega^2}\,\min(L_{\rm heat}/z,1)
\end{equation}
\begin{equation}
{\dd\ln T\over\dd\ln P}={{\cal H} P\over K\rho T\Omega^2}\,\min(L_{\rm heat}/z,1)
\end{equation}
Use $P=\Rgas T\rho/\mu$, and define
\begin{equation}
\nabrad={{\cal H} \Rgas/\mu\over K\Omega^2}\,\min(L_{\rm heat}/z,1).
\end{equation}
So we would need to solve $\dd\ln T/\dd\ln P=\nabrad$,
but with convection we have mixing, so the  actual temperature
gradient becomes equal to the adiabatic one and therefore
\begin{equation}
{\dd\ln T\over\dd\ln P}=\min(\nabrad,\nabad)
\end{equation}
and $\nabad=1-1/\gamma$.

\begin{figure*}[t!]\begin{center}
\includegraphics[width=\textwidth]{pSummary_alp_radiative}
\end{center}\caption[]{
S curve for $\alpha_{\rm SS}=5\times10^{-3}$ (red),
$2\times10^{-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
$\Sigma=3\times10^{-7}\g\cm^{-3}\Mm=0.3\g\cm^{-2}$.
}\label{pSummary_alp_radiative}\end{figure*}

\begin{figure*}[t!]\begin{center}
\includegraphics[width=\textwidth]{pSummary_alp}
\end{center}\caption[]{
S curve for $\alpha_{\rm SS}=5\times10^{-3}$ (red),
$2\times10^{-3}$ (black), and
with convection and $\alpha_{\rm SS}=5\times10^{-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.
}\label{pSummary_alp}\end{figure*}

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{pstrat_comp}
\end{center}\caption[]{
%Comparison of solutions for $\Sigma=3.6\times10^{-6}$.
%AB: must have been a mistake, as is evident from the first panel
Comparison of solutions for $\Sigma=3.6\times10^{-7}$.
$\kappa_0=2\times10^{12}$, with Kramers/$\Hminus$ opacity
using $a_{\Hminus}=1$, $b_{\Hminus}=4$, and $T_0=2\times10^4\K$,
$\rho_0=10^{-5}\g\cm^{-3}$.
Note that the models on the upper and middle branches are
convectively unstable (top-heavy and negative vertical entropy gradient).
}\label{pstrat_comp}\end{figure}

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{pstratification_comp}
\end{center}\caption[]{
Comparison between analytic solution (dashed, black) and
{\sc Pencil Code} solution (solid, red) for the H state solution
with ${\cal H}_0=1.4\times10^{-5}$.
}\label{pstratification_comp}\end{figure}

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{pstratification_comp_strat_H6.0e-06.ps}
\end{center}\caption[]{
Same as \Fig{pstratification_comp}, but for the M solution
with ${\cal H}_0=6\times10^{-6}$.
Note that the parameters $\Sigma$ and $\alpha_{\rm SS}$ are the same
in both cases.
}\label{pstratification_comp_strat_H6.0e-06}\end{figure}

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{pnablas}
\end{center}\caption[]{
Comparison of different $\nabla$ as a function of $z$
for ${\cal H}_0=10^{-3}$, $\dot{M}=5.8\times10^{-5}$,
and $\Sigma=6.3\times10^{-6}$.
}\label{pnablas}\end{figure}

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{pradconvflux}
\end{center}\caption[]{
Ratio $\nabla/\nabrad$, which corresponds to the
ratio of radiative to total flux when $\nabla/\nabrad<0$.
(back) and ratio of convective to total flux (red)
in the improved mixing length model.
}\label{pradconvflux}\end{figure}

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{pcomp}
\end{center}\caption[]{
Growth of Alfv\'en speed in a turbulent MRI simulation.
(Different solutions are patched together.)
}\label{pcomp}\end{figure}

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{peps}
\end{center}\caption[]{
Energy dissipation (magnetic red, and kinetic blue)
as a function of $t$.
No external heating has been applied.
}\label{peps}\end{figure}

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{pepsKMz}
\end{center}\caption[]{
Energy dissipation (magnetic red, and kinetic blue)
as a function of $z$.
}\label{pepsKMz}\end{figure}

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{puz_slice_MCV288_h1em5d4_S10_nu5em2_eta5em2a_noheat2}
\end{center}\caption[]{
Velocity and entropy fluctuation in a vertical plane.
The $\tau=1$ surface is shown as a dashed line.
}\label{puz_slice_MCV288_h1em5d4_S10_nu5em2_eta5em2a_noheat2}\end{figure}

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{Bz60}
\end{center}\caption[]{
Visualization of $B_z$ on the periphery of the computational domain.
}\label{Bz60}\end{figure}

In \Fig{Bz60} 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 \citep{refs}
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.

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{pTTmz_MCV288_h1em5d4_S10_nu5em2_eta5em2a_noheat2}
\end{center}\caption[]{
Mean temperature versus $t$ and $z$.
}\label{pTTmz_MCV288_h1em5d4_S10_nu5em2_eta5em2a_noheat2}\end{figure}

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{pstressKM}
\end{center}\caption[]{
Total, magnetic (red), and kinetic (blue) stresses.
}\label{pstressKM}\end{figure}

\section{Disk parameters}

\EQ
H=\cs^{\max}/\Omega
\EN
\EQ
\dot{M}=3\pi\;{2{\cal H}_0L_{\cal H}\over(\threehalf\Omega)^2}
\EN
factor 2?
\EQ
\alpha_{\rm SS}=\dot{M}/3\pi\;\Sigma\Omega H^2
\EN
\EQ
\sigmaSB\Teff^4=(\Mdot/3\pi)\,	(\threehalf\Omega)^2
\EN

\section{Comparison with time integration}

In \Figs{pstratification_comp}{pstratification_comp_strat_H6.0e-06}
we verify solutions with Pencil 1D.

$\Delta P=\Delta z^2/2\HP^2$
and $\HP=\sqrt{{\cal R}T/\mu}/\Omega$.

\section{Problems with vertical disk integration}

\EQ
-{16\sigmaSB T^{3-b}\over3\kappa_0\rho^{1+a}}\,{\dd T\over\dd z}
=\int{\cal H}\,\dd z
\EN
\EQ
{4\sigmaSB T^{4-b}\over3\kappa_0\Sigma^{1+a}}H^{a}=\nut\Sigma (\threehalf\Omega)^2
\EN

\section{Mixing-length treatment}

A relatively simple treatment of convection is to set $\nabla$ equal to
$\nabad$ when $\nabrad>\nabad$ and equal to $\nabrad$ 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 the convective
energy flux is adopted.
In standard mixing length theory (MLT), the convective energy flux is
proportional to the superadiabatic gradient, $\nabla-\nabad$, which would
obviously be zero in the primitive approach of setting $\nabla=\nabad$
in the convection zone.
Instead, we now set the total flux equal to the sum of radiative and
convective energy fluxes.
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
\EQ
F_{\rm tot}=F_{\rm rad}+F_{\rm conv}.
\EN
The radiative flux is still given by \Eq{XX} and $\nabla$ is still
the actual logarithmic temperature gradient, but now it is no longer
equal to $\nabrad$.
Instead, ...

\begin{itemize}
\item
what is the KH time scale for given $\kappa_0$?
How does it depend on $\rho$?
\item
How is $T_{\rm mid}$ related to $\Sigma$?
\item
Discuss strategy of fixing ${\cal H}$.
\item
Plot $T$ vs $\tau$.
\end{itemize}

%r e f
\begin{thebibliography}{}

\bibitem[Abramowicz et al.(1996)]{ABL96}
Abramowicz, M. A., Brandenburg, A., \& Lasota, J.-P.\ymn{1996}{281}{L21}

\bibitem[Balbus \& Hawley(1991)]{BH91}
Balbus, S. A., \& Hawley, J. F.\yapj{1991}{376}{214}

\bibitem[Balbus \& Hawley(1992)]{BH92}
Balbus, S. A. \& Hawley, J. F.\yapj{1992}{400}{610}

\bibitem[Bath \& Pringle(1981)]{BP81}
Bath, G. T., \& Pringle, J. E.\ymn{1981}{194}{967}

\bibitem[Bath \& Pringle(1982)]{BP82}
Bath, G. T., \& Pringle, J. E.\ymn{1982}{199}{267}

\bibitem[Brandenburg et al.(1995)]{BNST95}
Brandenburg, A., Nordlund, \AA., Stein, R. F., \& Torkelsson, U.\yapj{1995}{446}{741}

\bibitem[Cannizzo et al.(1982)]{Cannizzo82}
Cannizzo, J. K., Ghosh, P., \& Wheeler, J. C.\yapjl{1982}{260}{L83}

\bibitem[Cannizzo(1992)]{Cannizzo92}
Cannizzo, J. K.\yapj{1992}{385}{94}

\bibitem[Hawley \& Balbus(1991)]{HB91}
Hawley, J. F. \& Balbus, S. A.\yapj{1991}{376}{223}

\bibitem[Hawley \& Balbus(1992)]{HB92}
Hawley, J. F., Balbus, S. A.\yapj{1992}{400}{595}

\bibitem[Hawley et al.(1999)]{HBW96}
Hawley, J. F., Balbus, S. A., \& Winters, W. F.\yapj{1999}{518}{394}

\bibitem[Hawley et al.(1995)]{HGB95}
Hawley, J. F., Gammie, C. F., \& Balbus, S. A.\yapj{1995}{440}{742}

\bibitem[Hawley et al.(1996)]{HGB96}
Hawley, J. F., Gammie, C. F., \& Balbus, S. A.\yapj{1996}{464}{690}

\bibitem[Hirose et al.(2014)]{Hirose14}
Hirose, S., Blaes, O., Krolik, J. H., Coleman, M. S. B., \& Sano, T.\yapj{2014}{787}{1}

\bibitem[Meyer \& Meyer-Hofmeister(1981)]{MMH81}
Meyer, F. \& Meyer-Hofmeister, E.\yana{1981}{104}{L10}

\bibitem[Meyer \& Meyer-Hofmeister(1982)]{MMH82}
Meyer, F. \& Meyer-Hofmeister, E.\yana{1982}{106}{34}

\bibitem[Pessah et al.(2007)]{PCP07}
Pessah, M. E., Chan, C.-k., \& Psaltis, D.\yapj{2007}{668}{L51}

\bibitem[Shakura \& Sunyaev(1973)]{SS73}
Shakura, N. I., \& Sunyaev, R. A.\yana{1973}{24}{337}

\bibitem[Shaviv \& Wehrse(1991)]{SW91}
Shaviv, G., \& Wehrse, R.\yana{1991}{251}{117}

\bibitem[Stone et al.(1996)]{SHGB96}
Stone, J. M., Hawley, J. F., Gammie, C. F., \& Balbus, S. A.\yapj{1996}{463}{656}

\end{thebibliography}

\vfill\bigskip\noindent\tiny\begin{verbatim}
$Header: /var/cvs/brandenb/tex/helena/SSdisk/notes.tex,v 1.37 2020/03/18 18:01:56 brandenb Exp $
\end{verbatim}

\end{document}
