Simulations with a simplified response matrix

I present a simplified response for Arcus simulations, consisting of just four files: ARF/RMF for zero order and ARF/RMF summed over all dispersed orders. For simplicity, I sum the effective area of the relevant dispersed spectra and average the resolving power. This ignores problems with order-sorting, which will be relevant for detecting weak absorption lines. Also, in reality, different orders will have different spectral resolution. If order-sorting (continuum sources!) or very high spectral resolution (> 3000) is important, simulations should be done for individual orders, not with this merged response.

from nbtemplate import display_header, get_path
display_header('SimSimple.ipynb', status='Arcus Probe')

Revision status: Arcus Probe

Last revision in version control:

• Author: H. Moritz Günther hgunther@mit.edu
• commit 807a8cb61a0b79d0cb60b4ab3204e24c0afa0837
• Date: Sat Mar 4 07:50:31 2023 -0500

This document is git version controlled. The repository is available at https://github.com/hamogu/arcus. See git commit log for full revision history.

Code was last run with:

• MARXS ray-trace code version 1.3.dev109 (commit hash: 675e82d) (Date of dirty version: 20230301)
• ARCUS CALDB version 3bec96a (2023-02-15 18:14:20)

If you are reading a static version of this document (e.g. as a readme distributed with arf/rmf files), check https://space.mit.edu/home/guenther/ARCUS/SimSimple.html for updates.

The need for a simple simulations approach

Arcus has four different channels, and in each channel we detect > 15 dispsered orders. Since Arcus has two cameras, some of the signal is detected quite close to the zero order and has resolving powers $R$ of just a few hundred, while most of the signal is dispersed to much larger distances, where we expect to reach resolving powers > 3500. For each order, $R$ changes with wavelength, and at any given wavelength, we will detect signal in different orders simultaneously at different values of $R$.

To make matters more complicated, order-sorting is not perfect and significant order confusion will appear, i.e. in order -5, we will also see some signal from order -4 and -6, which is dispersed in a different way. Thus, the full ARF/RMF set to simulate all this needs > 180 ARF/RMF pairs (4 channels x 15 orders x 3, where the last factor acounts for the extracted order and the interloper orders from insufficient order sorting). This clearly is too complicated for most simulations. Thus, I provide here a simplified set of ARF/RMF that averages over all grating orders.

Caveats

By averaging, we loose some information and brush over some important science implications. However, I expect that the ARF and RMF files presented here will be sufficient to study the majority of science cases.

Reduced peak resolving power, but overestimate of effective area at intermediate resolving power

I add up effective areas and average $R$ over all relevant orders (using the effective area in each order to weight the average). Before averaging, some orders provide a better $R$, sometimes significantly so. For example, at 20 Å, Arcus will see several orders with different $R$. For a science investigation that needs to separate two close emission lines, the average resolving power might be insufficient, but it maybe could still be done by using only the order that provides the best $R$ (at the cost of a lower effective area). Such a science case needs to be simulated using "per order" arf and rmf files. On the other hand, not all orders have the average $R$. A science investigation that requires exactly $R = R_\mathrm{average}$ will find that in practice some orders have higher and some orders have lower $R$. If the orders with lower $R$ are insufficient and only some orders can be used, then the exposure time needs to be longer to account for that.

Order confusion

Dispersed Arcus orders will land at the same physical location of the detector and orders need to be sorted using the CCD energy resolution. Since orders are closely spaced in energy space, this order sorting is not easy. Either significant effective area is lost by extracting only very narrow regions in energy space (up to factors of a few in the most extreme scenarios!) or a significant fraction (20% for some scenarios with broad order-sorting regions) of the photons are sorted into the wrong order. This is not as much a problem for emission line sources since comparison of different orders can help to identify lines seen in the wrong order, but for sources with absorption lines on a strong continuum, this will effectively act as an extra background, and increases in Poisson noise in the continuum measurement.

What is included in the ARF and RMF?

