\chapter{The Physics Behind XSTAR}
\label{sec:physics}


\section{Assumptions}

In this section we describe the computational procedure,
assumptions, free parameters, and the quantities which are calculated.  
Chief among the assumptions is that each model consists 
of a spherical gas cloud with a point  
source of continuum radiation at the center.  Therefore it
implicitly assumes spherical symmetry and radially 
beamed incident radiation.  In principle, more complicated geometries
can be mimiced by adding the local emission from  various 
spherical sections with appropriately chosen conditions. 
Also important is the assumption that all physical processes 
affecting the state of the gas are in a steady-state, i.e. that 
the the timescales for variation in the gas density and 
illuminating radiation are long compared with timescales affecting 
all atomic processes and propogation of radiation within the gas.
The validity of this assumption in any given situation depends on 
the conditions there, such as the gas density, temperature, and 
degree of ionization, and can be evaluated by using a model 
assuming steady-state and then calculating atomic rates which can
(hopefully) justify the steady-state assumption a posteriori.

The primary difference between these models and atmospheric models lies
in the treatment of the radiation field.
In an optically thick atmosphere the state of the gas
at any point in the cloud is coupled to the state of the gas in a
large part of the   rest of the cloud by the continuum 
radiation field and, in the limit of very large optical depth, 
can affect the excitation and ionization by suppressing radiative 
free-bound (recombination) transitions.  We attempt to mimic some of 
these effects by assigning to each recombination event an escape 
probability, using an expression given in the following section.
We also calculate the transfer of radiation by assuming that diffuse
radiation emitted at each radius is directed radially outward or 
inward.  These assumptions will be described in more detail later in this
section.

A further assumption governs the treatment of the transport of 
radiation in spectral lines.  Over a  wide range of
plausible situations large optical depths occur in the cores   of
lines of abundant ions, which may be important in cooling
the gas.  In treating the transfer of these photons we make 
the (conventional) assumption of complete redistribution 
in the scattering, which assumes that the transfer of the 
line photons occurs in a spatial region very close to the 
point where the photons are emitted.  Therefore the line 
emission rates are multiplied by an escape probability 
using an expression given in the following section.
This factor is intended to simulate the line scattering in the immediate vicinity 
of the emission region, and it assumes that escape from this region 
occurs when the photon scatters into a frequency where the optical depth is 
less than unity.  Following escape from the local region, the line photon
is assumed to be subject to absorption by continuum processes which 
are treated using the same 2-stream transfer equation as for the 
continuum.

\section{Input}

The input parameters  are the source spectrum, the gas composition, and the gas density 
or pressure.  The spectrum of the central source of radiation is described by  the spectral 
luminosity, $L_{0\varepsilon}=Lf_\varepsilon$, where $L$ is the total
luminosity (in erg s$^{-1}$).  The spectral function, $f_\varepsilon$,
is normalized such that $\int_0^\infty f_\varepsilon d\varepsilon=1$  and 
may be of one   of a variety of types, including:  Thermal bremsstrahlung, 
$f_\varepsilon\sim {\rm exp}(-\varepsilon/kT)$;  blackbody,  
$f_\varepsilon\sim\varepsilon^3/[{\rm exp}(\varepsilon/kT) -1]$; 
or   power law, $f_\varepsilon\sim\varepsilon^{\alpha}$; or the user
may define the form of the ionizing continuum by providing a table of 
energies and fluxes.   The gas consists of the elements H, He, C, N, O, Ne, Mg, 
Si, S, Ar, Ca, Fe, and Ni  with  relative abundances specified by the user.  
The default abundances are the solar values given by \cite{Grevesse1996}.

\section{Elementary Considerations}

When the gas is optically thin, the
radiation field at each radius is determined simply by   geometrical
dilution of the given source spectrum $f_\varepsilon$.  
Then, as shown by \cite{Tarter1969},
the state of the gas depends only on the ionization parameter  $\xi=L/nR^2$, 
where $L$ is the (energy) luminosity of the incident radiation integrated from 1 to 1000 Ry,
$n$ is the gas density, and $R$ is the distance from the radiation source.  This
scaling law allows the results of one model   calculation to be
applied to a wide variety of situations.  For a given choice of
spectral shape this parameter   is
proportional to the  various other customary ionization parameter 
definitions, i.e.  $U_H=F_H/n$
(\cite{Davidson1979}), where $F_H$ is the incident photon number 
flux above 1 Ry; $\Gamma=F_\nu(\nu_L)/(2hcn)$, where 
$F_\nu(\nu_L)$ is incident (energy) flux at 1 Ry; and  
$\Xi=L/(4\pi R^2 cnkT)$ (e.g. \cite{Krolik1981}).

In the optically thick case, 
 \cite{Hatchett1976}, and \cite{Kallman1983}
