XSTAR2XSPEC
Introduction
The previous chapter XSTAR output described the various output files that XSTAR produces for a specific set of parameters. Many users will also want to compare XSTAR results with observational data and in particular fit observed X-ray spectra with XSTAR models, i.e., vary XSTAR parameters until a good match between observed and synthetic spectra is reached. To facilitate this process, XSTAR2XSPEC was developed. XSTAR2XSPEC is a perl script which calls XSTAR multiple times and generates table models from the results of these simulations which can then be utilized for model fitting in the xspec spectral fitting program. Table models have advantages and disadvantages and the user is strongly recommended to read the following sections completely which contain several important caveats.
General Notes on Table Models in XSPEC
As indiviudal XSTAR runs can take anywhere between seconds and hours, calling XSTAR from within xspec to evaluate and optimize models is not feasible in many cases. An easy way around this is evaluating XSTAR on a pre-defined grid of parameters and tabulating the results of the individual calculations. This table model can then be loaded into and xspec will interpolate the model spectrum for any given parameter combination within the grid.
XSPEC has a specific table model format (see OGIP memo 92-009 for details) that allows for very efficient searching and interpolating among the tabulated spectra. XSPEC distinguishes between additve table models (\(\texttt{atable}\)), multiplicative table models (\(\texttt{mtable}\)), and exponential table models (\(\texttt{etable}\)) that behave differently regarding normalization and rebinning. Per default, XSTAR2XPEC is designed to create only \(\texttt{atable}\) and \(\texttt{mtable}\) models. Users who want to use \(\texttt{etable}\) models will need to build them from the existing \(\texttt{mtable}\) (make thread to show this).
The advantage of table models is that the calculation of the table can easily be parallelized so even computationally expensive models can be produced in a reasonable amount of time. While table models are very straightforward to calculate and run very fast, there are a few things to consider when using tabulated models with observational data.
First, the number of tabulated spectra grows exponentially with the number of free parameters and linearly with the range and sampling of each parameter. In practice, this usually limits the number of free parameters to about one to four which is oftern rather small for complex physical models and requires additional assumptions or simplifications. The ranges and stepsizes for the parameter grid has to be choosen appropriately to avoid interpolation artifacts in the final model spectrum. This is especially problematic when varying abundances in absorption models.
Second, properties like spectral resolution or energy ranges need to be defined at the time of creation of the table model and cannot be changed afterwards. Therefore, a model that was calculated for a specific instrument or dataset cannot be easily tranfered to another and often needs to be recalculated from scratch.
In order to reduce table sizes, xspec table models consider two types of free parametes: interpolated and additive parameters. The idea behind additive parameters is that they are suffciently independent of all other parameters so that tabulated spectra do not need to stored for every single parameter combination. Interpolated parameters are treated as one might expect: if the model is represented by a vector \(M_i(x_j)\), corresponding to the model flux at various energies \(\varepsilon_i\) and stored at various values of the free (intepolated) parameter \(x_j\), then for a value of the parameter not on the tabulated grid xspec calculates the model value as
where \(\omega_j\) are suitably chosen weights for, e.g. linear interpolation. This is in contrast to the treatment of additive parameters: for an \(\texttt{mtable}\), xspec needs the value of the model tabulated at only two values of the free additive parameter \(y\), 0 and \(y_\text{max}\). Then the model is calculated for some arbitrary \(y\) as
i.e. it is assumed that the model scales linearly with the value of the additive parameter \(y\). This formalism was developed for emission models, where emissivities might be expected to scale linearly with elemental abundance. The priniple is illustrated below.
Illustration of the concept of interpolated versus additive parameters. In this example \(\log\xi\) and \(N_\text{H}\) are interpolated parameters, the iron abundance \(A_\text{Fe}\) and sulfur abundance \(A_\text{S}\) are additive parameters. At each gridpoint in the \(\log\xi\)-\(N_\text{H}\) plane, three spectra are stored. For a given combination \((\log\xi, N_\text{H}, A_\text{S}, A_\text{Fe},)_i\), the spectra are interploted on the \(\log\xi\)-\(N_\text{H}\) grid and for each additive abundance. This way only a total of 60 spectra need to be stored. If all parameters were interpolated the model would contain 80 tabulated spectra.
Variable Abundances in Table Models
If you are using this approach to calculate absorption models (\(\texttt{mtables}\)) with variable abundances, then you will likely get totally unphysical results unless you note the following. In the case of absorption models, the value of \(M\) which xspec uses is a transmission coefficient, and this does not scale linearly with abundance. Rather, the transmission coefficient is related to the optical depth by:
and the optical depth \(\tau_i(y)\) does scale (approximately)
linearly with abundance. For this reason, xspec has incorporated the
\(\texttt{etable}\) models, in which the model builder supplies
optical depths rather than transmission coefficient, linear
interpolation is used to calculate the optical depth for a given
abundance, and then the transmission coefficient is calculated using
the above expression. If so, users should use the model
etable{xout_etable.fits} command in xspec instead. Entering
parameter, etc., in xspec is the same for \(\texttt{etables}\) as
for \(\texttt{mtables}\).
The problem with \(\texttt{mtables}\) models is that the multiplicative factor per spectral bin is calculated as
where \(M_i^0\) is the ‘zero abundance’ version of the model at that ionization parameter and column, \(x_j\) is the abundance of element \(j\), and \(M_i^j\) is the model calculated with that element set to unity (relative to cosmic). The first problem with this approach is that tranmission does not scale this way with abundance. If you double the abundance of an element, you should double the optical depth. Doubling the transmittivity has both the wrong qualitative and quantitative behavior. But an even bigger problem with this is that the whole thing, \(M_i\), multiplies your continuum. It is almost guaranteed to bigger than 1 if some of the \(x_j\) are non-zero. So this is totally the wrong formulation to use when the abundances vary.
In the case of an \(\texttt{etable}\) the transmittivity is:
where \(E_i^0=-{\rm ln}(M_i^0)\) etc.
This is better, since it predicts the correct dependence on abundance, i.e. that the optical depth scales approximately linearly with abundance.
Below are six figures to illustrate these problems. All use the \(\texttt{mtable}\) or \(\texttt{etable}\) from a table model “grid18” that we will also use for other exmaples. Instructions on how to calculate this table can be found at Running XSTAR2XSPEC. They show the multiplicative transmission coefficient between 0.6 and 0.8 keV for a \(\log(\xi)=1\) and \(N_\text{H}=10^{20}\,\text{cm}^{-2}\).
Transmission coefficients for \(\texttt{mtable}\) and \(\texttt{etable}\) models, once for all metal abundances set to zero (only H and He) and once with oxygen abundance set to 0.5 solar and once to solar. The transmission coefficient should always between 0 and 1, with 1 indicating no absorption and 0 full absorption. By construction of the \(\texttt{mtable}\) model, however, the transmission coefficient becomes larger than 1, because transmission coefficients for additive parameters (like elemental abundances) are summed up. Additionally, the line depth is not calculated correctly. For \(\texttt{etable}\) table models, the optical depth \(\tau\) is summed up over additive parameters and the model returns the total transmission coefficient \(\exp\left(-\tau\right)\).
In addition to these obvious problems, it is also important to remember that what you are trying to simulate is a ‘real’ photoionized gas with some abundance set {\(x_j\)} that fits your data. What you are using to fit to your data are \(M_i^0\), which in the case of xstar are models calculated with only hydrogen and helium, and the set of \(M_i^j\), which are supposed to be calculated with pure element \(j\).
The first problem is that the cooling, and therefore the temperature, depends on the element abundance in the gas, and this is not linear. The temperature couples back into the ionization balance, and that affects the opacity, transmittivity, etc.
The second problem is that, while it is straightforward to calculate a model with no elements heavier than H and He, it is not straightforward to calculate a model with pure iron, or calcium, say. That is because H and He have nice simple cooling properties at low temperature, and their opacity is simple, etc. Pure iron models at best are likely to be very different than H+He+Fe models, or from a cosmic mix. So what I do is make the \(M_i^j\) using H+He+element \(j\) in cosmic ratios. But the problem with this is that then the opacity due to H and He are present in every :math:`M_i^j and \(E_i^j\) , particularly at low ionization parameters. So if you sum over cosmic abundances using an etable you are including the H and He opacity with every one of the :math:`E_i^j. At high ionization parameter (\(\log(\xi)>1\)?) this should not be a problem because H and He are ionized sufficiently that they do not contribute to the opacity. But at low ionization parameter, the model spectra calculated using etables will have significantly more low energy opacity due to the multiple counting of H and He than they would in a real XSTAR model with the corresponding parameters.
The strategy to use at low ionization parameters therefore is to use the \(\texttt{etable}\) to find a reasonably good fit to the spectrum with xspec, and then run xstar2xspec and make a very small table, i.e. just 1 column and 2 ionization parameters (xspec needs tables to have at least 2 entries), with constant abundances, set to the values which come from the fit. Constant abundance grids when used as \(\texttt{mtable}\) do not suffer from any of these problems.
Energy Resolution and Line Widths in Table Models
Users may find that turbulent velocity \(v_\text{turb}\) can make a significant difference in the results of xstar modeling of absorption spectra, and the reason is at least partially due to numerics. In calculating the energy dependent opacity xstar tries to be clever: it calculates the profile function on a very fine temporary internal energy grid and then integrates this quantity over the current energy bins. In each bin, the opacity is calculated such that the equivalent width in that bin is equal to the equivalent width that comes from the fine grid calculation. This provides a fairly accurate opacity distribution, which is better than eg. evaluating the line profile function only on grid boundaries.
Close-up of \(\text{O}\,\textsc{viii}\,\text{Ly}\alpha\) lines in \(\texttt{etable}\) models calculated for different turbulent velocities. The solid lines indicate the energy grid (that the model gets evaluatued on), the dotted lines the model grid (that the model gets calculated on).
Notes on Units and Normalization in Table Models
The problem with creating a flexible tool for modelling emission and absorption is that there have several free parameters affecting real spectra, including: source luminosity, distance, reprocessor column density, ionization parameter, and geometry. By geometry we mean covering fraction around the source, which affects emission, and covering fraction across our line of sight to the source, which affects absorption. With the xstar2xspec tables you should be able to model a wide range of choices for this, but there is not a unique one-to-one mapping between the values of these physical parameters and the values used in running xstar and constructing the tables.
The free parameters which can be varied when running xstar2xspec include the abundances, column density, gas density, and ionization parameter. These all have a straightforward intepretation as physical parameters. The emitter normalization and geometry are not uniquely determined, owing to the ambiguity between source luminosity and distance.
It is helpful here to be very specific, at the risk of being repetitive. The procedure followed by xstar2xspec in making tables is:
Generate a sequence of command line calls to XSTAR (this is done by the tool xstinitable)
Step through the calls to XSTAR
Calculate each appropriate XSTAR model
Take the spectrum output of the model (the file xout_spect1.fits) and convert it to the right units and append it to the fits table (xout_aout.fits, xout_ain.fits, or xout_mtable.fits). The last step is done by the tool xstar2table.
The conversion is as follows: For additive models (\(\texttt{atable}\)), xspec expects a bin-integrated spectrum in units of \(\text{photons}\,\text{s}^{-1}\,\text{cm}^{-2}\). We denote this quantity \(F_n^\text{mod}\) for bin \(n\). In XSTAR, however, the user specifies a bolometric source luminosity \(L_\text{tot}^\text{xstar}\) in units \(10^{38}\,\text{erg}\,\text{s}^{-1}\) and XSTAR returns the specific luminosity \(\mathcal{L}_\varepsilon^\text{xstar}\) in units of \(10^{38}\,\text{erg}\,\text{s}^{-1}\,\text{erg}^{-1}\) at each energy \(\varepsilon\) of the model grid. At distance \(D\) from the source, the observed flux in bin \(n\) is therefore
where \(\Delta\varepsilon/\varepsilon\) is the fractional energy bin size. xstar has no information about the distance so assuming a standard distance of \(D=1\,\text{kpc}\) and renormalizing to an incident luminosity of \(10^{38}\,\text{erg}\,\text{s}^{-1}\) gives
The meaning of this quantity and the normalization can be better understood if we consider how xspec works in more detail. In order to calculate the model counts in bin \(m\), xspec multiplies the \(F_n^\text{mod}\) vector with the response matrix \(A_{nm}\) and a normalization factor \(\kappa\):
Note that xspec will always add such a normalization factor to each additive model so that absolute values distance and intrinsic luminosity do not need to be known at the time of the model calculation in order to fit an observed spectrum. xspec will renormalize the model internally and the \(\kappa\) can be made a free fit parameter like any other. The normalization factor \(\kappa\) will also account or effects like partial covering geometries with covering fractions smaller than unity.
The two equations above illustrate how \(\kappa\) is connected to physical parameters such as distance and luminosity. So, for example, if your \(\texttt{atable}\) model fits to the data with a normalization of \(\kappa=1\), and you used a luminosity in creating the table, then this would imply that your data was consistent with a full shell illuminated by a luminosity of \(10^{38}\,\text{erg}\,\text{s}^{-1}\) at a distance of \(1 \,\text{kpc}\), or it also could be a shell illuminated by a luminosity of \(10^{44}\,\text{erg}\,\text{s}^{-1}\) at a distance of \(1000\,\text{kpc}\).
For another example, let’s say you know the distance to the source is \(1 \,\text{kpc}\) and the luminosity is \(10^{38}\,\text{erg}\,\text{s}^{-1}\), but the best fit has an emitter normalization of \(\kappa=0.1\). This would suggest that rather than a full sphere, the emitter only subtends \(10\%\) of the solid angle around the source.
Obviously, the column density of the emitter is important also. If you do not know the column density, then a shell of column density \(10^{19}\,\text{cm}^{-2}\) illuminated by a \(10^{38}\,\text{erg}\,\text{s}^{-1}\) source will probably have very similar emitted X-ray spectrum to a shell of column density \(10^{20}\,\text{cm}^{-2}\) illuminated by a \(10^{37}\,\text{erg}\,\text{s}^{-1}\) source. If there are absorption features in the spectrum, then they may constrain the column density.
Running XSTAR2XSPEC
After this lengthy introduction, we now turn to how to actually use XSTAR2XSPEC. The parameter handling for XSTAR2XSPEC is designed to be as flexible as possible, in principle, limited only by the physical resources (RAM, disk space & CPU time) available on your machine. You have the choice of varying any of the physical parameters used as input in an XSTAR model. XSTAR2XSPEC is called from the command line very much like xstar and prompts the user in a similar way to enter parameters.
xstar2xspec [options]
The available options are currently:
\(\texttt{-save}\): Save the spectral FITS files, modifying the file name to include the value of the loopcontrol variable for better identification. Note this can use GBs of disk space.
\(\texttt{-verbose}\): Generate more diagnostic messages (applies to the XSTAR2XSPEC script only).
\(\texttt{-restart}\): Continues the XSTAR2XSPEC run using the previous run. Note that it does NOT check the integrity of the table files from the terminated run. This is the user’s responsibility.
The interface for entering the prompted parameters is very similar to
the one of XSTAR and should be already familiar. However, when
inspecting the $PFILES directory, user won’t find a parameter
file associated with XSTAR2XSPEC. This is because XSTAR2XSPEC is a
Perl script that calls a FTOOL programs that each have their own
parameter files. The relevant tools are:
xstinitable: This routine processes most of the input parameters and creates the parameter grid. It then produces a list of all XSTAR calls needed to calculate the full model. It also creates the (empty) FITS files that will store the model with appropriate dimensions. Changes to hidden xstar parameters need to be made through the parameter file “xstinitable.par”.
xstar2table: This routine takes care of the post-processing of each XSTAR run. It converts the information in the output spectral file “xout_spect1.fits” into emission spectra, absorption coeficients, optical depths etc., and fills them in into the table model. xstar2table is called every time an individual XSTAR calculation has finished.
Input Parameters
This section gives an overview over the input parameters to XSTAR2XSPEC and how they can be varied. We distinguish between parameters that are directly connected to physical properties of the photoionzed plasma and control parameters that affect the computational behavior of the code. Obviously only the physical parameters can be made variable to produce a spectral model that can be fitted to data.
Physical Parameters
XSTAR models are based on 21 physical parameters described in detail in Chapter User Input to XSTAR. Each of these parameters can be varied in different ways or held constant for all XSTAR runs. Due to the limitations of table models mentioned in earlier sections, user will usually only vary two to four interpolated parameters and possibily a few additional parameters, while all other parameters are held constant. The levels of variability for each parameter are:
Constant (variation type = 0): This parameter is held constant in all the XSTAR runs.
Additive (variation type = 1): The additive class of parameters provides a simple method for varying parameters that are reasonably independent of the others. For more info on how additive parameters function, see the XSPEC manual and see OGIP memo 92-009.
Interpolated (variation type = 2): Interpolated parameters provide the greatest accuracy in building table models. They also require the most processing time. For each interpolated parameter, you can define some number of points between a maximum and minimum range. The placement of these intermediate points is determined by the interpolation type – linear (interpolation type = 0) or logarithmic (interpolation type = 1).
It is useful to take a look at the parameter file xstinitable.par that controls the variation of each parameter. This file is many entries so we only show a few exerpts below for illustration. While for single XSTAR runs, it is not recommended to manipulate the parameter file, for XSTAR2XSPEC it can be a convient method to configure the table model calculation.
[...]
density,r,a,1.E+12,1.e+4,1.E21,"density soft maximum (cm**-3)"
densitytyp,i,h,0,0,2,"density variation type"
densityint,i,a,1,0,1,"density interpolation type"
densitysof,r,a,0.,1.e+4,1.e+18,"density soft minimum"
densitynst,i,a,1,1,20,"density number of steps"
[...]
Here, density is held constant (densitytyp=0). All XSTAR runs will use the (“soft maximum”) density of \(10^{12}\,\text{cm}^{-3}\).
[...]
column,r,a,1.e+24,1.e+10,1.E25,"column density soft maximum (cm**-2)"
columntyp,i,h,2,0,2,"column density variation type"
columnint,i,a,1,0,1,"column density interpolation type"
columnsof,r,a,1.E+21,1.,1.E25,"column density soft minimum"
columnnst,i,a,7,1,20,"column density number of steps"
rlogxi,r,a,5.,-10.,+10.,"log(ionization parameter) soft maximum (erg cm/s)"
rlogxityp,i,h,2,0,2,"log(ionization parameter) variation type"
rlogxiint,i,a,0,0,1,"log(ionization parameter) interpolation type"
rlogxisof,r,a,1.,-10.0,+10.0,"log(ionization parameter) soft minimum"
rlogxinst,i,a,9,1,40,"log(ionization parameter) number of steps"
[...]
Here, column density and ionization parameter are interpolated
parameters (variation type=2). The column density will be varied from
\(10^{21}\,\text{cm}^{-2}\) to \(10^{24}\,\text{cm}^{-2}\) in
seven steps. The ionization parameter will be varied from 1 to 5 in 9
steps. The column density will be varied in logarithmic steps
(interpolation type=1), since it spans several orders of
magnitude. The parameter rlogxi is already the log of the
ionization parameter and is therefore sampled in linear steps
(interpolation type=0).
[...]
feabund,r,h,1.,0.,100.,"iron abundance soft maximum"
feabundtyp,i,h,1,0,2,"iron abundance variation type"
feabundint,i,h,1,0,1,"iron abundance interpolation type"
feabundsof,r,h,0.,0.,1.,"iron abundance soft minimum"
feabundnst,i,h,1,1,20,"iron abundance number of steps"
[...]
In this example, the iron abundance is an additive parameter (variation type=1). The limits here are 0 and 1 with only one step. This means is all interpolated parameter combinations will be calculated once with the iron abundance set to 0 and once with the iron abundance set to 1.
Control Parameters
In addition to the parameters that define physical properties of the photoionized plasma, there is a number of non-physical parameters that control the behavior of the code. These parameters are never varied but are also included in the parameter file xstinitable.par so the user has full control over them. Additional parameters that are relevant to XSTAR2XSPEC but not to XSTAR are:
elow: This parameter determines the low energy end (in eV) of the spectrum selected from the XSTAR output files.ehigh: This parameter determines the high energy end (in eV) of the spectrum selected from the XSTAR output files.
[...]
elow,r,h,1.0E+2,0.,5.11E+5,"energy band low end (eV)"
ehigh,r,h,2.0E+4,0.,5.11E+5,"energy band high end (eV)"
[...]
As an example, spectral range here is \(100\,\text{eV--}20\,\text{keV}\).
Runtime Estimates
In estimating the running time of XSTAR2XPEC, the key factor is the number of times XSTAR is called. Consider a run with \(\text{N}_\text{I}\) interpolated parameters, where interpolated parameter \(i\) is evaluated at \(n_{i}\) points \((1 \leq i \leq {\rm N}_{\rm I})\). The total number of times XSTAR is called is then \(\prod\limits_{\text{i}=1}^{\text{N}_\text{I}} n_i\). However, if \({\rm N}_{\rm A}\) additive parameters are also defined, then for each set of interpolated variables, there is one run with all the additive parameters are zero and the remaining \(\text{N}_\text{A}\) runs have one of the additive parameters at it’s maximum value and the rest all zero. This means that the total number of times XSTAR must be run is given by
and this can provide you with a feel for how long a complete XSTAR2XSPEC run will require. As an example, if we defined 2 interpolated parameters (one evaluated at 5 points and the other evaluated at 4 points) and 8 additive parameters, XSTAR would be run a total of
which means that if the XSTAR runs for the appropriate model range averages five minutes, it will take approximately \(180\times5\) minutes = 900 minutes = 15 hours for a complete XSTAR2XSPEC run.
Output Files
XSTAR2XSPEC produces five output files:
xstar2xspec.log is a concatenation of all the xout_step.log files from all the xstar runs called by XSTAR2XSPEC.
xout_ain.fits is a FITS file containing the \(\texttt{atable}\) with the reflected emission spectrum produced by XSTAR.
xout_aout.fits is FITS file containing the \(\texttt{atable}\) with the emission spectrum in the forward (transmitted) direction produced by XSTAR.
xout_mtable.fits is a FITS file containing the \(\texttt{mtable}\) with the absorption spectrum in the forward (transmitted) direction produced by XSTAR.
xout_etable.fits is a FITS file containing the \(\texttt{etable}\) with the absorption spectrum in the forward (transmitted) direction produced by XSTAR.
All the FITS files are formatted that they can directly loaded in xspec. Here we load the forward emission \(\texttt{atable}\) model of grid18:
unix~> xspec
XSPEC version: 12.13.1
Build Date/Time: Fri Feb 9 12:12:40 2024
XSPEC12>model atable{xout_aout.fits}
[...]
========================================================================
Model atable{xout_aout.fits}<1> Source No.: 1 Active/Off
Model Model Component Parameter Unit Value
par comp
1 1 template column 1.00000E+21 +/- 0.0
2 1 template rlogxi 0.0 +/- 0.0
3 1 template cabund 1.00000 frozen
4 1 template nabund 1.00000 frozen
5 1 template oabund 1.00000 frozen
6 1 template neabund 1.00000 frozen
7 1 template mgabund 1.00000 frozen
8 1 template siabund 1.00000 frozen
9 1 template sabund 1.00000 frozen
10 1 template arabund 1.00000 frozen
11 1 template caabund 1.00000 frozen
12 1 template feabund 1.00000 frozen
13 1 template niabund 0.0 frozen
14 1 template z 0.0 frozen
15 1 template norm 1.00000 +/- 0.0
----------------------------------------------------------------------
***Warning: Magnitudes of parameters in model for data group 1
have been found to differ by more than 1e+10.
This can lead to significant roundoff errors during fit calculation.
XSPEC12>