ARF and RMF are calculated based on MARXS ray-tracing. Details of this ray-tracing are described at Moritz home page. In short, the procedure is as follows: A ray-trace simulation is run for a monoenergetic point source at the nominal aimpoint position for a large number of photons. These photons are then propagated through a simple mirror model and the gratings to the CCDs. The ratio of the number of input photons vs. the number of detected photons gives the ARF. The spatial distribution can be fitted and gives the input to the RMF. Simulations are run using a large number of photons on a relatively coarse grid (1 Å). That grid is sufficient to determine the factors that depend mostly on the geometry (e.g. the shadows cast by the grating support structure). The spot sizes are fitted with a combination of Gaussian and Lorentzian components and stored in LSFPARM files in exactly the same format that the Chandra CALDB uses for this purpose. When generating the RMF on an arbitary energy grid, these factors are interpolated. For the ARF, the "geometric" factors from the ray-trace are interpolated on the new energy grid, the position of the CCD edges in wavelength space is calculated analytically, and then multiplied by the factors that depend only on the photon energy (such as the QE of the CCDs). The latter factors are given at a higher energy resolution. Thus, the resulting ARF will properly describe sharp features such as edges in the CCD response, but not have the same random noise that would happen if the ARF and RMF files were constructed from Monte-Carlo simulations are that run independently for every wavelength bin. (Also, this approach requires several orders of magnitude less computing time.)

The follwing effects are included in ARF and RMF:

• Sum of effective area for order 2, 1, and -1 to -11 (0 order is treated separately)
• Arcus Jan 2023 SPO size and location
• SPO coating
• SPO in-plane and off-plane Gaussian scatter
• 5.7 $\mu\mathrm{m}$ deep CAT gratings with nominal 1.8 deg blaze (CAT locations from Ed Hertz)
• 4 channels, offset such that chip gaps never overlap (channel placement from optics telecon in Jan 2023),
• newest filter and CCD QE curves by Eric Miller
• contamination estimate by Meg
• all geometric shadowing (from CAT grating mount under SPO)
• CAT grating L1 and L2 blockage
• L1 cross-dispersion and L2 dispersion
• nominal residual pointing jitter (=uncertainty in pointing knowledge), but no jitter beyond that (i.e. constant aspect solution, so chip-gaps are not filled in by jitter)
• mechanical misalignment according to Arcus alignment budget

What it not included?

• order sorting and losses due to OSIP filtering
• choice of spatial extraction regions. Currently, I'm using all photons that are not cross-dispersed, but that is at the same time too restrictive (for low background, one can get more signal by including a wider region), and not ambitious enough ($R$ could be enhanced by extracting a narrower region).
• chip edges are included, but treatment of the exact edge is not fully correct (should not be this sharp due to finite RMF)
• cross-channel contamination
• off-axis / extended sources
• CCD details like read-out streaks, hot pixels etc.

How do ARF and RMF look?

import numpy as np
import astropy.units as u

import matplotlib.pyplot as plt
import matplotlib as mpl

from sherpa.astro import ui

%matplotlib inline

# Improve labelling of log scales
mpl.rcParams['axes.formatter.min_exponent'] = 2

arfrmfpath = '../../arfrmf_20230228/'

date = '2021-03-30'  # Date in filename
arf0 = ui.unpack_arf(arfrmfpath + 'chan_all_+0.arf')
arfd = ui.unpack_arf(arfrmfpath + 'merged/chan_all_dispersed.arf')
rmf0 = ui.unpack_rmf(arfrmfpath + 'chan_all_+0.rmf')
rmfd = ui.unpack_rmf(arfrmfpath + 'merged/chan_all_dispersed.rmf')

fig, ax = plt.subplots(ncols=2, figsize=(10, 4))
ax[0].semilogx(arf0.x, arf0.y, label='order 0')
ax[0].semilogx(arfd.x, arfd.y, label='dispersed')
ax[0].set_xlabel(arf0.get_xlabel())
ax[0].set_ylabel('$A_\mathrm{eff}$ (' + arf0.get_ylabel() + ')')
ax[0].get_xaxis().set_minor_formatter(mpl.ticker.LogFormatterMathtext(labelOnlyBase=False,
                                                                   minor_thresholds=(2, .5)))
ax[0].tick_params(axis='x', labelsize=mpl.rcParams['xtick.labelsize'], which='both')
out = ax[0].legend()
ax[0].grid(True)

