% This is the aspauthor.tex LaTeX file
% Copyright 2014, Astronomical Society of the Pacific Conference Series
% Revision:  14 August 2014

% To compile, at the command line positioned at this folder, type:
% latex aspauthor
% latex aspauthor
% dvipdfm aspauthor
% This will create a file called aspauthor.pdf.

\documentclass[11pt,twoside]{article}
\usepackage{./asp2014}

\aspSuppressVolSlug
\resetcounters

\bibliographystyle{asp2014}

\markboth{T.\ A.\ A.\ Sigut}{The \texttt{Bedisk} and \texttt{Beray} Circumstellar Disk Codes}

\begin{document}

\title{The \texttt{Bedisk} and \texttt{Beray} Circumstellar Disk Codes}

\author{T.\ A.\ A.\ Sigut
\affil{The University of Western Ontario, London, Ontario, Canada; 
\email{asigut@uwo.ca}}}

\paperauthor{T.\ A.\ A.\ Sigut}{asigut@uwo.ca}{0000-0002-0803-8615}
{The University of Western Ontario}{Physics and Astronomy}
{London}{Ontario}{N6C 3K9}{Canada}

\begin{abstract}

The \texttt{Bedisk} and \texttt{Beray} radiative transfer codes can be
used to compute synthetic spectra and images for stars surrounded by
circumstellar disks. The operation and atomic data requirements of
these codes are discussed, as well as applications to the classical 
Be stars.

\end{abstract}

\section{Classical Be Stars and Related Objects}

Circumstellar disks are ubiquitous; newly forming stars are surrounded
by accretion disks that act to transport material to the star and angular
momentum away. Such disks dissipate by the time these stars arrive
on the core hydrogen-burning main sequence; however,  some
massive, main sequence stars, corresponding to spectral type
B with masses between $3$ and $25\,\rm M_{\sun}$, will find themselves 
again surrounded by
a circumstellar disk in the course of their main sequence evolution.
About one fifth of all main sequence B stars show  observational
evidence for such disks \citep{ZorecBriot1997}, such as emission lines in their
spectra, an excess of infrared radiation compared to isolated stars,
and polarization of their light. Such stars are classified as ``B
emission-line" stars or Be stars. Extensive interferometric observations
of Be stars in the solar neighbourhood confirm that the circumstellar
material is in the form of a thin, equatorial disk in Keplerian rotation
about the central star \citep{Quirrenback1997}.

The Be stars are the most rapidly-rotating stars on the main sequence,
often with rotation periods of just a day or two. This rapid
rotation is thought to be the key driver of the process that ejects gas from the stellar
photosphere into a disk \citep{RiviniusReview2013}. Such
out-flowing disks are usually referred to as {\it decretion\/} disks.
Models of stellar evolution with rotation demonstrate that internal
redistribution of angular momentum via rotational mixing can produce
episodes of {\it critical rotation\/} in which the effective gravitational
acceleration at the stellar equator vanishes \citep{Granada2013}. Although
the detailed physics of the gas ejection remains unknown,
the outward transport of angular momentum via a decretion disk is
an efficient mechanism to help the central star rid itself of its excess angular
momentum \citep{Krticka2011}. After ejection, the subsequent evolution
of the disk is well described by hydrodynamical evolution governed by viscous torques
\citep{RiviniusReview2013}.

Be stars range from those visible to the unaided eye (such as
$\gamma\,$Cassiopeia) to those detected in nearby galaxies. Large
fractions of Be stars occur in the Large and Small Magellanic
Clouds \citep{Martayan2007}, likely due to the low
metallicities in these galaxies. For this reason, the Be stars are excellent laboratories
for studying both stellar evolution with rotation and the complex hydrodynamical
processes occurring in circumstellar disks. In addition to these Be stars,
often called {\it Classical\/} Be stars,  there are several other
populations of closely related stars: Herbig Ae/Be stars exhibit many of
the properties ascribed to the Be stars, with the added presence of
dust emission and association with newly forming B stars still surrounded by their
accretion disks. At the other end of the evolutionary scale are the
sgB[e] stars, which may be examples of disk ejection via critical rotation in post
main sequence evolution \citep{Granada2016}.

