# Simulate a High Energy Transmission Grating Spectrum of an X-ray Binary with ISIS

Proposer Threads (Cycle 16)

## Contents

**1. Thread Overview****2. Preliminary Considerations**-
**3.Simulate Source Spectrum** -
**4. Assessing the Simulations** **5. Assessing Pile-up****6. Complete Target Form****7. Further Resources****8. Thread Summary**-
**Images**

## 1. Thread Overview

This thread shows how to simulate a simple high energy transmission gratings (HETG) spectrum of an X-ray binary, and to assess whether it is significantly affected by pile-up. The simulations are based upon a successful Chandra proposal to study a low mass X-ray binary (LMXB). This source was detected with XMM-Newton, with the CCD spectra revealing flourescent Fe Kalpha emission at 6.4 keV, absorption from highly ionized Fe, as well as strong absorption from the interstellar medium. The high energies of these features are well-beyond the upper cutoff energies of the Reflection Gratings Spectrometer (RGS) on-board XMM-Newton, but are well-suited to the capabilities of HETG. The goal of the feasibility study is to determine the constraints that can be placed on the Fe features, as well as on the neutral column, using the high spectral resolution of the HETG instrument, when integrating over the X-ray bright phases of the binary orbit. Additionally, we wish to find out whether or not we must make any instrumental modifications to mitigate the possibility of photon pile-up.

During the bright phases of the binary orbit, which lasts approximately 65 ksec, the source was observed to have a relatively stable spectrum. We therefore only need to simulate a single spectrum. Similar procedures, however, could be followed for simulating the X-ray faint phases of the binary orbit. XMM-Newton observations show that the unabsorbed luminosity, presuming a distance of 15 kpc, is 4x10^37 ergs/s, yielding an unabsorbed flux of 1.5e-9 ergs/cm^2/s. This is sufficiently bright that the question of photon pile-up must be considered; however, observations also show NH=7.5x10^22 atoms/cm^2. This high column, coupled with the high energy features, implies that the High Energy Gratings (HEG - covering 0.7-8 keV) will be the prime Chandra instrument, although the Medium Energy Gratings (MEG - covering 0.4-7 keV) will be important for studies of the absorption column, and may even directly measure absorption edges (e.g., Si at 6.67 Angstrom, S at 5 Angstrom).

This thread uses an ISIS simulation. It you prefer to work with Sherpa please refer to the companion Sherpa thread.

## 2. Preliminary Considerations

Here are some things to consider when following this thread.- You need to have the Interactive Spectral Interpretation System
( ISIS) installed. This can be
installed from source code, or from binary. Similar simulation
procedures can be followed using Sherpa or XSPEC; however, ISIS
incorporates in a programmable manner access to atomic databases
that could be useful in subsequent analysis, it allows for
convenient in memory addition and regrouping of datasets, and should
photon pile-up be deemed an issue, it allows for modeling of
gratings pile-up. (See the Chandra ABC Guide to
Pileup and the CIAO Pileup AHELP Page for
further details.)