showed that the state of the gas could be  
parameterized in terms of an additional parameter which is a
function of the product of $L$ and either $n$ (the number density) 
or $P$ (the pressure), depending on which quantity is held fixed. 
In the case $n$ = constant, 
this second parameter is simply $(Ln)^{1/2}$ (\cite{McCray1977}).
This parameter does not allow easy scaling of model results from 
value of $Ln$ to another, since the dependence on this parameter 
is non-linear, but it does provide a useful indicator of which 
combinations of parameter values are likely to yield similar results
and vice versa.

When the electron scattering optical depth, $\tau_e$, of the cloud
becomes significant, the outward-only approximation used here breaks
down, and different methods of describing the radiative
transfer must be used (e.g. \cite{Ross1979}). Therefore, the range of 
validity of the models presented here is restricted to $\tau_e\le0.3$,
or electron column densities $\le10^{24}$ cm$^{-2}$.

\section{Algorithm}

The construction of a model consists
of the simultaneous determination of the state of the   gas and the
radiation field as a function of distance from the source. 

\subsection{Atomic Level Populations}

The state of the gas is defined by its temperature and by the ionic level
populations.   As a practical matter, we maintain the distinction 
between the total abundance of a given ion relative to its parent 
element (the ion fraction or fractional abundance) and the relative 
populations of the various bound levels of that ion (level populations), 
although such distinctions are somewhat arbitrary given the presence of 
transitions  linking non-adjacent ion stages and excited levels of
adjacent ions.  

Calculation of level populations proceeds in 2 steps.
First, a calculation of ion fractions is performed using total ionization 
and recombination rates into and out of each ion analogous to those used in XSTAR v.1 (KM82).
Then we eliminate ions with abundance less than a fixed fraction $\epsilon$ relative to 
hydrogen from further consideration.  Experimentation has shown that 
$\epsilon=10^{-8}$ (as parameterized by the input parameter critf) 
yields gas temperatures within 1$\%$ of those calculated 
using a larger set of ions for most situations when the density is low.   
This criterion leads to ion 
sets which can include up to 10 stages for heavy elements such as iron and nickel.
We also make sure that the selected ions are all adjacent. i.e. we force the 
inclusion of ions which fall below our threshold if they are bracketed by 
ions which satisfy the abundance crterion.

The second step consists of solving the full kinetic equation matrix linking the 
various levels of the ions selected in step 1.  We include all processes in the 
database which link the bound levels of any ion in our selected set with any 
other level, and also including the bare nucleus as the continuum level for the 
hydrogenic ion, if indicated.  This results in a matrix with dimensionality
which may be as large as 2400.  The equations may be  written
schematically as $(rate\ in) = (rate\ out)$ for each level.  In place of the 
equation for the ground level of the most abundant ion we solve the number conservation 
constraint.

We include collisional and radiative bound-bound transitions
(with continuum photoexcitation), collisional
ionization, photoionization and  recombination for
all the levels of every ion for which the required atomic rate data is available.  
The  effects of line scattering in all transitions are accounted for 
by taking into account the fact that 
line scattering reduces the net decay rate by repeated 
absorption and reemission of the line photon.  
An analogous procedure is used for free-bound 
(recombination) transitions.

\subsection{Atomic Levels}

A large fraction of recombinations occur following cascades 
from a very large number of levels close to the continuum.  Since 
explicit treatment of these levels is not feasible, we treat this process 
as follows (this procedure, along with detailed descriptions of 
other aspects of the database and 
the multilevel scheme are described in detail in \cite{bautista2000}):  
For every ion we choose a set of
spectroscopic levels starting with the ground level, which are responsible 
for the identifiable emission lines and recombination continua; there are 
typically 10 -- 50 of these for most ions, although for a few ions we include 
$\geq$100 such levels.  In addition
we include one or more superlevels and continuum levels. The continuum levels represent
bound levels of more highly ionized species (in practice at most only a few such levels 
are of importance).  The  superlevel is an artificial level used to account for 
recombination onto the infinite series of levels that lie above the spectroscopic levels. 
In H and He-like ions the superlevels also account for the 
recombination cascades of these high lying levels onto the spectroscopic
levels, and the rates for such decays are calculated by fitting to the results of 
population kinetic calculations for individual ions which explicitly include $\geq$1000 levels.
For these isoelectronic sequences we explicitly include 
excited levels with a spectator electron, which give rise to satellite lines, 
excitation of these levels accounts for excitation-autoionization and radiative deexcitation, 
and recombination accounts for the dielectronic recombination
process.  For other iso-electronic sequences, the superlevels are 
assumed to decay directly to the ion's ground level, and the rates into and out of the
superlevel are calculated in order to fit to the total recombination rates 
for the various ions \cite{bautista2000}.  This approach   
allows us to simultaneously account for the contributions of excitation, ionization, and 
recombination to the ion's level populations.   In this way we solve 
ionization and excitation balance without the use of total recombination rates
which is customary in many nebular calculations.