To study all of these systems, it is important to have computational
tools that can predict spectra, overall spectral energy distributions (SEDs),
and interferometric phases and visibilities, preferably in a unified way,
in order to test disk formation and evolution mechanisms. In this
contribution, I will describe one such code suite, \texttt{Bedisk/Beray}, 
and give examples of its application.


\section{\texttt{Bedisk}}

The \texttt{Bedisk} code \citep{SigutJones2007} requires as input a
user-specified density structure for the disk, $\rho(R,Z)$, and
the fundamental parameters of the central B star, including its
photoionizing radiation field. Then \texttt{Bedisk} enforces radiative
equilibrium to determine the gas temperature at each point in the
disk.  The main heating source of the disk
gas is photoionization by UV photons from the central star. At
each disk location, the radiative transfer equation is solved along rays
directed back to the central star, with the star's photospheric spectrum
taken as the ``upwind" boundary condition on the transfer equation. The
mean intensity at each disk location can then be used to compute the
photoionization heating. Cooling processes include escaping radiation
formed by recombination and collisionally-excited line radiation. Currently
\texttt{Bedisk} includes atomic models for the nine most abundant
elements, each over several ionization stages. As a by-product of
the radiative equilibrium calculation, atomic level populations are
obtained which can be used to construct the monochromatic opacity and
emissivity for use in further radiative transfer solutions.

A commonly-used disk density structure assumes an axisymmetric disk in
which the total gas density in the equatorial plane falls as a power-law
with distance. The vertical structure of the disk can be determined by
enforcing hydrostatic equilibrium with vertical component of the star's
gravitational acceleration. If we assume for the purposes of the density
model only that the temperature in the disk can be represented by some 
average temperature, $T_{\rm HE}$, and
that the disk is geometric thin, we obtain
\begin{equation}
\rho(R,Z) = \rho_0 \left(\frac{R_*}{R}\right)^{n} \, e^{-\left(\frac{Z}{H}\right)^2} \;.
\label{eq:rho}
\end{equation}
Here $R$ is the radial distance in the equatorial plane, $R_*$ is
the stellar radius, and $Z$ is the vertical height above or below the
plane. The base density of the disk, $\rho_0$, and the radial
power-law index, $n$, are free parameters. The disk
scale height is given by $H=H_0(T_{\rm
HE})\,[R/R_*]^{3/2}$. The ``hydrostatic" disk temperature, $T_{\rm HE}$, is
commonly taken to be 60\% of the central star's effective temperature.
However, \texttt{Bedisk} can also
integrate the hydrostatic equilibrium equation at each disk location $R$ in
a manner consistent
with the vertical temperature distribution, eliminating the assumption of a constant $T_{\rm HE}$ 
\citep{Sigut2009}. The simple scale
height model with a constant $T_{\rm HE}$ produces a flared-disk with $H\propto R^{3/2}$.
The more realistic treatment considerable alters this prediction, particularly near
the star.

%
% Figure set one.
%
%\articlefiguretwo{Sigut_Fig1.eps}{Sigut_Fig2_v3.eps}{fig1}
\articlefiguretwo{Sigut_Fig1a.eps}{Sigut_Fig1b.eps}{fig1}
{{\it Left}: The temperature distribution $T(R,Z)$ in a disk surrounding a B2V
star computed with \texttt{Bedisk}. The colourbar indicates temperature
in Kelvin. 
{\it Right}: Selected hydrogen, helium and iron emission lines
produced by this disk as seen by a distant observer at two different
inclination angles, $i=65^{\circ}$ (solid lines) and $i=85^{\circ}$
(dotted lines), computed by \texttt{Beray}. Note that for H$\beta$, the
wings of the underlying photospheric absorption line are still visible.}