ISIS can be programmed using the S-lang interpretted language (http://www.s-lang.org) . This thread is saved as a script which a user can simply edit and customize, and rerun for their own source - Download the latest
instrument responses. The example and figure presented
below used Cycle 13
responses and ISIS 1.6.1-52 for illustration purposes. There
might be some minor numerical differences when using the latest
responses and some slight cosmetic differences in the ISIS prompts.

**Note**: the linked script references Cycle 15 response files - e.g., the script has a commandhma=load_arf("aciss_heg-1_cy15.garf")

as opposed to the examplehma=load_arf("aciss_heg-1_cy13.garf")

As a reminder, pre-made responses are provided for proposal purposes only and THEY SHOULD NOT BE USED FOR DATA ANALYSIS. - Check for previous observations of this source. You will need to justify additional observations if any are found. For detailed instructions, please refer to the "Proposing to Observe an On-Axis Point Source with ACIS" thread.
- Check for other bright sources in the field of view. The Observation Visualizer, ObsVis, available as part of CIAO, is used for this task. For further information, please refer to the the "Proposing to Observe an On-Axis Point Source with ACIS" thread. Note, however, that due to photon order sorting - i.e., the concept that for a dispersed gratings spectrum a specific location along the dispersion arms must correspond to a narrow range of energies in order to be associated with the spectrum from a known source location - gratings observations (especially of bright point sources) can be fairly tolerant of background and other sources within the field of view. Complicated regions with many sources, might require the use of MARX simulations (http://space.mit.edu/ASC/MARX/).
- An estimate of the neutral hydrogen column to the source can be obtained from colden if it is not already known. For this source, however, previous observations indicate an Nh of 7.5e22 cm^-2.

## 3.Simulate Source Spectrum

### Simulation Script

The steps in this thread have been incorporated into a script,
lmxb_sim.sl. *The easist way to follow
this thread is to cut and paste from the script, or edit the script
and run at the ISIS prompt:*

unix% isis Welcome to ISIS Version 1.6.1-52 Copyright (C) 1998-2011 Massachusetts Institute of Technology Isis web page: http://space.mit.edu/cxc/isis/ Mailing list archive: http://space.mit.edu/cxc/isis/archive/ Send questions to the mailing list: <isis-users@space.mit.edu>. For a summary of recent changes, type: "help changes" isis> require("xspec"); isis> .load lmxb_sim.sl

### Run ISIS and read in responses

The ISIS binary contains all the current XSPEC models, and if one builds ISIS from source, ISIS also can have access to any XSPEC local models (see the ISIS web page, http://space.mit.edu/ASC/ISIS, for details.) In this thread, we will only use standard XSPEC spectral models.

unix% isis Welcome to ISIS Version 1.6.1-52 Copyright (C) 1998-2011 Massachusetts Institute of Technology Isis web page: http://space.mit.edu/cxc/isis/ Mailing list archive: http://space.mit.edu/cxc/isis/archive/ Send questions to the mailing list: <isis-users@space.mit.edu>. For a summary of recent changes, type: "help changes" isis> require("xspec");

The CXC provides a set of instrument responses for use in proposal planning. For this proposal, we use a set of responses appropriate for the (plus and minus) first order HEG and MEG instruments. This consititutes eight files in all: two ancillary response files (ARF) each for HEG and MEG, and two redistribution matrix files (RMF) each for HEG and MEG. The ARF files are the most crucial for estimating count rates, but the RMF files are necessary to accurately estimate line widths.

Load the ARFs and RMFs, assign them to data set numbers 1-4, and then set the ARF exposure times to 650 ksec. (We also set the data exposure values to 650 ksec, although this is not strictly necessary - the fake data file integration time are set by the ARF exposures.) We choose a factor of 10 higher integration time than desired, at first, so that we have sufficiently good statistics to verify our settings of the model normalizations.

isis> hma = load_arf("aciss_heg-1_cy13.garf"); isis> hpa = load_arf("aciss_heg1_cy13.garf"); isis> mma = load_arf("aciss_meg-1_cy13.garf"); isis> mpa = load_arf("aciss_meg-1_cy13.garf"); isis> hmr = load_rmf("aciss_heg-1_cy13.grmf"); isis> hpr = load_rmf("aciss_heg1_cy13.grmf"); isis> mmr = load_rmf("aciss_meg-1_cy13.grmf"); isis> mpr = load_rmf("aciss_meg1_cy13.grmf"); isis> assign_arf(hma,1); % Dataset 1 = HEG-1 isis> assign_arf(hpa,2); % Dataset 2 = HEG+1 isis> assign_arf(mma,3); % Dataset 3 = MEG-1 isis> assign_arf(mpa,4); % Dataset 4 = MEG+1 isis> assign_rmf(hmr,1); % Dataset 1 = HEG-1 isis> assign_rmf(hpr,2); % Dataset 2 = HEG+1 isis> assign_rmf(mmr,3); % Dataset 3 = MEG-1 isis> assign_rmf(mpr,4); % Dataset 4 = MEG+1 isis> set_arf_exposure([hma,hpa,mma,mpa],6.5e5); isis> set_data_exposure(1,6.5e5); isis> set_data_exposure(2,6.5e5); isis> set_data_exposure(3,6.5e5); isis> set_data_exposure(4,6.5e5);

### Define the Model

ISIS for the most part uses model component names identical to those found in XSPEC. ISIS, however, does not allow the names to be abbreviated, and appends to each name "(#)", where # is an integer. By using different values of #, multiple model components (e.g., lines) of the same type can be introduced. If each model component is viewed as a vector of values (i.e., counts in each wavelength bin), then nearly any "mathematically sensible" combination of model components will be properly treated by ISIS as a complete model.

For this spectrum, we choose the sum of a blackbody plus a power law plus a gaussian (6.4 keV fluorescence) minus two gaussian absorption components (ionized Fe absorption), all multiplied by neutral column absorption. This is done with the ISIS "fit_fun" command:

isis> fit_fun("phabs(1)*(bbodyrad(1)+powerlaw(1)+gaussian(1)-gaussian(2)-gaussian(3))");

Invoking the ISIS "edit_par" command will open an editor window (e.g., using emacs or vi) wherein individual parameter values can be entered. Alternatively, one can enter parameter values on the command line with the ISIS "set_par" command. We take parameter values from the fits to the XMM-Newton data, and set the parameter minimum and maximum values to reasonable limits for subsequent fits (in order to prevent, for example, gaussians becoming very broad and fitting the continuum):

isis> set_par("phabs(1).nH", 0., 0, 0, 20); isis> set_par("bbodyrad(1).norm", 41., 0, 0, 100); isis> set_par("bbodyrad(1).kT", 1.18, 0, 0.5, 3); isis> set_par("powerlaw(1).PhoIndex", 2.02, 0, 1, 3); isis> set_par("gaussian(1).norm", 3.e-5, 0, 0, 0.01); isis> set_par("gaussian(1).LineE", 6.4, 0, 6, 6.7); isis> set_par("gaussian(1).Sigma", 0.1, 0, 0, 1); isis> set_par("gaussian(2).norm", 3.e-5, 0, 0, 0.01); isis> set_par("gaussian(2).LineE", 6.72, 0, 6.4, 6.9); isis> set_par("gaussian(2).Sigma", 0.02, 0, 0, 0.1); isis> set_par("gaussian(3).norm", 3.e-5, 0, 0, 0.01); isis> set_par("gaussian(3).LineE", 7.00, 0, 6.7, 7.3); isis> set_par("gaussian(3).Sigma", 0.02, 0, 0, 1);

Note that the paper describing the XMM-Newton spectra only give the
line equivalent widths, thus we initially set the line normalizations
to low values, and discuss below how we set the normalizations based
upon the published equivalent widths. Additionally, the XMM-Newton
paper quotes the unabsorbed flux. Thus, we begin by setting the
absorption to 0, to verify that we have properly understood the model
normalization described in the XMM-Newton paper. The XMM-Newton
spectra only give 90% confidence value *upper limits* to the
absorption line physical widths; we choose physical line widths 40% of
those values.

### Run fakeit, and Set Model Normalizations

Having defined a model (which, by default, is the same for all four gratings arms), we now create a fake dataset (complete with Poisson noise) with the ISIS "fakeit" command. The exposure time has already been set by assigning that exposure time to the respective ARFs. One then simply has to type "fakeit", followed by a "list_data" to verify that the datasets have indeed been created.

isis> fakeit; isis> list_data; Current Spectrum List: id instrument m prt src use/nbins A R totcts exp(ksec) target 1 -1 0 0 8192/ 8192 1 1 2.9792e+06 650.000 2 1 0 0 8192/ 8192 2 2 2.8209e+06 650.000 3 -1 0 0 8192/ 8192 3 3 5.5118e+06 650.000 4 1 0 0 8192/ 8192 4 4 5.5037e+06 650.000

The HEG detectors (datasets 1 and 2) have less effective area, and hence fewer counts, than the MEG detectors (datasets 3 and 4).

ISIS allows a great deal of customization. Here to help us set the normalizations, we define a function, "flux", that will return the flux (both as a print to the screen, and as a variable) in units of photons and ergs per cm2/sec. We also define a function, "eqw", that similarly returns line equivalent widths (for a given detector and model component normalization) in units of mAngstrom and eV.

isis> public define flux (id, e_lo, e_hi) { variable l_lo, l_hi, m, e_avg, pflux, eflux; l_hi = _A(e_lo); % Convert keV to A, or visa versa l_lo = _A(e_hi); % Convert keV to A, or visa versa m = get_model_flux(id); pflux = rebin (l_lo, l_hi, m.bin_lo, m.bin_hi, m.value)[0]; e_avg = _A(1.)*(1./m.bin_hi+1./m.bin_lo)/2.*1.602e-9; eflux = rebin (l_lo, l_hi, m.bin_lo, m.bin_hi, m.value*e_avg)[0]; () = printf("\n photons/cm2/sec: %11.3e, ergs/cm2/sec: %11.3e \n \n", pflux, eflux); return pflux, eflux; % Photon flux & ergs flux (per cm^2/sec) } isis> public define eqw(indx,par) { variable eva, eve; variable fv_hold = Fit_Verbose, pmin=0; Fit_Verbose = -1; () = eval_counts; variable a = get_model_flux(indx); variable phold = get_par_info(par); set_par(par,pmin,0,0,phold.max); () = eval_counts; variable b = get_model_flux(indx); set_par(par,phold.value,phold.freeze,phold.min,phold.max); () = eval_counts; variable iw = [0:length(a.bin_lo)-1]; variable diff = a.value - b.value; variable bdena = b.value/(a.bin_hi-a.bin_lo); variable bdene = reverse(reverse(b.value)/(_A(a.bin_lo)-_A(a.bin_hi))); variable fdena = diff/(a.bin_hi-a.bin_lo); variable fdene = reverse(reverse(diff)/(_A(a.bin_lo) - _A(a.bin_hi))); variable flux = sum( diff[iw] ); variable iwa = min(where( abs(fdena[iw]) == max(abs(fdena[iw])) )); variable iwe = min(where( abs(fdene[iw]) == max(abs(fdene[iw])) )); eva = flux/bdena[iw[iwa]]*1000.; eve = flux/bdene[iw[iwa]]*1000.; Fit_Verbose = fv_hold; () = printf("\n EqW (mA): %11.3e, EqW (eV): %11.3e \n \n", eva, eve); return eva, eve; } isis> (pflux,eflux) = flux(3,0.1,8); % Meg -1 flux, from 0.1 to 8 keV photons/cm2/sec: 6.530e-01, ergs/cm2/sec: 1.566e-09 isis> (g1_a,g1_e) = eqw(1,"gaussian(1).norm"); % Heg -1 Equivalent width EqW (mA): 8.094e-01, EqW (eV): 2.677e+00 isis> (g2_a,g2_e) = eqw(1,"gaussian(2).norm"); % Heg -1 Equivalent width EqW (mA): -8.533e-01, EqW (eV): -3.112e+00 isis> (g3_a,g3_e) = eqw(1,"gaussian(3).norm"); % Heg -1 Equivalent width EqW (mA): -8.960e-01, EqW (eV): -3.541e+00

We see that the unabsorbed flux is very close to the quoted value, therefore we do not alter the normalizations of the continuum components. The line equivalent widths, however, are all too small. The three features should have equivalent widths of 95 eV, -14 eV, and -23 eV (from low to high energy). We therefore scale the line normalizations by the ratio of desired to measured equivalent widths, run fakeit again, and then verify that the new equivalent widths are those that we desire.

isis> set_par("gaussian(1).norm",get_par("gaussian(1).norm")*95/g1_e); isis> set_par("gaussian(2).norm",get_par("gaussian(2).norm")*(-14)/g2_e); isis> set_par("gaussian(3).norm",get_par("gaussian(3).norm")*(-23)/g3_e); isis> fakeit; isis> (g1_a,g1_e) = eqw(1,"gaussian(1).norm"); % Heg -1 Equivalent width EqW (mA): 2.828e+01, EqW (eV): 9.500e+01 isis> (g2_a,g2_e) = eqw(1,"gaussian(2).norm"); % Heg -1 Equivalent width EqW (mA): -3.571e+00, EqW (eV): -1.302e+01 isis> (g3_a,g3_e) = eqw(1,"gaussian(3).norm"); % Heg -1 Equivalent width EqW (mA): -5.604e+00, EqW (eV): -2.215e+01

Although not exact matches to the desired equivalent widths, these values are well-within the XMM-Newton error bars, so we deem the normalizations sufficiently accurate for simulation purposes.

Finally, we set the neutral column to the value actually measured by XMM-Newton, we reduce the exposure time to 65 ksec, and we run fakeit once more. Furthermore, we save a parameter file, lmxb.par, with the final model.

isis> set_par("phabs(1).nH",7.5); isis> save_par("lmxb.par"); isis> set_arf_exposure([hma,hpa,mma,mpa],6.5e4); isis> set_data_exposure(1,6.5e4); isis> set_data_exposure(2,6.5e4); isis> set_data_exposure(3,6.5e4); isis> set_data_exposure(4,6.5e4); isis> fakeit; isis> list_data; Current Spectrum List: id instrument part/m src use/nbins A R totcts exp(ksec) target 1 -1 0 8192/ 8192 1 1 6.9668e+04 65.000 file: 2 +1 0 8192/ 8192 2 2 7.1021e+04 65.000 file: 3 -1 0 8192/ 8192 3 3 5.7348e+04 65.000 file: 4 +1 0 8192/ 8192 4 4 5.7378e+04 65.000 file:

## 4. Assessing the Simulations

### Fe lines viewed with HEG data

We begin by considering the Fe line features in the 6-7.5 keV region of the HEG spectrum. Since we are mostly concerned with just the line features, we perform a local fit for the continuum model in just the 1.65-2.3 Angstrom (5.4-7.5 keV) range. Given the restricted wavelength/energy range, we freeze the values of the neutral column, the blackbody temperature, and the power law photon index (i.e., we only allow continuum normalizations and line parameters to vary in the fit.) Furthermore, we at first will only consider fits to the higher spectral resolution HEG data.

To increase the signal-to-noise, we group the HEG data sets by a
factor of 2 and the MEG data by a factor of 4 (i.e., to about twice
and one times, respectively, the intrinsic resolution of the gratings
detectors). Additionally, for fitting purposes, we will combine the
two HEG datasets, and we will combine the two MEG datasets. * This
is an "in memory" combination only.* ISIS allows such regroupings of
data and combinations of data during analysis in order to facilitate
easier exploration of the optimal groupings.

isis> group_data([1,2]; min_chan=2); % Group HEG data by a factor of 2 isis> group_data([3,4], min_chan=2); % Group MEG data by a factor of 4 isis> heg = combine_datasets(1,2); % In memory combination of HEG data isis> meg = combine_datasets(3,4); % In memory combination of MEG data isis> notice_values([1,2],1.65,2.3; unit="a"); % Only notice the 1.65-2.3 A data for HEG isis> exclude(3,4); % Ignore MEG data for now isis> freeze([1,3,5]); % Freeze Nh, blackbody kT, Photon Index isis> fit_counts; Parameters[Variable] = 14[11] Data bins = 130 Chi-square = 118.9635 Reduced chi-square = 0.9996931 0 isis> list_par; phabs(1)*(bbodyrad(1)+powerlaw(1)+gaussian(1)-gaussian(2)-gaussian(3)) idx param tie-to freeze value min max 1 phabs(1).nH 0 1 7.5 0 20 10^22 2 bbodyrad(1).norm 0 0 40.88613 0 100 3 bbodyrad(1).kT 0 1 1.18 0.5 3 keV 4 powerlaw(1).norm 0 0 0.1447135 0 10 5 powerlaw(1).PhoIndex 0 1 2.02 1 3 6 gaussian(1).norm 0 0 0.0009578358 0 0.01 7 gaussian(1).LineE 0 0 6.394854 6 6.7 keV 8 gaussian(1).Sigma 0 0 0.1082421 0 1 keV 9 gaussian(2).norm 0 0 0.0001736271 0 0.01 10 gaussian(2).LineE 0 0 6.717939 6.4 6.9 keV 11 gaussian(2).Sigma 0 0 0.03438081 0 0.1 keV 12 gaussian(3).norm 0 0 0.0001914165 0 0.01 13 gaussian(3).LineE 0 0 7.000227 6.7 7.3 keV 14 gaussian(3).Sigma 0 0 0.02443937 0 1 keV isis> save_par("heg_fit.par"); % Save the HEG fit parameters isis> putenv("PGPLOT_BACKGROUND=white"); isis> putenv("PGPLOT_FOREGROUND=black"); % Make plots black on white isis> set_line_width(3); set_frame_line_width(3); charsize(1.1); % Nicer plot defaults isis> xrange(1.65,2.3); % Plot x-range to only show what we've noticed isis> ylog; % Logarithmic y-axis isis> rplot_counts(1); % HEG-1 data (although fit includes HEG+1); isis> rplot_counts(2); % HEG+1 data (although fit includes HEG-1); isis> hd = get_combined(heg, &get_data_counts); % Get the summed data isis> hm = get_combined(heg, &get_model_counts); % Get the summed model isis> label("Angstrom","Counts/Bin","Combined HEG Data"); isis> hplot(hd.bin_lo,hd.bin_hi,hd.value); isis> ohplot(hm.bin_lo,hm.bin_hi,hm.value);

The resulting plots are seen here (Figure 1), here (Figure 2), and here (Figure 3). Note that the ISIS commands for generating hardcopies of plots follows the format of:

isis> id = open_plot("file_name.ps/vcps"); % Create a color postscript file isis> resize(21,0.85); % 21 cm, with an aspect ratio of 0.85:1 ... plot commands ... isis> close_plot(id);

We now turn to determining the confidence bounds for the line parameters found above. For this, we use the ISIS "conf" command, which defaults to 90% confidence limits. 68% or 99% can also be specified. (Arbritrary values of the confidence level can be specified with the "fconf" command.) The conf command will return the upper and lower parameter values for these confidence intervals.

isis> conf("gaussian(1).norm"); 0.00114008 0.000802636 isis> conf("gaussian(1).LineE"); 6.41943 6.37659 isis> conf("gaussian(1).Sigma"); 0.135983 0.0904234 isis> conf("gaussian(2).norm"); 0.0003705 9.82817e-05 isis> conf("gaussian(2).LineE"); 6.73847 6.68033 isis> conf("gaussian(2).Sigma"); 0.0971562 0.0146643 isis> conf("gaussian(3).norm"); 0.000263035 0.000129759 isis> conf("gaussian(3).LineE"); 7.01259 6.98854 isis> conf("gaussian(3).Sigma"); 0.0396256 0.0138938

We find that the line positions are all determined with 90% confidence level uncertainties of 0.01-0.02 keV, and likewise the line widths are determined with confidences of 0.01-0.02 keV. In terms of velocities, this corresponds to approximately 500-1000 km/sec. Note that for one of the error bars searches, the fit failed, so we switched to a different (slower) fit method, and searched for the error bar again.

Above, we also find the uncertainties in the line normalizations. To
*properly* cast that in terms of equivalent width uncertainties, one
should individually change and freeze the line normalizations to their
confidence limits, re-fit the model, and then recalculate the
equivalent width; we show this below.

isis> set_fit_method("subplex"); % Slow fit method, but often useful for gratings data isis> (g1_a,g1_e) = eqw(1,"gaussian(1).norm"); EqW (mA): 2.616e+01, EqW (eV): 8.653e+01 isis> set_par("gaussian(1).norm",0.000803,1); () = fit_counts(); Parameters[Variable] = 14[10] Data bins = 130 Chi-square = 121.6589 Reduced chi-square = 1.013824 isis> (g1_a,g1_e) = eqw(1,"gaussian(1).norm"); EqW (mA): 2.176e+01, EqW (eV): 7.177e+01 isis> set_par("gaussian(1).norm",0.00114,1); isis> () = fit_counts(); Parameters[Variable] = 14[10] Data bins = 130 Chi-square = 121.6729 Reduced chi-square = 1.013941 isis> (g1_a,g1_e) = eqw(1,"gaussian(1).norm"); EqW (mA): 3.127e+01, EqW (eV): 1.037e+02 isis> load_par("heg_fit.par"); % Restore the HEG fit parameters isis> (g2_a,g2_e) = eqw(1,"gaussian(2).norm"); EqW (mA): -4.976e+00, EqW (eV): -1.810e+01 isis> set_par("gaussian(2).norm",0.00000983,1); isis> () = fit_counts(); Parameters[Variable] = 14[10] Data bins = 130 Chi-square = 137.6527 Reduced chi-square = 1.147106 isis> (g2_a,g2_e) = eqw(1,"gaussian(2).norm"); EqW (mA): -2.868e-01, EqW (eV): -1.043e+00 isis> set_par("gaussian(2).norm",0.00037,1); isis> () = fit_counts(); Parameters[Variable] = 14[10] Data bins = 130 Chi-square = 277.286 Reduced chi-square = 2.310717 isis> (g2_a,g2_e) = eqw(1,"gaussian(2).norm"); EqW (mA): -1.036e+01, EqW (eV): -3.767e+01 isis> (g3_a,g3_e) = eqw(1,"gaussian(3).norm"); EqW (mA): -5.799e+00, EqW (eV): -2.292e+01 isis> set_par("gaussian(3).norm",0.000130,1); isis> () = fit_counts(); Parameters[Variable] = 14[10] Data bins = 130 Chi-square = 121.6507 Reduced chi-square = 1.013756 isis> (g3_a,g3_e) = eqw(1,"gaussian(3).norm"); EqW (mA): -3.984e+00, EqW (eV): -1.574e+01 isis> set_par("gaussian(3).norm",0.000263,1); isis> () = fit_counts(); Parameters[Variable] = 14[10] Data bins = 130 Chi-square = 121.6708 Reduced chi-square = 1.013923 isis> (g3_a,g3_e) = eqw(1,"gaussian(3).norm"); EqW (mA): -7.861e+00, EqW (eV): -3.107e+01

In the above, we find the fitted line equivalent width to be: 87+17/-15 eV, -18-20/+8 eV, and -23-8/+7 eV, respectively. The original input values are encompassed within these error bars, but as expected from a Monte Carlo simulation, they are not exact matches.

### Spectrum Viewed with MEG data

We now consider the spectrum, specifically the fitted neutral column, by fitting the combined MEG+/-1 data. First, we exclude the HEG data from consideration, and only consider the 2.1-7.2 Angstrom (1.7-5.9 keV) range of the MEG data. This range has more than 20 counts per bin everywhere, and avoids the Fe line region (which is better viewed with the HEG data).

To facilitate the fitting process, we load the saved simulation parameters (which includes the lines), evaluate the model once (without fitting), and then rewrite the model without the line parameters. In this manner, we restore the unfrozen, starting fit values for the neutral column, blackbody, and power-law. (ISIS will always remember the last used parameter values for any model component, which makes it easy to subtract and add these components in a useful manner.)

isis> exclude([1,2]); % Ignore HEG isis> include(3,4); % Notice MEG isis> notice_values([3,4],2.1,7.2; unit="a"); % Only notice 1.7-7.2 Angstrom (1.7-5.9 keV) isis> load_par("lmxb.par"); % Simulation parameters, including lines isis> () = eval_counts; Parameters[Variable] = 14[14] Data bins = 509 Chi-square = 444.8406 Reduced chi-square = 0.898668 isis> fit_fun("phabs(1)*(bbodyrad(1)+powerlaw(1))"); % Don't fit lines isis> () = fit_counts; Parameters[Variable] = 5[5] Data bins = 509 Chi-square = 429.7703 Reduced chi-square = 0.8527188 isis> list_par; phabs(1)*(bbodyrad(1)+powerlaw(1)) idx param tie-to freeze value min max 1 phabs(1).nH 0 0 8.107716 0 20 10^22 2 bbodyrad(1).norm 0 0 34.06908 0 100 3 bbodyrad(1).kT 0 0 1.217122 0.5 3 keV 4 powerlaw(1).norm 0 0 0.3114823 0 10 5 powerlaw(1).PhoIndex 0 0 2.341737 1 3 isis> save_par("meg_fit.par"); % Save the MEG fit parameters isis> conf("phabs(1).nH"); 8.930898611813703 6.947107830366841 isis> title(""); % We set the plot title in "label" above, so remove it isis> xrange(2.1,7.2); % Only plot 2.1-7.2 Angstroms isis> rplot_counts(3); % MEG-1 data (although fit combined with MEG+1) isis> rplot_counts(4); % MEG+1 data (although fit combined with MEG-1) isis> flux_corr(3); % Calculated the "unfolded" spectrum for MEG-1 isis> plot_data_flux(3); % Plot the MEG-1 unfolded spectrum

Plots of the fits can be found here (Figure 4) and here
(Figure 5). Note that most of the structure seen in this
counts/bin spectrum is *instrumental*, as can be seen by comparing to
the "unfolded" MEG-1 spectrum shown here (Figure 6). Note,
however, that real structure is seen at the Si (6.67 A) and S (5 A)
edges.

With these fits, we find that the neutral column is found to be 7.0(+1.4)(-0.3)x10^22 cm^-2. There are fair number of counts shortward of the Si edge; therefore, a proposer could additionally explore what kinds of limits can be set on higher energy lines (e.g., Sulfur lines, using the MEG and HEG data).

## 5. Assessing Pile-up

Photon pile-up occurs when two or more photons land in the same, or neighboring, detector pixels in the same readout frame (see the Chandra ABC Guide to Pileup and the CIAO AHELP Pileup Page.) Pile-up is a concern for point sources viewed with the Chandra CCDs; however, it can also be a concern when using the gratings to view bright objects, such as X-ray binaries. Here we come up with a simple estimate for the peak pile-up fraction in these proposed observations.

Roughly speaking, the peak pile-up fraction (which translates to the
fraction by which the observed count rate *will be reduced* in the
observed spectrum) goes as the peak counts per 3 pixel region (along
the gratings arm) per frame readout time. This can be determined
graphically using ISIS.

The MEG detectors have the highest area, and the coarsest spectral resolution, of the HETG detectors, and hence typically have the greatest pile-up. Often this peak occurs close to the peak effective area of MEG, near 2 keV. However, given the very large neutral column to this particular source, the peak pile-up occurs at somewhat higher energies/shorter wavelengths (near 3 A/4 keV, as we show below).

The simulations indicate that the MEG+1 detector (dataset 4) has the
highest count rate of the two MEG detectors; therefore, we assess the
pile-up first for that detector as a "worst case" scenario. We use
the ISIS "plot_bin_density" plotting option to see the data plotted as
counts/sec/Angstrom, we plot with a linear y-scale, and we use the
ISIS "cursor" function, which draws interactive cross hairs on the
plot. The cursor functions lets us click on the peak of the plot, and
then prints the corresponding wavelength and y-value. (Note: the
cursor function should be used on a linear-linear scale plot in this
case.) We then multiply the y-value by a typical frame time (3.2
sec), 3 pixels (events are read in 3X3 pixel detector cells), and
0.011 (the number of Angstroms per MEG pixel), to arrive at an
estimate of the *peak* pile-up fraction.

isis> ylin; % Linear y-axis isis> plot_bin_density; % Create plots as counts/sec/Angstrom isis> plot_data_counts(4); % MEG+1 is the brightest arm isis> cursor; % Measure with the pgplot cursor Click for cursor coordinates, type 'q' in plot window to quit x= 2.917383e+00 y= 5.821421e-01 ch=D isis> 0.58*3.2*3*(0.011); % cts/s/A * 3.2 sec frame * 3 pixels * 0.011 A/pixel 0.061248

The plot used for this assessment can be seen here
(Figure 7). Thus, we find that for a *full readout frame
time of 3.2 sec*, the peak pile-up fraction is expected to be 6%, and
occurs near 2.9 Angstroms. Although not very severe, one can almost
always reduce the frame time for gratings observations down to 1.7 sec
by only using half of the chips along the chip readout direction.
This loses none of the dispersed spectra, and further reduces the peak
pile-up fraction to 3.4%. Given the very large neutral column, and
hence the very low flux at long wavelengths (i.e., far along the
dispersed arms, away from the center of the dispersed image), one
could consider using an even smaller sub-array and further reducing
the readout time and pile-up fraction.

## 6. Complete Target Form

For general instructions on how to submit a proposal, please see the Using RPS to Prepare and Submit a Chandra Proposal thread.

The above simulations provide details to use in the science justification (e.g., expected limits on Fe line wavelengths, physical widths, and strengths; estimates of the accuracy of the fitted neutral column), and can be used to fill out the proposal target forms.

For this observation, we will use the HETG detector with the ACIS-S
chips. To minimize pile-up, we will use only 1/2 of the array. The
1st order count rate can be found by summing the counts from the
simulations (use the "list_data" command, as shown above) and dividing
by 65 ksec. This yields 4.2 counts/sec. For the total field rate, we
note that approximately an additional 20% count rate will appear in
the higher orders, and for a source dominated by counts > 2 keV, the
0th order count rate will be approximately the same as the *total*
first order rate. (The 0th order image will suffer severely from
pile-up, hence reducing the detected rate. Calculation of that effect
is beyond the scope of this thread. By ignoring pile-up in the 0th
order image, we obtain a "worst case" estimate of the count rate. For
successful proposals, the Chandra gratings contact scientists will
work with the proposer to optimize settings to minimize the count rate
in the 0th order image.) Thus, as a very conservative estimate, we
use a total field rate of 10 counts/sec.

This count rate is sufficiently low such that we can use "Faint Mode" for the telemetry format. If the source were several times brighter, "Graded Mode" would be required. Furthermore, pile-up would increase. Sufficiently bright sources would require the use of Continuous Clocking mode. Currently, the calibration is better understood for Timed Event mode, as opposed to Continuous Clocking mode. It is the recommendation of the gratings team that so long as estimates place the peak pile-up fraction at less than 10%, Timed Event mode should be used instead of Continuous Clocking mode (unless sub-1.7 sec time resolution is required to meet the science goals of the observation). For brighter sources, one should consider using Continuous Clocking mode. Users can contact the gratings team for further advice.

Given the low effective area of the frontside chips at low energies/long wavelengths, the S0 chip, covering the low energies of HEG-1/MEG-1, can often be designated as an "optional chip". This is especially true in this particular case where the large neutral column already ensures very little flux at low energy.

The RPS form should have the following parameter values for this observation. If a parameter isn't listed here, use the default RPS value or leave the field blank. Note that the detector and SIM offsets have been chosen to place a custom 1/2 subarray on the half of the chips nearest the readout. This not only reduces the readout time, but minimizes Charge Transfer Inefficiency (CTI). The gratings contact scientists can help optimize these choices for successful proposals.

- SIM Translation Offset -- -6.14
- Total Observing Time (ksec) -- 65 ksec
- Count Rate -- 4.2
- Total Field Count Rate -- 10
- Exposure Mode -- TE
- Event Telemetry Format -- Faint
- CCDs On -- S0 = O1, S1-S5 = Y
- Subarray Type -- Use Custom Subarray
- Custom Subarray Start Row -- 1
- Custom Subarray End Ros -- 512

Use this link to view the completed Target Form.

## 7. Further Resources

All of the steps taken above have been placed in this S-lang script found here (lmxb_sim.sl). Save the script as lmxb_sim.sl in the same directory as the Cycle 13 response files, start ISIS in that directory, and then type:

isis> .load lmxb_sim.sl

The script can be edited to change the values of the parameters, and to save the plots as postscript files.

For further information about installing and running ISIS, see: http://space.mit.edu/ASC/ISIS.

For useful scripts, specifically utilities that will aid those more familiar with XSPEC, and useful plotting scripts (including ones that can handle both combined data and unfolded data), start at: http://space.mit.edu/ASC/ISIS/xspec_comparison.html.

## 8. Thread Summary

This thread shows how to simulate a simple gratings spectrum of an X-ray binary, wherein the spectrum shows narrow Fe emission and absorption lines at energies > 6 keV. Furthermore, it shows how to estimate the constraints that the MEG place upon measurements of the neutral column. Simulating a 65 ksec spectrum, we found that HEG constrains the line positions and widths to within 500-1000 km/sec. The neutral column is constrained to within 20%. We looked at the estimated level of pile-up in the dispersed spectrum, and concluded that with a 1/2 sub-array, the peak pile-up fraction will be less than 4%. Accordingly, we choose a custom 1/2 sub-array to perform the 65 ksec observation.