% $Id: notes.tex,v 1.28 2018/04/17 14:03:36 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,url}
\graphicspath{{./fig/}{./png/}}
\def\blue{\textcolor{blue}}
\input gafd_journ
\input apj
\title{Radiative accretion disk models}
\author{Notes and calculations assembled by
A.\ Brandenburg and U.\ Das}
\date{\today,~ $ $Revision: 1.28 $ $}
\begin{document}
\maketitle

\section{Background}

\subsection{Cataclysmic variables}

[The notes in this subsection are excerpts from an excellent wikipedia
article\footnote{\url{https://en.wikipedia.org/wiki/Cataclysmic_variable_star}}.]
Cataclysmic variables (CV) consist of a white dwarf (WD) primary and a
mass transferring secondary (donor star).
Strong UV and X-ray emission is often seen from the accretion disc,
powered by the loss of gravitational potential energy from the infalling
material.
Material at the inner edge of disc falls onto the surface of the white
dwarf primary.
A classical nova outburst occurs when the density and temperature at
the bottom of the accumulated hydrogen layer rise high enough to ignite
runaway hydrogen fusion reactions.
The accretion disc may be prone to an instability leading to dwarf nova
outbursts, when the outer portion of the disc changes from a cool, dull
mode to a hotter, brighter mode for a time, before reverting to the cool
mode. Dwarf novae can recur on a timescale of days to decades.

\subsection{Relevant literature}

In a recent paper, \cite{Hirose14} presented shearing box simulations
of magneto-rotational instability (MRI) turbulence in a regime relevant
to dwarf novae and soft X-ray transient outbursts.
They find two stable thermal equilibria in the effective temperature
versus surface mass density diagram, which is consistent with the
hypothetical S-curve dependence adopted since the 1980s
\citep{BP81,BP82,MMH81,MMH82,Cannizzo82}.
In their new work, \cite{Hirose14} find that convection strengthens the
dynamo and enhances the Shakura-Sunyaev $\alpha$ parameter \citep{SS73}
by generating vertical magnetic field that seed the axisymmetric MRI
and by increasing cooling.
Our goal is to provide an independent confirmation of the results
of \cite{Hirose14} and to arrive at a better understanding of the
transition between the two branches that help parameterizing
the relevant physics, which can then be used in simpler one-dimensional
accretion disk models.

\section{Polytropic solution}
\subsection{Hydrostatic balance}

We expect accretion disks to be turbulent, e.g., because of MRI or
because of convection in the outer layers.
In both of these cases, the cause of turbulence is related to the
occurrence of an instability of the hydrostatic base state,
Once turbulence develops, the original base state would no longer
be directly relevant, because it would become strongly modified.
To understand this transition, we need to establish first the
hydrostatic reference state.

Hydrostatic balance implies\footnote{see the manual to the
{\sc Pencil Code} on \url{https://github.com/pencil-code}}
\begin{equation}
0=-\nab h+T\nab S-\nab\phi,
\end{equation}
where $T$ is temperature, $S$ is specific entropy, $h=\cp T$
is the enthalpy for a perfect with fixed ionization with
$\cp$ being the specific heat at constant pressure, and
$\phi=\half(z_\infty-z)\,\Omega^2$ is the gravitational potential
with $z$ being the height from the midplane, $z_\infty$ is the
height of the surface, and $\Omega$ is the angular velocity
at the position where the our local coordinate system is situated.

For a polytropic solution, the $T\nab S$ can be written underneath
a gradient and can be combined with the enthalpy to give what we
might call the pseudo-enthalpy, so
\begin{equation}
0=-\nab(\tilde{h}+\phi)
\end{equation}
so
\begin{equation}
\tilde{h}={1-{1\over\gamma}\over1-{1\over\Gamma}}\,c_pT
={1\over2}\left(z_\infty^2-z^2\right)\Omega^2,
\end{equation}
where $\gamma=\cp/\cv$ is the ratio of specific heats,
$\cv$ is the specific heat at constant volume, and
$\Gamma$ is a parameter characterizing the polytropic
solution and is related to the polytropic index;
see \App{AppA} for details.
For monatomic gas we have $\gamma=5/3$.
If we assume a Kramers-like opacity, then $n=3.25$,
so $\Gamma=1+1/n=1.308$; see \cite{BB14} for details.
Thus, we have
\begin{equation}
\tilde{h}=(n+1)\left(1-{1\over\gamma}\right)\,c_pT
={1\over2}\left(z_\infty^2-z^2\right)\Omega^2,
\end{equation}
and therefore
\begin{equation}
T={\textstyle{5\over17}}\left(z_\infty^2-z^2\right)\Omega^2/c_p.
\end{equation}

\subsection{Thermal balance}

To sustain thermal equilibrium, we have to include a heating term
${\cal H}$, so the time-dependent entropy equation takes the form
\begin{equation}
\rho T {\DD s\over\DD t}=\nab\cdot K\nab T+{\cal H}
\end{equation}
Here we choose
\begin{equation}
{\cal H}=5\times10^{-3}\g\cm^{-3}\Km^3\s^{-3}\Mm^{-1},
\end{equation}
Using $\sigmaSB=5.67\times10^{-20}$ \citep[in code units; see][]{BB14},
we have
\begin{eqnarray}
T_{\rm surf} \!&=&\!
\left({\textstyle{3\over4}} z_{\rm heat}{\cal H}/\sigmaSB\right)^{1/4}
\nonumber \\ \!&=&\!
\left(0.75\times5\times10^{-3}/5.67\times10^{-20}\right)^{1/4}
\nonumber \\ \!&=&\!
23980\K.
\end{eqnarray}
We find
\begin{equation}
K=7.14\times10^{-6}
\end{equation}
For details about radiation transfer, see \cite{BB14} and references therein.
At some point, ionization will need to be included; see \cite{BB16}.

\begin{figure}[h!]\begin{center}
\includegraphics[width=\columnwidth]{pcomp}
\end{center}\caption[]{
$\kappa_0=10^4$, $z_{\max}=12$,
$z_{\rm heat}=5$, ${\cal H}=5\times10^{-3}$
$a=1$, $b=-3.5$. 144 meshpoints.
}\label{pcomp}\end{figure}

\begin{figure}[h!]\begin{center}
\includegraphics[width=\columnwidth]{pvar_HSgzDO_144c3}
\end{center}\caption[]{
$\kappa_0=10^5$, $z_{\max}=8$,
$z_{\rm heat}=5$, ${\cal H}=10^{-4}$
$a_{\rm Kr}=1$, $b_{\rm Kr}=-3.5$,
$a_\Hminus=1$, $b_\Hminus=4$.
144 meshpoints.
}\label{pvar_HSgzDO_144c3}\end{figure}

\subsection{Sketchy notes}

Tall box in \Fig{pcomp_tall}.
Time step is less stringent when the opacity floor is
allowed to be smaller, but it is quadratic
in the invert mean free path, $\delta t\propto\ell^2\propto(\kappa\rho)^{-2}$.
Unclear.

\begin{figure}[h!]\begin{center}
\includegraphics[width=\columnwidth]{pcomp_tall}
\end{center}\caption[]{
Floor values 0.2 and 0,02, kappa0=1e5
H=1e-4, time steps 4e-7 and 4e-6, respectively.
288 meshpoints, 1D, but time still short because of short timestep.
Should plot s profile.
}\label{pcomp_tall}\end{figure}

2-D convection, $u_\downarrow=1.2\Km\s^{-1}$,
$u_\uparrow=0.8\Km\s^{-1}$.
Note the extreme aspect ratio.

\begin{figure}[h!]\begin{center}
\includegraphics[width=\columnwidth]{pcomp_strat_144f}
\end{center}\caption[]{
Stratification, 144f
}\label{pcomp_strat_144f}\end{figure}

\begin{figure}[h!]\begin{center}
\includegraphics[width=\columnwidth]{puz_slice_288a}
\end{center}\caption[]{
Vertical slice, 2D, 288x288, time step 2e-3.
Initial amplitude 0.1 km/s.
$t=4\ks$.
$\nu=0.02\Mm\Km\s^{-1}$.
}\label{puz_slice_288a}\end{figure}

%\begin{figure}[h!]\begin{center}
%\includegraphics[width=\columnwidth]{pcomp_strat_288a}
%\end{center}\caption[]{
%Stratification, 288a
%}\label{pcomp_strat_288a}\end{figure}

\begin{figure}[h!]\begin{center}
\includegraphics[width=\columnwidth]{puz_slice_288c_t200}
\end{center}\caption[]{
Vertical slice, 2D, 288x288, time step 8e-4.
$t=200\ks$
$\nu=0.01\Mm\Km\s^{-1}$.
$\kappa_0=3\times10^5$, ${\cal H}=10^{-4}$.
}\label{puz_slice_288c_t200}\end{figure}

\begin{figure}[h!]\begin{center}
\includegraphics[width=\columnwidth]{puz_slice_288d_H1em5_nu5em3_t280}
\end{center}\caption[]{
Vertical slice, 2D, 288x288, time step 1e-3.
$t=3600\ks$.
$\nu=0.005\Mm\Km\s^{-1}$.
$\kappa_0=10^6$, ${\cal H}=10^{-5}$.
}\label{puz_slice_288d_H1em5_nu5em3_t280}\end{figure}

\begin{figure}[h!]\begin{center}
%\includegraphics[width=\columnwidth]{puz_slice_288d_H1em5_t10}
%\includegraphics[width=\columnwidth]{puz_slice_288d_H1em5_t50}
\includegraphics[width=\columnwidth]{puz_slice_288d_H1em5_t70}
\includegraphics[width=\columnwidth]{puz_slice_288d_H1em5_t443}
\end{center}\caption[]{
Vertical slice, 2D, 288x288, time step 8e-4.
$t=700\ks$, and $4200\ks$.
$\nu=0.01\Mm\Km\s^{-1}$.
$\kappa_0=10^6$, ${\cal H}=10^{-5}$.
}\label{puz_slice_288d_H1em5}\end{figure}

%\begin{table}[htb]\caption{
%}\vspace{12pt}\centerline{\begin{tabular}{lccccccc}
%\hline
%\label{Ttimescale}\end{tabular}}\end{table}

\begin{verbatim}
nu     kappa0  omeg  Run
5e-2   1e5           144f
5e-2   3e5           144g  never ran
2e-2   1e5           288a
1e-2   1e5     2.47  288b
\end{verbatim}

\subsection{New runs}

$\Sigma=2\times10^{-12}$, ${\cal H}=7.1\times10^{-7}$,
two stable disk states; see \Fig{pTT_last}.

\begin{figure}[h!]\begin{center}
\includegraphics[width=\columnwidth]{pTT_last}
\end{center}\caption[]{
$\Sigma=2\times10^{-12}$, ${\cal H}=5\times10^{-5}$,
two stable disk states.
}\label{pTT_last}\end{figure}

\section{CV disks}

The following parameters have been assembled from the
paper by \cite{SW91}.

$M=1\,M_\odot$, $R=10^8\m$, so $\Omega=\sqrt{GM/R^3}=12\ks^{-1}$,
$\Sigma=10^2\g\cm^{-2}=10^{-6}\Mm\g\cm^{-3}$.

We start with an isothermal model with $\cs=36\Km\s^{-1}$,
$\rho_0=10^{-6}\G\cm^{-3}$.
\\

\begin{verbatim}
kappa0 Heat ztau1 z [Mm]  taumax chimin run
1e10   1e-5 3.65  0 ... 4.5 289   3.2 CV288n
4e10   1e-5 5.37  0 ... 6.0 270   28  CV288q
1.6e11 5e-6 6.72  \pm 7.5   140  401  CV288q5
\end{verbatim}

\begin{verbatim}
M_WD=0.8 Msun
R_WD=10 Mm
R=(6.67e-11*2e30/12e-3^2)^.3333333333333
0.9e8 = 90 Mm.
\end{verbatim}

\begin{figure}[h!]\begin{center}
\includegraphics[width=\columnwidth]{formula_ii}
\end{center}\caption[]{
$\kappa_0=2\times10^{12}\cm^3\g^{-1}\Mm^{-1}$,
$a_\Hminus=1$, $b_\Hminus=4$.
}\label{formula_ii}\end{figure}

\begin{figure}[h!]\begin{center}
\includegraphics[width=\columnwidth]{ppxyaver_CV1152q9}
\end{center}\caption[]{
Run with
$\kappa_0=2\times10^{12}\cm^3\g^{-1}\Mm^{-1}$,
$a_\Hminus=1$, $b_\Hminus=4$.
$H=2\times10^{-5}$, $z_H=0.5\Mm$.
}\label{ppxyaver_CV1152q9}\end{figure}

\begin{figure}[h!]\begin{center}
\includegraphics[width=\columnwidth]{ppxyaver_CV1152q9_noheat}
\end{center}\caption[]{
Run with
$\kappa_0=2\times10^{12}\cm^3\g^{-1}\Mm^{-1}$,
$a_\Hminus=1$, $b_\Hminus=4$.
$H=0$
}\label{ppxyaver_CV1152q9_noheat}\end{figure}

\begin{figure}[h!]\begin{center}
\includegraphics[width=\columnwidth]{ppxyaver_CV288q10}
\end{center}\caption[]{
Run with
$\kappa_0=2\times10^{12}\cm^3\g^{-1}\Mm^{-1}$,
$a_\Hminus=1$, $b_\Hminus=4$.
${\cal H}=3\times10^{-6}$, $z_{\cal H}=1\Mm$.
}\label{ppxyaver_CV288q10}\end{figure}

\begin{figure}[h!]\begin{center}
\includegraphics[width=\columnwidth]{pstrat_CV288q10}
\end{center}\caption[]{
Stratification
}\label{pstrat_CV288q10}\end{figure}

%\begin{itemize}
%\item
%\end{itemize}

\appendix
\section{Relation between $\Gamma$ and $n$}
\label{AppA}

A polytropic stratification is one where the
stratifications of density and pressure, $\rho(z)$
and $P(z)$ are related to that of temperature $T(z)$
in a powerlaw fashion, i.e.,
\begin{equation}
\rho(z) \propto T(z)^n,\quad P(z) \propto T(z)^{n+1}.
\end{equation}
where the latter relation is a consequence of the perfect
gas equation, $P\propto\rho T$.
Thus, $n=\dd\ln\rho/\dd\ln T$ and $n+1=\dd\ln P/\dd\ln T$.
It is also useful to define
\begin{equation}
\Gamma=\dd\ln P/\dd\ln\rho={n+1\over n}=1+{1\over n},
\end{equation}
which is convenient as it can easily be compared with $\gamma$
(which is a property of the gas, not of the stratification,
as is $Gamma$).
To see that, let us calculate the vertical gradient of the
specific entropy $s$ using the formula
\begin{equation}
\dd s/\cp={1\over\gamma}\dd\ln P-\dd\ln\rho.
\end{equation}
Thus, we have
\begin{equation}
{\dd s/\cp\over\dd\ln P}={1\over\gamma}-{1\over\Gamma},
\end{equation}
which is also known as the superadiabatic gradient, $\nabla-\nabad$. where
$\nabla=\dd\ln T/\dd\ln P$ is the double-logarithmic temperature gradient.
Since $\dd\ln P/\dd z=-\HP$ is a negative quantity ($\HP$ is the pressure
scale height), a positive superadiabatic gradient corresponds to a negative
specific entropy gradient and thus to a Schwarzschild-unstable stratification.
Thus, $\nabla-\nabad>0$ means $\Gamma>\gamma$.

We have thus three closely related quantities,
$n=\dd\ln\rho/\dd\ln T$, $\Gamma=\dd\ln P/\dd\ln\rho$, and
$\nabla=\dd\ln T/\dd\ln P$, which are related to each other via
\begin{equation}
\Gamma=1+{1\over n}={1\over1-\nabla},
\end{equation}
\begin{equation}
n={1\over\Gamma-1}={1\over\nabla}-1,
\end{equation}
\begin{equation}
\nabla=1-{1\over\Gamma}={1\over n+1}.
\end{equation}
For $\gamma=5/3$, marginal stability corresponds to
$n=3/2=1.5$, $\Gamma=\gamma=5/3$, and $\nabla=5/2=2.5$.

%r e f
\begin{thebibliography}{}

\bibitem[Barekat \& Brandenburg(2014)]{BB14}
Barekat, A., \& Brandenburg, A.\yanaN{2014}{571}{A68}
{Near-polytropic stellar simulations with a radiative surface}

\bibitem[Bath \& Pringle(1981)]{BP81}
Bath, G. T., \& Pringle, J. E.\ymn{1981}{194}{967}
{986}{The evolution of viscous discs. I - Mass transfer variations}

\bibitem[Bath \& Pringle(1982)]{BP82}
Bath, G. T., \& Pringle, J. E.\ymn{1982}{199}{267}
{280}{The evolution of viscous discs. II - Viscous variations}

\bibitem[Bhat \& Brandenburg(2016)]{BB16}
Bhat, P., \& Brandenburg, A.\yanaN{2016}{587}{A90}
{Hydraulic effects in a radiative atmosphere with ionization}

\bibitem[Cannizzo et al.(1982)]{Cannizzo82}
Cannizzo, J. K., Ghosh, P., \& Wheeler, J. C.\yapjl{1982}{260}{L83}
{L86}{Convective accretion disks and the onset of dwarf nova outbursts}

\bibitem[Hirose et al.(2014)]{Hirose14}
Hirose, S., Blaes, O., Krolik, J. H., Coleman, M. S. B., \& Sano, T.\yapjN{2014}{787}{1}
{Convection Causes Enhanced Magnetic Turbulence in Accretion Disks in Outburst}

\bibitem[Meyer \& Meyer-Hofmeister(1981)]{MMH81}
Meyer, F. \& Meyer-Hofmeister, E.\yana{1981}{104}{L10}
{L12}{On the elusive cause of cataclysmic variable outbursts}

\bibitem[Meyer \& Meyer-Hofmeister(1982)]{MMH82}
Meyer, F. \& Meyer-Hofmeister, E.\yana{1982}{106}{34}
{42}{Vertical structure of accretion disks}

\bibitem[Shakura \& Sunyaev(1973)]{SS73}
Shakura, N. I., \& Sunyaev, R. A.\yana{1973}{24}{337}
{355}{Black holes in binary systems. Observational appearance}

\bibitem[Shaviv \& Wehrse(1991)]{SW91}
Shaviv, G., \& Wehrse, R.\yana{1991}{251}{117}
{132}{Continuous energy distributions of accretion discs}

\end{thebibliography}

\vfill\bigskip\noindent\tiny\begin{verbatim}
$Header: /var/cvs/brandenb/tex/upasana/raddisk/notes.tex,v 1.28 2018/04/17 14:03:36 brandenb Exp $
\end{verbatim}

\end{document}