Figure~\ref{fig1} shows $T(R,Z)$
for a disk with $\rho_0=7.5\,\cdot\,10^{-11}\rm g\,cm^{-3}$ and $n=2.0$
surrounding a (assumed spherical) B2V star with $M=9\,M_{\sun}$ and
$R=5\,R_{\sun}$. Note the very inhomogeneous temperature structure and the
appearance of a inner, cool region near the star.
At this location at UV frequencies, the optical depth
along all rays back to the star are optically thick and
the heating rate is reduced. Moving above (or below) this region
reduces the optical depths for some rays and the temperature
rises. Further away, the dense gas in the equatorial plane
cannot shield the entire star and this cool zone disappears. Overall,
the density-weighted average temperature in the disk is approximately
60\% of the central star's $T_{\rm eff}$; however, this does depend on
the disk density parameters and the assumed chemical composition of the
gas. Disks for Be stars in the Small Magellanic Cloud (SMC), for example, are
hotter than this scaling relation due to the lower metallicities in the
SMC \citep{Ahmed2012}.

An unrealistic aspect of Figure~\ref{fig1} is that the central star is
spherical. The rapid rotation of Be stars distorts their shape and leads
to the phenomena of gravitational darkening in which there is a temperature
gradient with latitude across the stellar surface, with a hotter
pole and cooler equator \citep{Espinosa2011}. This
effect is implemented in both \texttt{Bedisk} and \texttt{Beray}
\citep{McGill2011,McGill2013} and must be included in any realistic
analysis of Be star spectra \citep{AhmedSigut2017}.

\subsection{Atomic Data Needs of \texttt{Bedisk}}

In order to determine the microscopic rates of heating and cooling, the
rate of each atomic process must be explicitly computed, and the required atomic level
populations are found by enforcing statistical equilibrium.
Included in these rate equations are both collisional
processes that occur at the thermodynamic equilibrium rate set by
the local temperature and density, and radiative processes that couple
to the radiation field. Radiative processes typically do not
occur at the thermodynamic rate as the radiation field is non-local due
to photon loses from the boundaries and photon gains from other parts
of the disk. This overall treatment is usually referred to as a ``non-LTE"
calculation, although thermodynamic equilibrium is recovered when
non-local photons are unimportant. To compute the statistical equilibrium
equations, all radiative bound-bound transitions require Einstein
$A_{ji}$ values and all bound-bound collision transitions require thermally-averaged
collision strengths. Bound-free process require photoionization cross
sections and collisional ionization rates. Heating and cooling must also
include free-free processes (i.e.\ bremsstrahlung) but these do not enter into the
rate equations.  Of particular importance is photoionization as the ejected electrons are the principle source
of heating. The rate of photoionization of level $i$
of an atom or ion is given by
\begin{equation}
R_{i\kappa} = 4\pi\,\int_{\nu_0}^{\infty}\;J_{\nu}\,\sigma_{i\kappa}(\nu)\,\frac{d\nu}{h\nu} \;;
\end{equation}
therefore, accurate calculation requires a good estimate of both the mean
intensity, $J_{\nu}$, and the photoionization cross section. While the
central star in Be systems are typically hot, $T_{\rm eff}>10^4\;$K, the
disk gas can be considerable cooler, and photoionization of metals can make an important
contribution to the gas heating. This implies that autoionizing
resonances in $\sigma(\nu)$ can be important, although accurate resonant
positions are not as important as overall strengths.

\section{\texttt{Beray}}