ax[1].plot((arf0.x * u.keV).to(u.Angstrom, equivalencies=u.spectral()), arf0.y, label='order 0')
ax[1].plot((arfd.x * u.keV).to(u.Angstrom, equivalencies=u.spectral()), arfd.y, label='dispersed')
ax[1].set_xlabel('wavelength ($\AA$)')
out = ax[1].set_ylabel('$A_\mathrm{eff}$ (' + arf0.get_ylabel() + ')')
out = ax[1].grid(True)

This figure shows the effective area for the dispersed signal and the direct light (0 order) over energy (left) and wavelength (right). The 0 order has significant throughput at high energies. At low energies, the dispersed order dominates. Note that the wiggles as low energies for the dispersed order are not noise, but due to some order in some channel falling in a chip gap. However, channels are arranged such that the chip gaps occur at different wavelength for each channel. That way, none of the chip gaps is very deep in the summed effective area curve.

The next plots show five representative line spread functions for the 0 order and the dispersed spectrum and list some properties of the RMF like the energy range covered and the number of bins. Note that these plots use a logaritmic scale in energy, which makes the RMF look narrower at higher energies, although the FWHM actually increses with energy.

rmf0

<DataRMF data set instance '../../arfrmf_20230228/chan_all_+0.rmf'>

rmfd

<DataRMF data set instance '../../arfrmf_20230228/merged/chan_all_dispersed.rmf'>

The figure shows five exemplary line spread functions for the averaged RMF of the dispersed orders. Note that the y axis is logarithmic. The line spread function is narrow for low energies (large wavelength) that are seen in the higher orders and wider for higher energies, that are only see in order 1 or 2.

Check the spectral resolution

As a simple illustration of the spectral resolution, I fake a single, very narrow emission line. Then, I can check $R=\frac{\lambda}{\Delta\lambda}$ where we take the FWHM of the line as $\Delta\lambda$. For this and all following simulations, I'm using Sherpa, which is distrubuted with CIAO. Use the "toggle on/of code" button on the top of this page (interactive version) to see the commands. Of course, this is just an example. The ARF and RMF files are OGIP compliant and any other X-ray modelling package (XSPEC, ISIS, SPEX, ...) should work similarly.

Dispersed spectra

ui.set_source('fake_gauss', ui.gauss1d.g1)
g1.pos = 20.0
g1.fwhm = 0.0001
g1.ampl = 100000

# Set analysis only works if a dataset has been loaded already, so we make one first
# and then overwrite it again after calling "set_analysis"
ui.fake_pha('fake_gauss', arfd, rmfd, 5)
ui.set_analysis("wave")
ui.fake_pha('fake_gauss', arfd, rmfd, 1e4)

ui.ignore(None, None)
ui.notice(g1.pos.val-0.01, g1.pos.val+0.01)
#ui.ignore(None, g1.pos.val-0.01)
#ui.ignore(g1.pos.val-0.01, None)
ui.plot_fit('fake_gauss')
ax = plt.gca()
#out = ax.set_xlim(g1.pos.val-0.01, g1.pos.val+0.01)

dataset BPTau: 2.5:59.4995 Wavelength (Angstrom) -> no data

INFO:sherpa.ui.utils:dataset BPTau: 2.5:59.4995 Wavelength (Angstrom) -> no data

dataset fake_gauss: 2.5:59.4995 Wavelength (Angstrom) -> no data

INFO:sherpa.ui.utils:dataset fake_gauss: 2.5:59.4995 Wavelength (Angstrom) -> no data

dataset fake_gauss0: 1.24046:61.9921 Wavelength (Angstrom) -> no data

INFO:sherpa.ui.utils:dataset fake_gauss0: 1.24046:61.9921 Wavelength (Angstrom) -> no data

dataset BPTau: no data -> 19.9895:20.0105 Wavelength (Angstrom)

INFO:sherpa.ui.utils:dataset BPTau: no data -> 19.9895:20.0105 Wavelength (Angstrom)

dataset fake_gauss: no data -> 19.99:20.01 Wavelength (Angstrom)

INFO:sherpa.ui.utils:dataset fake_gauss: no data -> 19.99:20.01 Wavelength (Angstrom)

dataset fake_gauss0: no data -> 19.8375:20.16 Wavelength (Angstrom)

INFO:sherpa.ui.utils:dataset fake_gauss0: no data -> 19.8375:20.16 Wavelength (Angstrom)

