2.1.

Heasim Approach to Simulation

Heasim transforms a description of an astronomical source (or list of sources) into a flexible image transport system (FITS) file composed of a bintable with a list of “detected” X-ray events according to the characteristics of a particular instrument and details of the observation. The source is described by its location and distribution on the sky, the energy distribution and flux of its X-ray emission, and its variation in time. These properties are encapsulated in a source definition file (sdf) that may refer to auxiliary files provided by the user. Properties of the focal plane configuration [such as the detector field-of-view (FoV) and pixel size] are determined by the choice of X-ray mission and instrument and contained in a mission database file, whereas those having to do with the performance of the X-ray optics, detector response, and internal background are contained in supporting “calibration” files that are specified by heasim input parameter settings. The pointing direction in the sky and exposure time of the observation are selected using input parameter settings as well.

The steps that heasim takes in creating a set of simulated events for each input source, and for the background, are as follows. (1) The absorbed input spectrum (in $\text{photons}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}$ as a function of energy in keV) derived from one of several available spectral models, or from a user-defined spectrum file, is calculated in each input energy bin on a grid taken from the input ancillary response file (arf) that tabulates the instrument effective collecting area as a function of energy. (2) This is then converted to counts/bin by multiplying by the exposure time and by the effective area read from the arf. (3) For each source, the code loops over each input energy bin and over the counts in each bin, (4) probabilistically assigning a sky position to each event according to the source spatial distribution determined either by one of the available models or by a user-defined input image. (5) If a vignetting function, which quantifies the energy-dependent effective reduction of the X-ray telescope collecting area as a function of angle with respect to the telescope axis, is specified, events are discarded with a probability equal to one minus the vignetting function at the relevant energy and focal plane position. (6) Events are displaced from their incident locations according to a probability distribution determined by an input point-spread function (psf) or enclosed-energy function (eef), (7) discarded if located outside the FoV of the instrument, and (8) assigned final discrete sky coordinates that are resampled within detector pixels unless resampling is switched off. (9) Events in each energy bin not discarded due to vignetting or falling outside the FoV are redistributed into discrete output detector energy channels (referred to as PI or pulse-height invariant channels that account for detector gain and are related in a one-to-one manner with a range of photon energies) according to the input response matrix function (rmf) that defines the spectral redistribution from input energy to PI channel and the PI-to-energy mapping. These events are then (10) randomly assigned times according to the source temporal characteristics (constant, sinusoidal, or burst). The time assignment proceeds by mapping the light curve onto a uniform distribution using the integral of the light curve, employing a look-up table if necessary. For the sinusoidal option, this is conducted within a single cycle following selection of a random cycle within the exposure. For the burst, the look-up table spacing is adapted to resolve the steeply rising portion of the light curve. Sky background sources, i.e., from the output of skyback, are processed in the same manner. Internal background events from an input spectrum, following rescaling to the detector area and observation exposure time, are assigned random detector locations and times within the exposure. Events from all sources and from the internal background are merged and sorted on time, with an option to identify pileup based on a simple time interval criteria for events adjacent in time.

2.2.

Input

The heasim simulation input includes a number of required and optional parameters and supporting files that set the details of the simulated observation and simulation mode of operation, characterize the performance and properties of the instrument, and specify the source properties. We provide a brief description of these here; for more details, readers should consult the user guide.^10,11

The mission and instrument are specified in the heasim parameter file, as are the exposure time and pointing, and names of various optional and required files that contain the description of the source(s) (including any skyback output files describing background sources), and performance of the X-ray optics and detector (see below). Several operational options may set here as well, such as those needed for validation or observation splitting into subexposures (Sec. 2.3).

Heasim has two continuous, internal coordinate systems—one aligned with the plane of projection to celestial coordinates and centered on the pointing direction, and one rotated according to the roll angle to align with the instrument system. These are converted to discrete sky and focal plane coordinates according to the entries in the mission database file, heasim.mdb—an ASCII file that describes the detector orientation with respect to the telescope axis, and the criteria that determine how heasim resamples photons into pixels and discards those that end up outside the FoV or fall into gaps in the detector, for the supported instruments. The default heasim.mdb file is included in the refdata area of the HEASoft distribution.