By using the  approach described above and providing that every
transition process accompanied by its detailed balance inverse process 
we insure that the level populations will 
naturally converge to LTE under proper conditions. 


\subsection{Thermal Equilibrium}


The temperature is found by solving the equation of thermal equilibrium,
which may be  written schematically as $(Heating) = (Cooling).$
This is solved simultaneously with the condition of 
charge conservation.  
We treat heating and cooling by calculating the rate of removal or addition 
of energy to local radiation field associated with each of the processes 
affecting level populations (this is in contrast to the method
 where these were calculated via their effects on the 
electron thermal bath as in KM82).  Heating therefore includes photoionization heating 
and Compton heating.  The cooling term includes radiative
recombination, bremsstrahlung, and radiative deexcitation of bound levels.
Cooling due to recombination and radiative deexcitation is included 
only for the escaping fraction, as described elsewhere in this section.

In  the most highly ionized regions of our models, the dominant heating
process is electron  recoil following Compton scattering.  In the non-relativistic 
approximation the net heating rate may be written (Ross 1979)

$$n_e\Gamma_e={{\sigma_T}\over{m_ec^2}}\left(\int{\varepsilon 
J_\varepsilon d\varepsilon}
-4kT\int{J_\varepsilon d\varepsilon}\right) \eqno{(1)} $$

Here $\sigma_T$ is the Thomson
cross section, $n_e$  is the electron number density, T is the electron
temperature, and $J_\varepsilon$  is the local 
mean intensity in the radiation
field.  The first term in the brackets  represents the heating of
electrons by the X-rays, and the second term represents cooling of 
hot electrons  by scattering with low energy photons.
New in version 2.1e is a full relativistic treatment of Compton heating 
and cooling, and the addition of photon energy bins extending to 1 MeV.
The energy shift per scattering is calculated by interpolating in a
table (coheat.dat) which was calculated using the reults given by Guilbert (1986).

 The spectrum of photoelectron  energies for each ion is
found by convolving the radiation field, weighted by photoelectron
energy,  with the photoionization cross section (see, e.g., 
\cite{Osterbrock1974}).  The integral of this quantity provides the 
photoelectric heating rate.

The cooling rate due
to radiative recombination is calculated by explicitly evaluating the 
quadrature over the recombination continuum spectrum for 
each recombining level, weighted by the escape fraction for that 
transition.  The bremsstrahlung cooling
 rate is (\cite{Osterbrock1974})

$$n_e\Gamma_e=1.42\times 10^{-27} T^{1/2} z^2 n_e n_z 
{\rm ~ergs~cm^{-3} s^{-1}}, \eqno{(2)}$$
where T is the electron temperature,   
is the electron number density, z is the charge on the cooling  ion,
and $n_z$ is the ion density.


\section{Recombination Continuum Emission and Escape}

In analogy with the line emission, recombination emission and cooling rates
are calculated using the continuum level population $n_{\infty}$ and 
the quantities calculated from the photoionization cross section 
and the Milne relation. The spontaneous recombination rates are given by

$$\alpha_i=({{n_i}\over{n_{i+1} n_e}})^* 
\int_{\varepsilon_{th}}^{\infty}{{d\varepsilon}\over{\varepsilon}}
{{\varepsilon^3}\over{h^3 c^2}} \sigma_{pi} 
e^{(\varepsilon_{th}-\varepsilon)/kT}  \eqno{(3)} $$

\noindent where ${n_i}^*$ is the LTE 
density of ion $i$, and $n_e$ is the electron density.
The continuum emissivity due to this process is given by

$$j_\varepsilon=n_{upper} n_e ({{n_i}\over{n_{i+1} n_e}})^* 
{{\varepsilon^3}\over{h^3 c^2}} \sigma_{pi} 
e^{(\varepsilon_{th}-\varepsilon)/kT}  \eqno{(4)} $$

\noindent where $n_{upper}$ is the number density of ions in 
the recombining level  and $\sigma_{pi}$ is the photoionization cross 
section.  The cooling rate is given by the integral of this expression over
energy.  These rates are calculated separately for each level
included in the multilevel calculation.

In order to account for the suppression of rates due to emission and 
reabsorption of recombination continua, we multiply the rates and 
emissivities by an escape fraction given by:

$$P_{esc., cont.}={{1}\over{1000 \tau_{cont.}+1}}  \eqno{(5)} $$

\noindent where $\tau_{cont.}$  is the optical depth at the threshold
energy for  the relevant transition.  This factor is used to correct 
both the emission rate for the recombination events, and also the 
rates in the kinetic equations determining level populations etc, and has 
been found to give reasonably good fits to the results of more 
detailed  calculations for the case of H II region models 
in which the Lyman continuum of hydrogen is optically thick
(e.g. \cite{Harrington1989}).


