The Physics Behind XSTAR

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.

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 \(\text{erg}\,\text{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 Grevesse et al. [101].

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 Tarter et al. [187], 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\) [74], 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. Krolik et al. [126]).

In the optically thick case, Hatchett et al. [106], and Kallman [113] 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}\) [141]. 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. Ross et al. [172]). Therefore, the range of validity of the models presented here is restricted to \(\tau_e\le0.3\), or electron column densities \(\le10^{24}\,\text{cm}^{-2}\).

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.

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 (ratein) = (rateout) 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.

Radiative Excitation

Excitation due to absorption of continuum photons in bound-bound transitions can dominate over other mechanisms, depending on the conditions. The importance of this has been pointed out and explored by Kinkhabwala et al. [124] and Band et al. [25]. Accurate treatment of this process requires that the geometry be treated carefully. This is because the net excitation rate is also affected by the radiation emitted by the decays of the excited levels. If every radiative excitation leads to a radiative deexcitation and emission of radiation which escapes the model volune then there is no net effect of this process. In the case of a stationary spherical model in which we are only interested in the total spectrum seen by a distant observer, and if radiative deexcitation is the only decay mechanism for the excited levels, then emission following radiative excitation and absorption will balance exactly. This is the scenario envisioned by classical nebular models, and for this reason these models did not generally include the effects of radiative excitation.

Radiative excitation has a cross section which can be much greater than the photoabsorption cross sections associated with photoionization. If so, it is important to accurately treat the spatial dependence of the absorption and depletion of the incident radiation in lines. This necessitates spatial steps which are small enough to resolve the line absorption, i.e. \(\sim \kappa^{-1}\), where \(\kappa\) is the opacity in the line. Xstar accounts for this, but using the absorption cross section binned into the continuum grid. If the grid is coarse, then there will likely be no bin close to the centers of the strongest lines, and this process will not be accurately modeled. For this reason it is recommended that a fine energy grid be used for models where radiative excitation is of interest, and where the model cloud is optically thick in the line centers.

In order to preserve the classical result in which radiative excitation is completely ignored, xstar scales the rate for this process proportional to the value of the cfrac parameter. The rate of radiative excitation of a given transition is calculated using:

\[R_{lu}=\sigma_{line}F_\varepsilon\]