Of course, one does not directly observe gas density or temperature. The local thermodynamic 
conditions $(T,\rho)$ must be used along with the atomic level populations to construct the 
gas opacity ($\chi_\nu$) and emissivity ($\eta_\nu$) at each point in the disk. 
These can then be used in a formal solution of 
the radiative transfer equation along a series of rays directed at the observer to predict the specific 
intensity. This
calculation is performed by the \texttt{Beray} code \citep{Sigut2011}.
The specific intensity from a patch of sky, $I_{\nu}(x,y)$, can be computed via
\begin{equation}
I_{\nu}=I^{\circ}_{\nu}\,e^{-T_{\nu}} + \int S_{\nu}\, e^{-\tau_{\nu}} d\tau_{\nu} \;
\end{equation}
where the integral over optical depth $\tau_{\nu}$ is along a straight ray passing through the system 
directed at the observer. $I^{\circ}_{\nu}$ is
the ``upwind" boundary condition on the intensity, $T_{\nu}$ is the total optical depth along the ray, and 
$S_{\nu}$ is the monochromatic source function, defined as the ratio of the gas emissivity to opacity at frequency
$\nu$. For rays that 
originate on the stellar surface, $I^{\circ}_{\nu}$ is fixed by the photospheric spectrum, otherwise
$I^{\circ}_{\nu}=0$.
The fundamental quantity produced by \texttt{Beray} is the monochromatic intensity image of the star-disk 
system, $I_{\nu}(x,y)$, where $x$ and $y$ are Cartesian coordinates on the sky. Some Be stars
are close enough to be resolved by interferometers, and in this case, the predicted interferometric visibilities and
phases follow from the Fourier transform of $I_{\nu}(x,y)$.
For more distant and unresolved
objects, one can only measure the integrated intensity, or flux, as a function of frequency,
\begin{equation}
F_{\nu} = \int I_{\nu}(x,y) \,\frac{dx\,dy}{D^2} \;,
\end{equation}
where the integral is over the image on the sky and $D$ is the distance to the star. This equation can be used to compute
both spectral line profiles for individual elements and overall SEDs.