The plot shows a faked line at 20 Å. Orange is the input model (in the binning used in the ARF and RMF) and blue are simulated photons counts, where I have simulated so many photons, that the Poisson noise is negligible. I am now taking the x (wavelength) and y (count rate) numbers in the plot above and fit a Gaussian line to it without folding it through the RMF. That way, I can measure the observed line profile.

pl = ui.get_data_plot('fake_gauss')

from sherpa.models import NormGauss1D
from sherpa import stats
from sherpa import optmethods
from sherpa.data import Data1D
from sherpa.fit import Fit

d = Data1D('verify_R', pl.x, pl.y)
howwide = NormGauss1D('verify_R')
# Start close to final values to speed up convergence
howwide.fwhm = 0.01
howwide.pos = g1.pos.val
howwide.ampl = np.max(pl.y)
# Need to increase precision over default for this narrow line to make fit work well
mylevmar = optmethods.LevMar()
mylevmar.config['epsfcn'] = 1e-10
fit = Fit(d, howwide, method=mylevmar)

fit.fit()

<Fit results instance>

ui.plot_data('fake_gauss')
ax = plt.gca()
ax.plot(pl.x, pl.y, label='Input RMF model')
ax.plot(pl.x, howwide(pl.x), label='Gaussian fit')
ax.set_yscale('log')
ax.legend()
out = ax.set_ylim(500, 1e6)

The plot above shows the Gaussian fit to the simulated line (note that the y-axis is logarithmic now). The fit does a reasonable job at describing the line core and thus can be used to calculate $R$, but the line has extended wings. (The RMF is generated from a model that has both Gaussian and Lorentzian components because a Gaussian alone is insufficient to describe the wings seen in Arcus ray-traces.)

print (f'From the line shape, I calculate R= {howwide.pos.val / howwide.fwhm.val:5.0f}')

From the line shape, I calculate R=  2978

Zero-order spectra

I now repeat the same steps to check a line emission in the zero order.

ui.set_analysis('energy')
ui.set_source('fake_gauss0', ui.xsgaussian.g2)
g2.LineE = 1.
g2.Sigma = 0.001
g2.norm = 100000
ui.fake_pha('fake_gauss0', arf0, rmf0, 1e6)

pl = ui.get_data_plot('fake_gauss0')
d0 = Data1D('verify_R0', pl.x, pl.y)
howwide0 = NormGauss1D('verify_R0')
howwide0.fwhm = 0.1
howwide0.pos=1.1
howwide0.ampl = 1e4
fit = Fit(d0, howwide0)

fit.fit()

WARNING: fit failed: improper input parameters

WARNING:sherpa.optmethods:fit failed: improper input parameters

<Fit results instance>

out = plt.plot(pl.x, pl.y, lw=4)
out = plt.plot(pl.x, howwide0(pl.x))
out = plt.xlabel('Energy (kev)')
out = plt.ylabel('Counts/s/keV')
out = plt.xlim([0.9, 1.1])
out = plt.title('Fitting faked line with a Gaussian')

print (f'From the line shape, I calculate FWHM= {howwide0.fwhm.val * 1000:5.0f} eV.')

From the line shape, I calculate FWHM=   100 eV.

The input value for generating the RMF at this energy is 70 eV, so we can be confident that the ARF and RMF are working as intended.

A science simulation

As an example, I'm going to simulate a young star (BP Tau), applicable to the Arcus science case. For this star, we might want to use Arcus to check the density in the O VII triplet to see if the emission is coronal or comes from the accretion shock, where the densities are so high that the f/i ratio is reduced compared to the low-density coronal environment. This test has been done with XMM-Newton using a 130 ks exposure (Schmitt et. al, 2005), but with Arcus we could do it in a much shorter exposure time and possibly obtains a series of spectra to test for time variability. The input model for BP Tau is taken from Robrade & Schmitt (2006).

ui.set_xsabund('angr')  # Robrade and Schmitt use a mixuture of Anders & Grevesse and Grevesse & Sauval abundances
                        # in their paper, but for the purpose of this simulation, we ignore that detail.

 Solar Abundance Vector set to angr:  Anders E. & Grevesse N. Geochimica et Cosmochimica Acta 53, 197 (1989)