\subsection{Line Emission and Escape}


Since all level populations are calculated explicitly, line emissivities
and cooling rates are calculated as a straightforward product of the population of the 
line upper level, the spontaneous transition probability and 
an escape fraction.

Line optical depths may be large in some nebular situations. 
Photons emitted near the centers of these lines are likely to be
absorbed by the  transition which emitted them and reemitted at a new
frequency.  This line scattering will repeat  many times until the
photon either escapes the gas, is destroyed by continuum photoabsorption or 
collisional deexcitation, or is degraded into longer wavelength  photons 
which may then escape. Our treatment of resonance line transfer is based on the
assumption of complete  redistribution.  That is, we assume that there
is no correlation of photon frequencies before and  after each
scattering event. This has been shown to be a good approximation for a
wide variety of  situations, particularly when the line profile is
dominated by Doppler broadening.  In this case, more accurate numerical
simulations (e.g., \cite{Hummer1971}) have shown that line scattering
is  restricted to a small spatial region near the point where the
photons are emitted.  Line photons first  scatter to a frequency such
that the gas cloud is optically thin and then escape in a single long 
flight.  The probability of escape per scattering depends on the
optical depth, $\tau_0$  at the center of the  line.  For       
$1\leq\tau_0\leq 10^6$, the resonant trapping is effectively local.  For  
$\tau_0\geq 10^6$, the  lines become 
optically thick in the damping wings, and the
line escapes as a result of diffusion in  both space and frequency. 
Since the scattering in the Doppler core is always dominated by 
complete redistribution, and since most of the lines in our models are
optically thin in the wings,  we assume that all line scattering takes
place in the emission region.   

We use the following expression for escape probability (\cite{Kwan1981}):

$$P_{esc., line}(\tau_{line})={{1}\over{\tau_{line}\sqrt{\pi}(1.2+b)}} (\tau_{line}\geq 1)  \eqno{(6)} $$

$$P_{esc., line}(\tau_{line})={{1-e^{-2 \tau_{line}}}\over{2 \tau_{line}}} (\tau_{line}\leq 1)  \eqno{(7)} $$

\noindent where

$$b={{\sqrt{{\rm log}(\tau_{line})}}\over{1+\tau_{line}/\tau_w}}  \eqno{(8)} $$,

\noindent $\tau_{line}$  is the optical depth at line center, 
and $\tau_w=10^5$. 

The rates for line emission and the
probabilities for the various resonance line escape and  destruction
probabilities depend on the state of the gas at each point in the
cloud.  The cooling  function for the gas depends on the line escape
probabilities, and the effects of line trapping must  be incorporated
in the solution for the temperature and ionization of the gas. 
Once the state of the  gas at a given point has been determined, the emission
in each  line is calculated as the product of the upper level population
and the corresponding net decay rate, including the suppression due to 
multiple scattering.

\subsection{Continuum Emission}

Diffuse continuum
radiation is emitted by three processes:  thermal bremsstrahlung, 
radiative recombination, and two-photon decays of metastable levels. 
The thermal bremsstrahlung  emissivity is
given by \cite{Osterbrock1974}:

$$j_\varepsilon={{1}\over{4\pi}}n_zn_e{{32Z^2e^4h}\over{3m^2c^3}}
\left({{\pi h \nu_0}\over{3kT}}\right)^{1/2}e^{-h\nu/kT}g_{ff}(T,Z,r)  \eqno{(9)} $$

where T is the electron temperature,  $n_e$
 is the electron abundance, Z is the charge on the most  abundant ion,
$n_z$ is the abundance of that ion, and $g_{ff}$ is a Gaunt
factor (\cite{Karzas1966}). For two photon decays, we adopt the
distribution (\cite{Tucker1971}):

$$H\left({{\varepsilon}\over{\varepsilon_0}}\right)=
12\left({{\varepsilon}\over{\varepsilon_0}}\right)^2
\left(1-{{\varepsilon}\over{\varepsilon_0}}\right)  \eqno{(10)} $$
where $\varepsilon_0$ is the excitation
energy. 

\subsection{Continuum Transfer}

The continuum
radiation field is modified primarily by photoabsorption, for which
the  opacity, $\kappa(\varepsilon)$, 
is equal to the product of the ion abundance
with the total photoionization cross  section, summed over all levels. 

A model is constructed by dividing the cloud into a set of
concentric spherical shells.  The  radiation field incident on the
innermost shell is the source spectrum.  For each shell, starting with
the innermost one, the ionization and temperature structure is
calculated from the local balance  equations using the radiation field
incident on the inner surface.  The attenuation of the incident radiation field 
by the shell is then calculated.   The diffuse radiation emitted by the cloud 
is calculated using an expression of the formal solution if the equation of 
transfer:

$$L_\varepsilon=\int_{R_{inner}}^{R_{outer}}{4\pi R^2 j_\varepsilon(R)
e^{-\tau_{cont.}(R,\varepsilon)}dR}   \eqno{(11)} $$

\noindent where  $L_\varepsilon$ is the specific luminosity at the 
cloud boundary, $\tau_{cont.}(R,\varepsilon)$ is the optical depth
from $R$ to the boundary, and $j_\varepsilon$ is the emissivity at
the radius $R$.  Since our models in general have two boundaries, 
we perform this calculation for radiation escaping at both the 
inner and outer cloud boundaries.   This calculation is performed 
for each continuum energy bin, and 
separately for each line.  In the case of the continuum, we construct a 
vector of emissivities, $j_\varepsilon(R)$ which includes contributions 
from the escaping fraction  from all the levels which affect each 
energy.  For the lines, the emissivity used in this equation is the
escaping fraction for that line.

\subsection{Radiation Field Quantities and Transfer Details}

Equation (11) conceals a variety of important issues concerning the treatment of 
the radiation field and the values which are printed in the various output files
produced by xstar.  In an effort to clarify this we present here 
a complete description of the various radiation field quantities which are used
internally to xstar, and which are output to the user.  In this subsection,
all radiation fields are specific luminosity, $L_\varepsilon$, in units 
erg/s/erg for the continuum,  and luminosity, $L_i$, in units erg/s for lines.
We distiguish several different radiation fields.  First, the radiation 
field used locally by xstar for the calculation of photoionization rates and 
heating, we denote $L_{\varepsilon}^{(1)}$.  This is calculated during an outward
iteration using the transfer equation:

$$\frac{{\rm d}L_{\varepsilon}^{(1)}}{{\rm d}R}=-\kappa_{cont}(\varepsilon) L_{\varepsilon}^{(1)}
+ 4 \pi R^2 j_{\varepsilon}(R)   \eqno{(12)} $$

\noindent with the boundary condition that $L_{\varepsilon}^{(1)}=L_{\varepsilon}^{(inc)}$ at the 
inner radius of the cloud.  Here $\kappa_{cont}(\varepsilon)$ and $j_{\varepsilon}$ are the local continuum opacity 
and emissivity  and $L_{\varepsilon}^{(inc)}$ is the incident 
radiation field at the inner edge of the cloud.  

In addition we can define the various radiation fields of interest for use in fitting to 
observed data.  These include the spectrum transmitted by a model, i.e. the radiation which 
would be observed if the incident radiation field were subject to absorption alone:

$$L_\varepsilon^{(2)}=L_{\varepsilon}^{(inc)} e^{-\tau_{cont}^{(tot)}(\varepsilon)}   \eqno{(13)} $$

\noindent where $\tau_{cont}^{(tot)}(\varepsilon)$ is the total optical depth through the model 
cloud due to continuum photoabsorption,

$$\tau_{cont}^{(tot)}(\varepsilon)=\int_{R_{inner}}^{R_{outer}} \kappa_{cont}(\varepsilon) dR \eqno{(14)}$$

\noindent Also of interest is the total emitted continuum radiation in both the inward and outward
directions, which is given by equations similar to (11):

$$L_\varepsilon^{(3)}=\int_{R_{inner}}^{R_{outer}}{4\pi R^2 j_\varepsilon(R)
e^{-\tau_{cont}^{(in)}(\varepsilon)}
P_{esc, cont.}^{(in)}(R)}dR   \eqno{(15)} $$

$$L_\varepsilon^{(4)}=\int_{R_{inner}}^{R_{outer}}{4\pi R^2 j_\varepsilon(R)
e^{-\tau_{cont}^{(out)}(\varepsilon)}
P_{esc, cont.}^{(out)}(R)dR}   \eqno{(16)} $$

\noindent where the escape probabilities in the inward and outward directions
 are  $P_{esc, cont.}^{(in)}(R)=(1-C)/2$ and $P_{esc, cont.}^{(out)}(R)=(1+C)/2$, where
$C$ is the covering fraction, specified as an input parameter,
and $\tau_{cont}^{(out)}(\varepsilon)$ and $\tau_{cont}^{(out)}(\varepsilon)$ are the continuum
optical depths in the inward and outward directions.


Line luminosities are calculated separately, one at a time, according to an equation 
analogous to equation (12):

$$\frac{d L_i^{(1)}}{dR}=-\kappa_{cont}(\varepsilon) L_i^{(1)}
+ 4 \pi R^2 j_i(R) P_{esc,line}^{(in)}   \eqno{(16)} $$