Heasim users specify the spatial, spectral, and temporal characteristics of the sources to be simulated using the sdf, an ASCII file designed to represent a simple, compact, yet flexible, means of defining the spatial, temporal, and spectral properties of a source or ensemble of sources. Each source is described in the sdf by a single row that includes the source position, foreground hydrogen column density, flux (and bandpass over which that flux is defined), and describes the spectrum by one of five one-parameter models, or by a table provided by the user. Optionally, the user may also specify either one of the four possible extended-source functions and the corresponding set of parameters, or provide an image that heasim interprets as the true source distribution on the sky (otherwise the source is assumed to be point-like). Finally, sinusoidal or burst-like time variability may be specified. Currently, sources that are both time-variable and spatially extended are not supported. A fuller description of these source definition options may be found in Sec. 7.

Multiple-source simulations are realized simply by adding lines to the sdf. More complex sources with spatially varying spectra or multiple components with different temporal and/or spectral characteristics (such as a source with extended constant and point-like periodic components) are built up in a similar manner. For example, an extended source with a temperature gradient (in projection) is represented by an sdf with one row for each of a set of nonoverlapping, contiguous annuli.

Support for the input of user-defined image files and tabulated spectra enables the simulation of sources with arbitrary spatial and spectral characteristics; however, only a limited number of types of temporal variability are currently supported. The simx¹² (which shares some development heritage with heasim) and sixte¹³ simulation tools are designed to utilize simput files intended as a standard simulation data input format for X-ray telescopes;¹⁴ however, heasim currently is designed to be fully self-contained within HEASoft and does not support simput file input.

Supporting calibration files that describe the spectral response (effective area and energy redistribution) of the detector, the spatial response of the optics (redistribution of X-rays in the focal plane, vignetting function), and the internal background are specified by optional or required filenames in the heasim parameter file. These files are described in more detail in Sec. 8. A collection of standard prelaunch Hitomi calibration files, representing all of the standard formats, are included in a separate supporting files tarffile, ahsimfiles.tar.gz (Sec. 1).

2.3.

Operation and Output

As with any FTOOL, heasim may be run with parameters either prompted for, set on the command line, or preset using the pset command. Utilization of scripts streamlines the operation of repeated simulations for the same source using different detectors, and a set of example scripts are available for download (ahsimscripts.tar.gz; Sec. 1).

The simulation may be split into multiple subexposures by setting the heasim subexposure flag and duration parameters. In this case, continuity in phase or time offset from a burst is imposed, and the above steps are conducted for each subexposure, with the time-sorting, pileup calculation, and writing to output completed within each interval to reduce memory usage. If the subexposure duration parameter exceeds the total exposure time, the subexposure duration is “optimized” to include $<\mathrm{500,000}$ counts per subexposure based on the flux and bandpass in the sdf and the effective area. Time is assigned in blocks for each energy bin and events sorted in time, within each subexposure.

The output of heasim is a FITS event file¹⁵ with the following columns: TIME, the time of the event in seconds; $X$ and $Y$, the sky coordinates related to celestial coordinates by keywords written in the event file header; PI, the energy channel; and PILEUP, a flag set to 1 if the event is piled up (and 0 otherwise). Data products may be extracted from the simulated event file and analyzed using standard tools such as those in HEASoft. When XSELECT is used for extraction, the XSELECT environmental variable XSELECT_MDB should be set to xselect.mdb.heasim (included as part of the supporting files distribution). For extraction of spectra in particular, the coordinate used for the weighted map (“wmapname”) should be set to SKY coordinates $(x,y)$, and any region filtering must be applied using SKY or WCS coordinates.

3.1.

Input, Operation, and Output

No input is required other than setting the values in the parameter file. Skyback may be run with parameters either prompted for, set on the command line, or preset using the pset command. Utilization of scripts streamlines the operation of repeated simulations apropos to different detectors, and example scripts are available for download (ahsimscripts.tar.gz; Sec. 1).

Skyback output generally includes the following files: (1) a catalog of resolved point sources with source positions and spectra specifications for each source, in a format following that of heasim point sdfs; (2) a file that supplements the point source catalog with values of the redshift and intrinsic absorption; (3) a “catalog” with a single line representing the diffuse emission, in the heasim format for a flat extended source distribution; (4) the total diffuse (included unresolved point sources) background spectrum calculated on a user-specified energy grid in the form of a heasim user input spectrum ASCII file; and, (5) an XSPEC table model (FITS) file.

Because of the simplified treatment of the point source background, it is recommended to run skyback with different sets of parameters for different bandpasses. The selections in Table 1, where fabs0–5 represents the fraction of sources in each bin of intrinsic $N_H$, and fpar0–4 represents the fraction in each bin of source spectral index or opening angle (Sec. 9), roughly produce average values of the shape and normalization of the 0.5- to 2.0-keV and 2- to 10-keV astrophysical backgrounds, respectively. In combination with the Halo and LHB components, these are used to approximate the Hitomi SXI/SXS and HXI astrophysical backgrounds in the example in Sec. 5.2 below (though they are equally appropriate for any missions with similar bandpasses). The input spectra including all resolved, unresolved, and intrinsic diffuse components are shown in Fig. 1.