import numpy as np
import astropy.units as u

bptaumodel = ui.xsphabs.abs1 * (ui.xsvapec.v1 + ui.xsvapec.v2 + ui.xsvapec.v3)

abs1.nH = 0.15
v1.Fe = 0.28
v1.Si = 0.14
v1.O = 0.62
v1.Ne = 1.47
# couple all elements
for v in [v2, v3]:
    for p in v1.pars[1:-1]:  # leave out kT and norm
        setattr(v, p.name, getattr(v1, p.name))
        
v1.kT = 0.2
v2.kT = 0.63
v3.kT = 2.17

# XSPEC APEC normalization is 1e-14 / 4 pi d^2 with d in cm
normfac = 1e-14 / (4 * np.pi * ((130 * u.pc).to(u.cm).value)**2)

v1.norm = 3.48e52 * normfac
v2.norm = 5.28e52 * normfac
v3.norm = 10.3e52 * normfac

ui.notice(None, None)
ui.set_source('BPTau', bptaumodel)
ui.fake_pha('BPTau', arfd, rmfd, 25e3)
ui.set_analysis("wave")

dataset BPTau: 0.619596:0.620247 -> 0.208379:4.95937 Energy (keV)

INFO:sherpa.ui.utils:dataset BPTau: 0.619596:0.620247 -> 0.208379:4.95937 Energy (keV)

dataset fake_gauss: 0.619611:0.620231 -> 0.208379:4.95937 Energy (keV)

INFO:sherpa.ui.utils:dataset fake_gauss: 0.619611:0.620231 -> 0.208379:4.95937 Energy (keV)

dataset fake_gauss0: 0.615:0.625 -> 0.2:9.995 Energy (keV)

INFO:sherpa.ui.utils:dataset fake_gauss0: 0.615:0.625 -> 0.2:9.995 Energy (keV)

ui.group_width('BPTau', 3)
ui.plot_data("BPTau")
ax = plt.gca()
ax.set_xlim(21.55, 22.15)
out = ax.set_ylim(None, .5)

The plot shows an Arcus simulation for a 25 ks observation with no background. The signal is not great because the absorbing column denstiy really matters for the O VII lines, but the lines are clearly detected. The total count numbers are similar to the much longer XMM-Newton observation, but with Arcus, we obtain a much better wavelength resolution, allowing us to test for line shifts. In an accretion model, we might expcct the lines to the red-shifted, but the velocity resolution in XMM is insufficient to test that. Maybe we should re-run our simulation for a slightly longer exposure time to ensure that we won't be limited by counting statistics when the measure the line centroid?

However, that is beyond the scope of the current notebook.

Further checks

Below are a few extra plots that check consistency between different simulations. Since we use different systems to simulate Arcus that makes slightly different assumptions, disagreement between simulations is not necessarily a problem. For example, when extracting data for a different cross-dispersion width or different order sorting, then the total effective area will be somewhat different. So, these curves can give somewhat of an idea how slightly different choices for these parameters change the predicted outcome.

from marxs.missions.arcus.analyze_grid import orders_from_meta
from astropy.table import Table
import os
tarc = Table.read(os.path.join(get_path('arcus'), 'raygridRAeff.fits'))
tarc.orders = orders_from_meta(tarc.meta)


fig, ax = plt.subplots()
ax.plot((arf0.x * u.keV).to(u.Angstrom, equivalencies=u.spectral()), arf0.y, label='order 0')
ax.plot((arfd.x * u.keV).to(u.Angstrom, equivalencies=u.spectral()), arfd.y, label='dispersed ARF')

ax.fill_between(tarc['wave'], np.sum(tarc['Aeff4'][:, tarc.orders != 0], axis=1), 
                label='photon counting', color='0.7')

ax.set_xlabel('wavelength [$\AA$]')
out = ax.set_ylabel('$A_\mathrm{eff}$ (' + arf0.get_ylabel() + ')')
ax.get_xaxis().set_minor_formatter(mpl.ticker.LogFormatterMathtext(labelOnlyBase=False,
                                                                   minor_thresholds=(2, .5)))
ax.tick_params(axis='x', labelsize=mpl.rcParams['xtick.labelsize'], which='both')
out = ax.legend()