The right panel of Figure~\ref{fig1} shows some example line profiles computed by \texttt{Beray} using the
disk temperature/density model of the left panel. Using the level populations produced by \texttt{Bedisk}, 
\texttt{Beray} can compute
profiles for transitions between any of the included energy levels. While H$\alpha$ is the most commonly-modelled
emission line in Be stars, many other elements also produce emission lines, including He\,{\sc i}, Fe\,{\sc ii},
and Mg\,{\sc ii}. In the profiles shown in Figure~\ref{fig1}, the underlying photospheric absorption line
is consistently included. For rays that terminate on the stellar surface, an appropriately Doppler shifted (due to
the star's rotation), non-LTE absorption line profile is used as the upwind boundary condition on the 
transfer equation. Doppler shifts along
the ray due to the disk's rotation are included in the transfer equation in the observer's frame.
Figure~\ref{fig1} also shows that the inclination of the system, the angle between the line of sight and the
stellar rotation axis, has a strong controlling effect on the observed spectrum. Doubly-peaked emission
lines are predicted for intermediate inclinations ($i=65^{\circ}$ is shown), which become deep shell absorption
for higher inclinations ($i=85^{\circ}$ is shown). In a shell line, the line centre flux is reduced as
the rays pass through the densest part of the disk. In Figure~\ref{fig1}, this is particularly evident for H$\beta$.

\section{Application to the Classical Be Stars}


The \texttt{Bedisk/Beray} code suite has been used extensively to
model the emission line spectra, spectral energy distributions (SEDs), and
interferometric visibilities of Be stars. The focus of much of this
work is on the modelling of H$\alpha$ emission line profiles to determine
the range of disk density parameters $(\rho_0,n)$ appropriate to Be
star disks \citep{Silaj2010,Silaj2014}. The power law index $n$ seems to
be associated with the dynamical state of the disk \citep{Vieira2017}, either dissipating
($1.5\!\!<\!\!n\!\!<\!\!3$), steady, ($3\!<\!n\!<\!3.5$), or building ($n\!>\!3.5$). In addition,
$\rho_0$ and $n$ can be used to estimate the disk mass if an estimate
of the disk radius is available. \cite{Arcos2017} use the disk radius
containing 90\% of the H$\alpha$ light to determine the H$\alpha$ disk
masses for a large sample of southern-sky Be stars. \cite{Sigut2015}
and \cite{JonesSigut2017} analyzed interferometric observations for
the bright Be stars o~Aqr (B7IVe) and 48~Per (B3Ve),
which spatially resolve the disks, to determine disk masses
(Figure~\ref{fig2}, left panel). All of
these results are compared by \cite{Arcos2017} to the predictions of
critically-rotating stellar evolutionary models \citep{Granada2013}
computed under the assumption that disk eject acts to remove sufficient
angular momentum from the star to prevent super-critical rotation. As shown
in the right panel of Figure~\ref{fig2}, the
overall trend of disk mass with stellar mass is similar, but the disk
masses (which are lower limits) generally lie above the predicted masses.

%
% Figure group two...
%
%\articlefiguretwo{StarVis48Per3.eps}{DiskMass5.eps}{fig2}{
\articlefiguretwo{Sigut_Fig2a.eps}{Sigut_Fig2b.eps}{fig2}{
{\it Left}: Fit to the observed interferometric visibilities of the Be star
48~Per. The \texttt{Beray} model points are for $\rho_0=10^{-10}\;\rm g\,cm^{-3}$,
$n=3.0$, $i=45^\circ$ and $R_{\rm disk} = 50\,R_*$. The reduced $\chi^2$ of the fit
is 1.45. Adapted from \cite{JonesSigut2017}.
{\it Right}: Comparison of
the median H$\alpha$ disk masses in the BeSOS sample of southern Be stars to the stellar
evolution predictions of \cite{Granada2013}. The error bars represent the $1\sigma$ dispersion
in the observed disk masses. The individual Be stars 48~Per and o~Aqr are as indicated. Adapted from
\cite{Arcos2017}. 
}

As \texttt{Bedisk/Beray} can model lines for elements other
than hydrogen, precision abundance analysis of the
central B stars of Be systems can be attempted. Here ``corruption" of the photospheric
spectrum by circumstellar emission and gravitational darkening
is consistently modelled. \cite{AhmedSigut2017} searched for
nitrogen overabundances in a sample of Be stars from the MiMeS survey
\citep{Wade2014} as a diagnostic of rotational mixing. Nitrogen
overabundances were found in about $1/3$ of the sample, similar to
results found for normal B stars. This is a slightly puzzling result
as previous studies of normal B stars have focused on slow rotators. However,
it is consistent with \cite{Dunstall2010}, who found no difference in nitrogen abundance
between the B and Be stars of the SMC.

\newpage

\acknowledgements I would like to thank the University of Western Ontario's Faculty of Science and
Research Western for financial support.

\begin{thebibliography}{}

\bibitem[Ahmed \& Sigut (2017)]{AhmedSigut2017} Ahmed, A., \& Sigut, T.\ A.\ A., 2017 MNRAS, 471, 3398

\bibitem[Ahmed \& Sigut (2012)]{Ahmed2012} Ahmed, A., \& Sigut, T.\ A.\ A., 2012 ApJ, 744, 191

\bibitem[Arcos et al.\ (2017)]{Arcos2017} Arcos, C., Jones, C.\ E., Sigut, T.\ A.\ A., Kannan, S., \&
Cur\'{e}, M.\ 2017, ApJ 842, 48

\bibitem[Dunstall et al.\ (2011)]{Dunstall2010} Dunstall, P.\ R., Brott, I., Dufton, P.\ L., et al., 2011, A\&A. 536, A65

\bibitem[Espinosa Lara \& Rieutord (2011)]{Espinosa2011} Espinosa Lara, F., \& Rieutord, M., 2011, A\&A, 533, A43

\bibitem[Granada et al.\ (2016)]{Granada2016} Granada, A., Sigut, T.\ A.\ A., Jones, C.\ E., et al., 2017, ASP Conference Series Volume 508, 81

\bibitem[Granada et al.\ (2013)]{Granada2013} Granada, A., Ekstr\"{o}m, S., Georgy, C.\ et al., 2013, A\&A, 553, A25

\bibitem[Jones et al.\ (2017)]{JonesSigut2017} Jones, C.\ E., Sigut, T.\ A.\ A., Grzenia, B.\ J., Tycner, C., \& Zavala, R., 2017, ApJ, 843, 24

\bibitem[Krti\v{c}ka et al.\ (2011)]{Krticka2011} Krti\v{c}ka, J., Owocki, S.\ P., \& Meynet, G., 2011, A\&A, 527, A84

%\bibitem[Labeyrie et al.\ (2006)]{Labeyrie2006} Labeyrie, A., Lipson, S.\ G., \& Nisenson P., 2006, An Introduction to Optical Stellar Interferometry (Cambridge University Press: New York)

\bibitem[Martayan et al.\ (2007)]{Martayan2007} Martayan, C., Floquet, M., Hubert, A.\ M., et al., 2007, A\&A, 472, 577

\bibitem[McGill et al.\ (2013)]{McGill2013} McGill, M.\ A., Sigut, T.\ A.\ A., \& Jones, C.\ E., 2013, ApJS, 204, 2

\bibitem[McGill et al.\ (2011)]{McGill2011} McGill, M.\ A., Sigut, T.\ A.\ A., \& Jones, C.\ E., 2011, ApJ, 743, 111

\bibitem[Quirrenback et al.\ (1997)]{Quirrenback1997} Quirrenback, A., Bjorkman, K.\ S., Bjorkman, J.\ E., et al., 1997, ApJ, 479, 477

\bibitem[Rivinius et al.\ (2013)]{RiviniusReview2013} Rivinius, T., Carciofi, A.\ C., \& Martayan, C., 2013, A\&ARv 21, 69

\bibitem[Sigut et al.\ (2015)]{Sigut2015} Sigut, T.\ A.\ A., Tycner, C., Jansen, B., \& Zavala, R.\ T.\ 2015, ApJ, 814, 159

\bibitem[Sigut (2011)]{Sigut2011} Sigut, T.\ A.\ A., 2011, ASP Conference Series 272, 426

\bibitem[Sigut et al.\ (2009)]{Sigut2009} Sigut, T.\ A.\ A., McGill, M.\ A., \& Jones, C.\ E., 2009, ApJ, 699, 1973

\bibitem[Sigut \& Jones (2007)]{SigutJones2007} Sigut, T.\ A.\ A., \& Jones, C.\ E., 2007, ApJ, 668, 481

\bibitem[Silaj et al.\ (2014)]{Silaj2014} Silaj, J., Jones, C. E., Sigut, T. A. A., \& Tycner, C., 2014, ApJ 795, 82

\bibitem[Silaj et al.\ (2010)]{Silaj2010} Silaj, J., Jones, C. E., Tycner, C., Sigut, T.\ A.\ A., \& Smith, A.\ D., 2010, ApJS, 187, 228

\bibitem[Vieira et al.\ (2017)]{Vieira2017} Vieira, R. G., Carciofi, A. C., Bjorkman, J. E., et al., 2017, MNRAS, 464, 3071

\bibitem[Wade et al.\ (2014)]{Wade2014} Wade, G.\ A., Grunhut, J., Alecian, E.,  et al., 2014, IAU Symposium 302, 265

\bibitem[Zorec \& Briot (1997)]{ZorecBriot1997} Zorec, J., \& Briot, D., 1997, A\&A, 1997, 318, 443

\end{thebibliography}

\end{document}