Table 1

Skyback example parameters.

Fig. 1

Composite input background spectra (all sources) used for the Hitomi SXI/SXS (black) and HXI (red) as described in Table 1.

5.1.

Bright Point Source

A Hitomi SXS simulation of a bright point source provides a simple illustration of a joint heasim/sxsbranch application. The source is assumed to have a (constant) Crab Nebula-like flux ($2\times {10}^{-8}\mathrm{erg}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}$ in the 2- to 10-keV bandpass) and a spectrum characterized by an absorbed ($N_H=3\times {10}^{21}{\mathrm{cm}}^{-2}$) power law (index 2.1). The source is assumed to be on-axis, i.e., the pointing is set to be the same as the source position, with roll angle (267.7 deg) and exposure time (9735 s) the same as the nominal roll angle and cleaned exposure time of the actual Hitomi SXS Crab Nebula observation. Two simulations are conducted: one using prelaunch estimates of the relevant rmf, arf, eef, and nxb; and one using postlaunch versions of these constructed as explained below. Both are compared to the results of standard SXS analysis of the actual Hitomi Crab Nebula observation.^23,24

The postlaunch rmf is calculated using the sxsmkrmf FTOOL^25,26 applied to the cleaned Crab Nebula SXS event file for HP and MP grades, and then renormalized. The postlaunch on-axis arf for the entire SXS detector array is calculated using the aharfgen FTOOL²⁶ and then rescaled using the psf fraction that is part of the xrtraytrace screen/log output (the rescaling is needed since heasim requires an effective area that corresponds to all incident events—whether or not they land on the detector and are ultimately included in the final simulated event file). The postlaunch energy-dependent on-axis eef is calculated using the eeftable FTOOL²⁷ by inputting the ray-tracing event file created in the course of constructing the arf. The postlaunch nxb is calculated using the sxsnxb FTOOL.²⁵

The prelaunch simulation is internally subdivided into 86 subexposures and the postlaunch simulation into 5 subexposures. The simulated event files are input into sxsbranch (Sec. 4), and the resulting event file along with the additional ITYPE (grade) column is used to extract spectra for (1) all pixels and (2) for HP and MP grades only (the recommended selection for bright sources). Various count rates and computed branching ratios are compared in Table 2 for the cleaned observed and simulated event files, as are the Poisson statistical estimates based on the count rates and expected distribution over pixels for a point source (Sec. 4), which are in excellent agreement with the simulated values. Observed and postlaunch simulated spectra are compared in Fig. 2. For the spectral comparison, the appropriate files are applied to the grade-filtered simulated spectrum, i.e., nxb and rmf scaled by the HP + MP branching ratio and the original (not rescaled) arf that relates count rates in the region defined by the SXS array into intrinsic source flux. The large reduction in total count rate from pre- to postlaunch estimates is due to the effect of the SXS gatevalve^26,28 in the latter. Some of the discrepancies between real and simulated postlaunch event files, and between estimated and calculated branching ratios (note also that the former are based on a prelaunch estimate of the eef and the corresponding assumed pixel distribution of counts), may be partly attributed to the idealized treatment for the simulations that assigns grades based solely on the relative timing of cleaned events. However, the excess measure in the observed spectrum with respect to that for the simulated fiducial Crab Nebula flux and high-energy steepening of the spectrum are not understood and under active investigation.²⁸

Table 2

Counts and branching ratios, with Poisson statistical estimates of the latter in parentheses.

Fig. 2

Observed (black, upper) and simulated postlaunch (red, lower) HP + MP Crab spectra and their ratios with respect to the model input into the simulator (histogram).

In the examples that follow, we apply prelaunch calibration files that are also more suitable for XRISM simulations, unless stated otherwise.

5.2.

Including Astrophysical Background

In this example, the output from skyback with parameters as described in Table 1 is input into heasim with the point source flux set to a negligible value and the exposure time to 500 ks in order to emphasize the background component. The Hitomi SXI, HXI1, and HXI2 images and spectra extracted from the entire FoV are shown in Figs. 3 and 4, respectively. The spectra are compared (not fit) with a composite XSPEC model that includes one component for each contribution as follows: absorbed and unabsorbed thermal (Raymond–Smith) plasma models representing the Halo and LHB, and the cutoff power law model fit to the 3- to 60-keV spectrum from the High-Energy Astronomical Observatory 1 derived by Gruber et al.²⁹ representing the extragalactic background. The thermal component parameters are based on a joint SXI/SXS fit with an added power law representing the extragalactic background; the Gruber et al. model²⁹ (originally normalized per steradian) is scaled to the particular FoV of each instrument.