$$\frac{d L_i^{(2)}}{dR}=-\kappa_{cont}(\varepsilon) L_i^{(2)}
+ 4 \pi R^2 j_i(R) P_{esc,line}^{(out)}   \eqno{(17)} $$

\noindent where $L_i^{(1)}$ and $L_i^{(2)}$ are the luminosities of individual lines 
in the inward and outward directions, respectively.  The escape probabilities in the inward 
and outward directions are calculated using $P_{esc., line}(\tau_{line})$ from equations (6)-(8) and 
$P_{esc., line}^{(in)}=(1-C) P_{esc., line}(\tau_{i}^{(in)})$ and 
$P_{esc, line}^{(out)}=(1-C) P_{esc., line}(\tau_{i}^{(out)})
+ C  P_{esc., line}(\tau_{i}^{(out)}+\tau_{i}^{(in)})/2$, and $\tau_{i}^{(in)}$ and $\tau_{i}^{(out)}$
are the line scattering optical depths in the inward and outward directions:

$$\tau_{i}^{(in)}(R)=\int_{R_{inner}}^R \kappa_i {\rm dR} \eqno{(18)}$$

$$\tau_{i}^{(out)}(R)=\int_R^{R_{outer}} \kappa_i {\rm dR} \eqno{(19)}$$.

\noindent and $\kappa_i$ is the line center opacity.

None of the continuum luminosities defined in equations (12)-(16) have the effects of lines 
included, either in emission or absorption.  This is because lines scatter 
the radiation, while photoionization is true absorption.  The effects of lines on the 
continuum can be added to the continuum for the purposes of comparing with observed 
spectra by binning the lines, i.e. we can calculate  the binned specific luminosity and 
opacity:

$$L_{line, \varepsilon}^{(in)}=\Sigma_{i \ni \mid\varepsilon_i -\varepsilon\mid \leq \Delta\varepsilon}
\frac{L_i^{(1)}  \phi(\varepsilon-\varepsilon_i)}{\Delta \varepsilon} \eqno(20)$$

$$L_{line, \varepsilon}^{(out)}=\Sigma_{i \ni \mid\varepsilon_i -\varepsilon\mid \leq \Delta\varepsilon}
\frac{L_i^{(2)} \phi(\varepsilon-\varepsilon_i) }{\Delta \varepsilon} \eqno(21)$$

$$\kappa_{line}(\varepsilon)=\Sigma_{i \ni \mid\varepsilon_i -\varepsilon\mid \leq \Delta\varepsilon}
\kappa_i \phi(\varepsilon-\varepsilon_i) \eqno(22)$$

\noindent where $\varepsilon$ and $\Delta\varepsilon$ are  the energy and width, respectively, of the 
continuum bin closest to line $i$, and $\phi(\varepsilon-\varepsilon_i)$ is the profile function 
including the effects of broadening due to thermal Doppler motions, natural broadening, and turbulence.

Then we can define the total optical depth of the cloud

$$\tau^{(tot)}(\varepsilon)=\int_{R_{inner}}^{R_{outer}} 
(\kappa_{cont}(\varepsilon)+\kappa_{line}(\varepsilon)) {\rm dR} \eqno{(23)}$$

\noindent and the total transmitted specific luminosity

$$L_\varepsilon^{(5)}=L_{\varepsilon}^{(inc)} e^{-\tau^{(tot)}(\varepsilon)} \eqno{(24)}  $$

\noindent and the total emitted specific luminosity in the inward and outward directions:

$$L_\varepsilon^{(6)}=L_\varepsilon^{(3)}+L_{line, \varepsilon}^{(in)}  \eqno{(25)} $$

$$L_\varepsilon^{(7)}=L_\varepsilon^{(4)}+L_{line, \varepsilon}^{(out)}   \eqno{(26)} $$

The quantities $L_\varepsilon^{(inc)}$ 
$L_\varepsilon^{(5)}$, $L_\varepsilon^{(6)}$
and $L_\varepsilon^{(7)}$ are output in columns 2,3,4,5 of the file xout\_spect1.fits.
The quantities $L_\varepsilon^{(inc)}$ 
$L_\varepsilon^{(5)}$, $L_\varepsilon^{(3)}$
and $L_\varepsilon^{(4)}$ are output in columns 2,3,4,5 of the file xout\_cont1.fits.

In fact, the lines should be included in the continuum which is responsible for 
the local ionization and heating of the gas, since they can contribute to these processes.
So we define a modified version of equation (12):