where :math:sigma_text{line}` is the line mean photoabsorption cross section and \(F_\varepsilon\) is the specific flux at the energy of the line. This rate is then multiplied by a factor of 1-cfrac. So, if cfrac=1 the effective rate of radiative excitation is zero. If cfrac=0 then the full rate of radiative excitaiton is included.

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 : 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 \(\geq100\) 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 \(\geq1000\) 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 []. 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.

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) (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. The treatment of Compton heating and cooling in versions prior to 2.3 were not accurate for hard spectra with significant flux above 100 keV. This has been updated in version 2.3 using rates from I. Khabibullin (private communication), based on the expressions given by Shestakov et al. (1988 JQSRT 40 577) The energy shift per scattering is calculated by interpolating in a table (coheat.dat).

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., Osterbrock [160]). 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 (Osterbrock [160])

\[n_e\Gamma_e=1.42\times 10^{-27} T^{1/2} z^2 n_e n_z {\rm ~ergs~cm^{-3} s^{-1}}, (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.

Heating and Cooling: “Photon Point of View” vs. “Electron Point of View”

Consideration of heating and cooling depends on the definition of the energy scale. Two different natural ways to do this are the Photon Point of View and the Electron Point of View. In the former description we keep track of the destruction and creation of photons, and a heating or cooling event corresponds to the these two processes, respectively. In the Electron Point of View we keep track of energy deposited or extracted from the electron thermal bath. The two descriptions differ owing to the internal energy associated with ionization: processes which involoving bound-bound transitions have no effect on the electron thermal bath but do emit or absorb photons. Electron impact collisions (by themselves) have no effect on the photon distribution but do heat or cool the electrons. Photoionization heats the electrons with a rate proportional to the energy of the absorbed photon above threshold; it removes photons at a rate proportional to the energy of the absorbed photons above the ion ground level. Radiative recombination differs in a corresponding way between the two descriptions. Compton heating/cooling and bremsstrahlung cooling are the same in the two descriptions.

A key point is that in the ‘photon point of view’ calculation there is recombination cooling even at low temperature, since every recombination releases a continuum photon and also possibly line photons. Recombination is a monotonically decreasing function of temperature, and the cooling curves have this behavior at low temperature. In constrast, the ‘electron point of view’ cooling rates are dominated by bremsstrahlung, which goes as \(T^{1/2}\), at low and high temperatures.

It is easy to show that the two descriptions are equivalent when ionization equilibrium is achieved; the quantity heating-cooling is the same. Xstar uses the photon point of view for all its calculations of thermal balance. Traditional photoionization codes use the electron point of view. For ease of comparison, xstar (starting with version 2.54) also computes the rates per ion and total rates in the electron point of view and prints them out in the ascii file.

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} (3)\]

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} (4)\]

where \(n_\text{upper}\) is the number density of ions in the recombining level and \(\sigma_\text{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}} (5)\]

where \(\tau_\text{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. Harrington [105]).

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., Hummer and Rybicki [110]) 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 (Kwan and Krolik [128]):

\[ \begin{align}\begin{aligned}P_{esc., line}(\tau_{line})={{1}\over{\tau_{line}\sqrt{\pi}(1.2+b)}} (\tau_{line}\geq 1) (6)\\P_{esc., line}(\tau_{line})={{1-e^{-2 \tau_{line}}}\over{2 \tau_{line}}} (\tau_{line}\leq 1) (7)\end{aligned}\end{align} \]

where

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

\(\tau_\text{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.

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 Osterbrock [160]:

\[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) (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_\text{ff}\) is a Gaunt factor [114]. For two photon decays, we adopt the distribution [190]:

\[H\left({{\varepsilon}\over{\varepsilon_0}}\right)= 12\left({{\varepsilon}\over{\varepsilon_0}}\right)^2 \left(1-{{\varepsilon}\over{\varepsilon_0}}\right) (10)\]

where \(\varepsilon_0\) is the excitation energy.

Continuum Opacity and Thomson Scattering

Continuum photons can be absorbed by photoionization. This is regarded as a photon destruction process since the converse, radiative recombination, produces an ‘RRC’ spectrum which has a different energy dependnce than the absorption, and so the two are treated separately. Xstar includes the photoabsorption associated with every bound-free transition in the atomic database; there are \(\sim 10^5\) of them.

Photons can scatter via Thomson scattering or resonance scattering in lines. Xstar currently includes resonance scattering in lines in the final spectrum seen by a distant observer (as described in the following section). Thomson scattering does not strongly affect the shape of the spectrum unless the photon energy or temperature are a significant fraction of \(m_e c^2\). Furthermore, similar arguments apply to Thomson scattering as apply to the radiative excitation in section 9.D.2. In the case of a spherical stationary cloud where all the photons from the cloud are observed, Thomson scattering will have no net effect. For this reason, version of xstar prior to 2.54a omitted Thomson scattering. Beginning with verion 2.54a, Thomson scattering is included with the same multiplicative factor 1-cfrac as is used for radiative excitation. Thus, when cfrac=1 the classical results without Thomson scattering are recovered, and when cfrac=0 Thomson scattering is included in the continuum opacity. Of cours, Thomson scattering is unimportant when the cloud column is less than \(\sim 10^{23}\text{cm}^{-2}\). And, when the cloud is optically thick to Thomson scattering, xstar does not adequately treat the effects of multiple scatterings.

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} (11)\]

where \(L_\varepsilon\) is the specific luminosity at the cloud boundary, \(\tau_\text{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.

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) (12)\]

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)} (13)\]

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 (14)\]

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

\[ \begin{align}\begin{aligned}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 (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} (16)\end{aligned}\end{align} \]

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)} (16)\]
\[\frac{d L_i^{(2)}}{dR}=-\kappa_{cont}(\varepsilon) L_i^{(2)} + 4 \pi R^2 j_i(R) P_{esc,line}^{(out)} (17)\]

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:

\[ \begin{align}\begin{aligned}\tau_{i}^{(in)}(R)=\int_{R_{inner}}^R \kappa_i {\rm dR} (18)\\\tau_{i}^{(out)}(R)=\int_R^{R_{outer}} \kappa_i {\rm dR} (19)\end{aligned}\end{align} \]

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:

\[ \begin{align}\begin{aligned} 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} (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} (21)\\\kappa_{line}(\varepsilon)=\Sigma_{i \ni \mid\varepsilon_i -\varepsilon\mid \leq \Delta\varepsilon} \kappa_i \phi(\varepsilon-\varepsilon_i) (22)\end{aligned}\end{align} \]

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} (23)\]

and the total transmitted specific luminosity

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

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

\[ \begin{align}\begin{aligned}L_\varepsilon^{(6)}=L_\varepsilon^{(3)}+L_{line, \varepsilon}^{(in)} (25)\\L_\varepsilon^{(7)}=L_\varepsilon^{(4)}+L_{line, \varepsilon}^{(out)} (26)\end{aligned}\end{align} \]

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} (27)\]

\(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’

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).

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 \(\sim10000\) 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.

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 .

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 [71, 181], then averaged over resonances as in Bautista et al. [33] 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 Verner and Yakovlev [191]. 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 Kaastra and Mewe [112]. 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.

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 Raymond and Smith [170] for elements other than iron, and from Arnaud and Raymond [22] 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 [205]. 3-body recombination rates to all levels are calculated from the collisional ionization rates using the detailed balance principle.

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 [154] where available and Aldrovandi and Pequignot [17] for species other than iron ions and ions in the H and He isoelectronic sequences. For iron we use total rates from Arnaud and Raymond [22]. For H and He-like ions the total recombination rates were calculated by Bautista et al. [33] and .

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 [75] for X-ray and EUV lines, and the extensive R-matrix calculations by the Iron Project [109].

Charge Transfer

Rates for charge transfer reactions are taken from Butler et al. [63]. 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 [88].