Fig. 3

SXI, HXI1, and HXI2 images of the simulated X-ray background derived using skyback with the parameters displayed in Table 1.

Fig. 4

SXI and HXI2 (identical to HXI1) spectra of the simulated X-ray background derived using skyback with the parameters displayed in Table 1 compared to a model that includes LHB, Halo, and the extragalactic spectrum from Ref. 29.

5.3.

Constraining Nuclear Line Emission Extent

Figure 5 shows the FeK-region of the spectrum extracted from a simulation of a 200-ks SXS observation of a cluster of galaxies based on the Perseus cluster. The X-ray emitting components in the core include an active galactic nucleus (AGN) with an absorbed power law component and narrow 6.4-keV $\mathrm{FeK}\alpha$ emission feature, and an intracluster medium component (ICM) with spatial distribution modeled as a $\beta$-model with $\beta =0.53$ and core radius 1.26 arcmin.³⁰ The input ICM spectral model is constructed by fitting the actual Hitomi Perseus cluster SXS data with a model that includes a single temperature, variable abundance plasma with turbulent broadening,³¹ with the text output of the best-fit model specified in the sdf along with the $\beta$-model parameters. The sdf also includes two additional lines with power law and monoenergetic spectra representing the AGN continuum and redshifted Fe resonance line with fluxes based on the fit.

Fig. 5

6.2- to 7-keV spectrum extracted from prelaunch simulated Perseus cluster Hitomi SXS event file (black) along with total model (green) and contributions from the ICM (orange), AGN continuum (blue), and AGN $\mathrm{FeK}\alpha$ (red) components, highlighting the ICM Fe He-like $\mathrm{K}\alpha$ and H-like $\mathrm{Ly}\alpha$, and the AGN neutral Fe, features.

The spectra from the previous example could have also been produced from simple spectral-only simulations. An illustrative scientifically related application of heasim (in conjunction with sxsbranch) would be to assess the feasibility of constraining the extent of the nuclear $\mathrm{FeK}\alpha$ emission by comparing simulations with various exposure times, and assuming various spatial distributions for this component. For example, one may apply sxsbranch to heasim output events from a series of Perseus cluster-like simulations as described in Sec. 5.3, with the same $\mathrm{FeK}\alpha$ line flux but with different spatial distributions, selecting the high-energy resolution events in a narrow band, and then comparing the distributions of counts over pixels in the array. Figure 6 shows the results from two such simulations: one with a point-like $\mathrm{FeK}\alpha$ distribution, the other an identical simulation where the $\mathrm{FeK}\alpha$ emission is assumed to be extended in a constant surface brightness disk 1 arcmin in radius ($\sim 20\mathrm{kpc}$ at the redshift of the Perseus cluster). HP, MP, and MS grade events are selected in a 7-eV wideband centered on the redshifted centroid of the emission line. Figures 7 and 8 show these same distributions as heat maps. An $F$-test comparison of the two distributions yields $p$-values of $\sim 0.068$ for the full array and $\sim 0.021$ for the innermost 16 pixels; a Monte Carlo approach may be used for a more complete assessment as each simulation runs to completion in $<1\min$ on a MacBook Pro with 2.5 GHz quad-core Intel Core i7 processor and 16 GB 1600 MHz DDR3 of memory. A full spectral decomposition in each pixel to better isolate the line component may also be conducted, as was done in the analysis of the actual Hitomi analysis.³²

Fig. 6

Comparison of pixel (0–11, 13–35) fractional distributions of HP + MP + MS counts in a 7-eV wideband, centered on the redshifted centroid of the 6.4-keV $\mathrm{FeK}\alpha$ emission feature, for simulated Perseus cluster SXS event files (Sec. 5.3) This $\mathrm{FeK}\alpha$ component is including in addition to ICM and AGN continuum components and assumed to be either point-like (blue histogram, left) or extend to 1 arcmin in radius (green histogram, right).

Fig. 7

Heat map for the pixel distribution shown as the blue histogram in Fig. 6, with the narrow band HP + MP + MS counts indicated in each pixel.

Fig. 8

Heat map for the pixel distribution shown as the green histogramin Fig. 6, with the narrow band HP + MP + MS counts indicated in each pixel.