$$\frac{{\rm d}L_{\varepsilon}^{(1')}}{{\rm d}R}=-\kappa_{cont}(\varepsilon) L_{\varepsilon}^{(1')}
+ 4 \pi R^2 j_{\varepsilon}(R)  
+ 4 \pi R^2 \Sigma_{i \ni \mid\varepsilon_i -\varepsilon\mid \leq \Delta\varepsilon}
\frac{j_i(R) P_{esc,line}^{(out)} \phi(\varepsilon-\varepsilon_i) }{\Delta \varepsilon}  \eqno{(27)} $$

\noindent $L_{\varepsilon}^{(1')}$ is the quantity which is used by xstar to calculate the local
ionizing flux.  This is the quantity which is conserved by xstar when it calculates 
heating=cooling.

The quantities $L_i^{(1)}$, $L_i^{(2)}$, $\tau_i^{(in)}$ and $\tau_i^{(out)}$ 
are output in columns 6,7,8,9 of the file xout\_lines1.fits.

The quantities which contain the continuum only, before the lines are binned and added, 
are printed out to the log file, xout\_step.log, when the print switch lpri is set to 1 or greater.
Then they are in a table following the label 'continuum luminosities'.  
The quantities $L_\varepsilon^{(inc)}$, $L_\varepsilon^{(1')}$, $L_\varepsilon^{(3)}$, 
$L_\varepsilon^{(4)}$, $\tau_{cont}^{(in)}(\varepsilon)$ and $\tau_{cont}^{(out)}(\varepsilon)$ 
are output in columns 3,4,5,6,7 and 8.
Many other useful quantities are output to the log file when lpri=1.  This includes the 
quantities $L_i^{(1)}$, $L_i^{(2)}$, $\tau_i^{(in)}$ and $\tau_i^{(out)}$ 
are in columns 2-5 following the label 'line luminosities'

\subsection{Energy Conservation}

Energy conservation is imposed as a constraint when determining the temperature in xstar when the 
input parameter niter is non-zero.  If so, the temperature is iteratively improved until 
the heating and cooling rates are locally equal.  This is implemented by calculating the 
integral over the absorbed and emitted continuum energy in a given spatial zone, and also 
the sum over the energy emitted in the lines.  Compton heating and cooling are added 
analytically, since Comptonization of the radiation field is not treated. 
The error resulting from this procedure is tabulated in the log file 'xout\_step.log' in the 8th column
of the step-by-step output, in units of $\%$.  

Energy conservation locally should correspond to global energy 
conservation, i.e. that the total absorbed energy in the radiation field 
equals the total emitted energy in lines plus continuum.  This is 
tested at each spatial zone in xstar by calculating 
$\int{(L_{\varepsilon}^{(inc)} - L_{\varepsilon}^{(1)})}d\varepsilon - \Sigma_i(L_i^{(1)}+L_i^{(2)})$.
The error resulting from this procedure is tabulated in the log file 'xout\_step.log' in the 9th column
of the step-by-step output, in units of $\%$.  

It is important to point out that the specific luminosities in 
the file 'xout\_spect1.fits' 'xout\_cont1.fits' are not expected, in general, 
to show energy conservation.  This is primarily because the transmitted 
spectra in both of these files contain the effects of binned lines.
Line opacity is expected to produce a scattering event, i.e. the photon is 
likely to be reemitted near the same energy.  This  differs qualitatively from
photoelectric absorption, in which an absorbed photon is likely to be reemitted 
at a very different energy, with an accompanying net loss or gain of energy 
to the electron thermal bath.  Line opacity is not included in the radiative 
equilibrium integral used to calculate the gas temperature, and so 
the total absorbed energy in the radiation field  $L_{\varepsilon}^{(5)}$ will in 
general not equal the emitted energy in $L_{\varepsilon}^{(6)}+L_{\varepsilon}^{(7)}$.
Energy conservation can be checked using the quantities 
$L_{\varepsilon}^{(1)}$, $L_i^{(1)}$ and $L_i^{(2)}$ 
from the file xout\_step.log. 

Energy conservation checked using binned line spectra is also affected by the 
errors introduced by binning.  This is discussed at length in the section of this 
manual on table models for xspec, but we emphasize here that the binned 
spectrum cannot be accurately integrated to derive the total 
line absorption or emission unless the lines are broad compared 
with the energy grid spacing (requiring turbulent velocities $\geq$ 500 km/s currently).


\subsection{Algorithm}

Construction of a model of an X-ray
illuminated cloud consists of the simultaneous  solution of the local
balance equations.  The radiative transfer equation is
solved for both  the continuum and for the lines that escape the
region near the point of emission.  The large number  of ions in the
calculation results in many ionization edges that may affect the
radiation field.  We  solve the transfer equation on a frequency grid
that includes a total of 9999 continuum  grid points with even 
logarithmic spacing in energy from 0.1 eV to 20 keV resulting 
in a limiting resolution of  0.12 $\%$, corresponding to, e.g. 8.6 eV at 7 keV.
We calculate the luminosities of $\sim$10000 spectral lines and
solve the continuum transfer equation  individually for each of these.
The emissivity of each line at each point is the product of the 
emissivity  and the local escape fraction for
that line. The continuum opacity for each
line is the opacity calculated for the energy bin that contains the 
line. This procedure  is repeated for each successive shell with
increasing radius.

Calculation of the escape of the diffuse radiation 
field depends on a knowledge of the 
optical depths of the cloud from any point to both the inner and outer
boundaries.  Since these are not known a priori we iteratively 
calculate the cloud structure by stepping through the radial shells
at least 3 times.  For the initial pass through the shells we assume 
that the optical depths in the outward direction are zero.  This 
procedure is found to converge satisfactorily within 3-5 passes 
for most problems of interest.  This procedure is tantamount to the 
``$\Lambda$-iteration'' procedure familiar from stellar atmospheres, and
must suffer from the same convergence problems when applied to problems
with large optical depths.  These problems are reduced in our case 
by the use of escape probabilities rather than a full integration 
of the equation of transfer.

\section{Atomic Processes}

Here we summarize the most important data sources adopted for the 
calculations.  These are discussed in greater detail, along with 
a description of the fitting formulas and assumptions, in 
\cite{bautista2000}.

\subsection{Photoionization}

Photoionization rates are
obtained by convolving the radiation field with the photoionization
cross  section.  Cross sections are included for all levels of every 
ion for a 
wide range of photon  energies occurring in our model. 
The cross sections where taken from the Opacity Project (\cite{opacityp},
\cite{topbase}), then averaged over resonances as in
\cite{bautista1997} and split over fine structure according to
statistical weights.  

For inner shell photoionization not yet available from the 
opacity project we use the cross sections of 
\cite{Verner1995}.  Inner shell ionization in X-ray illuminated clouds
is enhanced by Auger cascades.  This process can result in the ejection of up to eight
extra  electrons (in the case of iron) in addition to the original
photoelectron.  We include this effect by treating each inner shell ionization/auger event
as a rate connecting the ground state of one ion with another level
of an ion (in general not adjacent to the initial ion).  The rates for 
inner shell ionization/auger processes are calculated using the relative 
probabilities of the various possible outcomes of 
an inner shell ionization event from \cite{Kaastra1989}.  These yields are 
multiplied by the appropriate inner shell photoionization cross section 
in order to calculate a rate for each inner shell ionization/Auger cascade 
individually.  Our level 
scheme also includes levels with inner shell vacancies, which are populated 
by inner shell/Auger events.  Populations of these levels are calculated 
in the customary way, and the decays from these levels produce inner shell 
flourescence lines.

\subsection{Collisional Ionization}

Ionization by electron collisions is important
if the gas temperature approaches a fraction of the ionization threshold energy of
the most abundant ions in the gas.   For ground states we include 
the rates from \cite{Raymond1976} for elements other than 
iron, and from \cite{Arnaud1992} for iron. 
Collisional ionization from excited levels may also be
important to the ionization balance. We include ionization rates 
for all excited levels of every ion using approximate formulae by 
D. Sampson and coworkers (\cite{sampsonzhang}).   
3-body recombination rates to all levels are calculated from  the
collisional ionization rates using the detailed balance principle. 

\subsection{Recombination}

Radiative and dielectronic
recombination rates to all spectroscopic levels are calculated from the
photoionization cross sections using the Milne
relation.   We include both spontaneous and stimulated recombination 
caused by the illuminating radiation.  Stimulated recombination by the 
locally emitted radiation is not treated explicitly, although its effect 
is taken into account in an approximate way by suppressing a fraction of the 
spontaneous recombinations using the escape probability described 
earlier in this section.
Recombination onto the superlevels is calculated in order to account
for the difference between the sum over all spectroscopic levels and 
the total ion recombination as given by Nahar and coworkers (\cite{Nahar1999},
\cite{Nahar2000}) where available and \cite{Aldrovandi1973} for species 
other than iron ions and ions in the H and He isoelectronic sequences.
For iron we use total rates from 
\cite{Arnaud1992}. For H and He-like ions the total recombination rates
were calculated by \cite{bautista1998} and \cite{bautista2000}. 

\subsection{Collisional Excitation and Radiative Transition Probabilities}

Collision strengths and A-values were collected from a large number
of sources.  Particularly important for this compilation were the 
CHIANTI data base (\cite{chianti}) for X-ray and EUV lines,  
and the extensive R-matrix calculations by the Iron Project (\cite{ironp}).

\subsection{Charge Transfer}

Rates for charge transfer reactions are taken from \cite{Butler1980}.
For highly charged ions, where accurate calculations
do not exist, we scale the rates  along isonuclear sequences,
assuming that the cross section is proportional to the square of the 
total residual charge transfer reaction of O II with H 
(\cite{Field1971}).

