Thermal Emission

The first example we’ll look at is that of thermal emission from a galaxy cluster. In this case, the gas in the core of the cluster is “sloshing” in the center, producing spiral-shaped cold fronts. The dataset we want to use for this example is available for download from the yt Project at this link.

First, import our necessary modules:

[1]:
import yt
import pyxsim
import soxs

Next, we load the dataset with yt. Note that this dataset does not have species fields in it, so we’ll set default_species_fields="ionized" to assume full ionization (as appropriate for galaxy clusters):

[2]:
ds = yt.load(
    "GasSloshing/sloshing_nomag2_hdf5_plt_cnt_0150", default_species_fields="ionized"
)
yt : [INFO     ] 2024-03-06 15:35:05,235 Parameters: current_time              = 1.1835090993823291e+17
yt : [INFO     ] 2024-03-06 15:35:05,235 Parameters: domain_dimensions         = [16 16 16]
yt : [INFO     ] 2024-03-06 15:35:05,236 Parameters: domain_left_edge          = [-3.70272e+24 -3.70272e+24 -3.70272e+24]
yt : [INFO     ] 2024-03-06 15:35:05,236 Parameters: domain_right_edge         = [3.70272e+24 3.70272e+24 3.70272e+24]
yt : [INFO     ] 2024-03-06 15:35:05,236 Parameters: cosmological_simulation   = 0

Let’s use yt to take a slice of density and temperature through the center of the dataset so we can see what we’re looking at:

[3]:
slc = yt.SlicePlot(
    ds, "z", [("gas", "density"), ("gas", "temperature")], width=(1.0, "Mpc")
)
slc.show()
yt : [INFO     ] 2024-03-06 15:35:05,719 xlim = -1542838790481162406985728.000000 1542838790481162406985728.000000
yt : [INFO     ] 2024-03-06 15:35:05,720 ylim = -1542838790481162406985728.000000 1542838790481162406985728.000000
yt : [INFO     ] 2024-03-06 15:35:05,721 xlim = -1542838790481162406985728.000000 1542838790481162406985728.000000
yt : [INFO     ] 2024-03-06 15:35:05,721 ylim = -1542838790481162406985728.000000 1542838790481162406985728.000000
yt : [INFO     ] 2024-03-06 15:35:05,726 Making a fixed resolution buffer of (('gas', 'density')) 800 by 800
yt : [INFO     ] 2024-03-06 15:35:06,005 Making a fixed resolution buffer of (('gas', 'temperature')) 800 by 800


Ok, sloshing gas as advertised. Next, we’ll create a sphere object to serve as a source for our investigations. Place it at the center of the domain with "c", and use a radius of 500 kpc:

[4]:
sp = ds.sphere("c", (500.0, "kpc"))

Now, we need to set up a source model. We said we were going to look at the thermal emission from the hot plasma, which in a galaxy cluster is in collisional ionization equilibrium (CIE), so to do that we can set up a CIESourceModel. The first argument specifies which model we want to use. Let’s use "spex". The next three arguments are the maximum and minimum energies, and the number of bins in the spectrum. We’ve chosen these numbers so that the spectrum has an energy resolution of about 1 eV.

CIESourceModel takes a lot of optional arguments, which you can investigate in the docs, but here we’ll do something simple and say that the metallicity is a constant \(Z = 0.3~Z_\odot\):

[5]:
source_model = pyxsim.CIESourceModel("spex", 0.05, 11.0, 1000, 0.3, binscale="log")
pyxsim : [INFO     ] 2024-03-06 15:35:07,167 kT_min = 0.025 keV
pyxsim : [INFO     ] 2024-03-06 15:35:07,168 kT_max = 64 keV

We can use this source_model object to compute X-ray fields for use in yt calculations. For example, we can create fields for emissivity, luminosity, and photon emissivity within a particular band, in this case 0.5-7 keV:

[6]:
xray_fields = source_model.make_source_fields(ds, 0.5, 7.0)
print(xray_fields)
pyxsim : [INFO     ] 2024-03-06 15:35:07,201 Using emission measure field '('gas', 'emission_measure')'.
pyxsim : [INFO     ] 2024-03-06 15:35:07,202 Using temperature field '('gas', 'temperature')'.
[('gas', 'xray_emissivity_0.5_7.0_keV'), ('gas', 'xray_luminosity_0.5_7.0_keV'), ('gas', 'xray_photon_emissivity_0.5_7.0_keV'), ('gas', 'xray_photon_count_rate_0.5_7.0_keV')]

From this we can compute the total luminosity in our sphere object:

[7]:
print(sp["gas", "xray_luminosity_0.5_7.0_keV"])
print(sp.sum(("gas", "xray_luminosity_0.5_7.0_keV")))
[8.08160771e+37 7.91112265e+37 7.94028711e+37 ... 1.36809959e+38
 1.37819810e+38 1.38508499e+38] erg/s
6.782865176868491e+44 erg/s

and we can make a projection of the emissivity:

[8]:
prj = yt.ProjectionPlot(
    ds, "z", ("gas", "xray_emissivity_0.5_7.0_keV"), width=(1.0, "Mpc")
)
prj.show()
yt : [INFO     ] 2024-03-06 15:36:07,181 Projection completed
yt : [INFO     ] 2024-03-06 15:36:07,182 xlim = -1542838790481162406985728.000000 1542838790481162406985728.000000
yt : [INFO     ] 2024-03-06 15:36:07,182 ylim = -1542838790481162406985728.000000 1542838790481162406985728.000000
yt : [INFO     ] 2024-03-06 15:36:07,183 xlim = -1542838790481162406985728.000000 1542838790481162406985728.000000
yt : [INFO     ] 2024-03-06 15:36:07,183 ylim = -1542838790481162406985728.000000 1542838790481162406985728.000000
yt : [INFO     ] 2024-03-06 15:36:07,184 Making a fixed resolution buffer of (('gas', 'xray_emissivity_0.5_7.0_keV')) 800 by 800

We’re almost ready to go to generate the photons from this source, but first we should decide what our redshift, collecting area, and exposure time should be. Let’s pick big numbers, because remember the point of this first step is to create a Monte-Carlo sample from which to draw smaller sub-samples for mock observations. Note these are all (value, unit) tuples:

[9]:
exp_time = (300.0, "ks")  # exposure time
area = (1000.0, "cm**2")  # collecting area
redshift = 0.05

So, that’s everything–let’s create the photons! We use the make_photons function for this:

[10]:
n_photons, n_cells = pyxsim.make_photons(
    "sloshing_photons", sp, redshift, area, exp_time, source_model
)
pyxsim : [INFO     ] 2024-03-06 15:36:07,489 Cosmology: h = 0.71, omega_matter = 0.27, omega_lambda = 0.73
pyxsim : [INFO     ] 2024-03-06 15:36:07,490 Using emission measure field '('gas', 'emission_measure')'.
pyxsim : [INFO     ] 2024-03-06 15:36:07,490 Using temperature field '('gas', 'temperature')'.
pyxsim : [INFO     ] 2024-03-06 15:37:45,605 Finished generating photons.
pyxsim : [INFO     ] 2024-03-06 15:37:45,605 Number of photons generated: 41406462
pyxsim : [INFO     ] 2024-03-06 15:37:45,608 Number of cells with photons: 4375313

Ok, that was easy. Now we have a photon list that we can use to create events using the project_photons function. Here, we’ll just do a simple projection along the z-axis, and center the photons at RA, Dec = (45, 30) degrees. Since we want to be realistic, we’ll want to apply foreground galactic absorption using the "tbabs" model, assuming a neutral hydrogen column of \(N_H = 4 \times 10^{20}~{\rm cm}^{-2}\):

[11]:
n_events = pyxsim.project_photons(
    "sloshing_photons",
    "sloshing_events",
    "z",
    (45.0, 30.0),
    absorb_model="tbabs",
    nH=0.04,
)
pyxsim : [INFO     ] 2024-03-06 15:37:45,620 Foreground galactic absorption: using the tbabs model and nH = 0.04.
pyxsim : [INFO     ] 2024-03-06 15:37:56,276 Detected 15389490 events.

Now that we have a set of “events” on the sky, we can use them as an input to the instrument simulator in SOXS. We’ll use a small exposure time (100 ks instead of 300 ks), and observe it with the as-launched ACIS-I model:

[12]:
soxs.instrument_simulator(
    "sloshing_events.h5",
    "evt.fits",
    (100.0, "ks"),
    "chandra_acisi_cy0",
    [45.0, 30.0],
    overwrite=True,
)
soxs : [INFO     ] 2024-03-06 15:37:56,281 Making observation of source in evt.fits.
soxs : [INFO     ] 2024-03-06 15:37:56,714 Detecting events from source sloshing_events.
soxs : [INFO     ] 2024-03-06 15:37:56,714 Applying energy-dependent effective area from acisi_aimpt_cy0.arf.
soxs : [INFO     ] 2024-03-06 15:37:57,688 Pixeling events.
soxs : [INFO     ] 2024-03-06 15:37:57,877 Scattering events with a multi_image-based PSF.
soxs : [INFO     ] 2024-03-06 15:37:58,116 1133856 events were detected from the source.
soxs : [INFO     ] 2024-03-06 15:37:58,158 Scattering energies with RMF acisi_aimpt_cy0.rmf.
soxs : [INFO     ] 2024-03-06 15:37:58,669 Adding background events.
soxs : [INFO     ] 2024-03-06 15:37:58,702 Adding in point-source background.
soxs : [INFO     ] 2024-03-06 15:37:59,057 Detecting events from source ptsrc_bkgnd.
soxs : [INFO     ] 2024-03-06 15:37:59,057 Applying energy-dependent effective area from acisi_aimpt_cy0.arf.
soxs : [INFO     ] 2024-03-06 15:37:59,062 Pixeling events.
soxs : [INFO     ] 2024-03-06 15:37:59,065 Scattering events with a multi_image-based PSF.
soxs : [INFO     ] 2024-03-06 15:37:59,097 13592 events were detected from the source.
soxs : [INFO     ] 2024-03-06 15:37:59,098 Scattering energies with RMF acisi_aimpt_cy0.rmf.
soxs : [INFO     ] 2024-03-06 15:37:59,204 Generated 13592 photons from the point-source background.
soxs : [INFO     ] 2024-03-06 15:37:59,205 Adding in astrophysical foreground.
soxs : [INFO     ] 2024-03-06 15:37:59,713 Adding in instrumental background.
soxs : [INFO     ] 2024-03-06 15:37:59,741 Making 4282 events from the galactic foreground.
soxs : [INFO     ] 2024-03-06 15:37:59,742 Making 113578 events from the instrumental background.
soxs : [INFO     ] 2024-03-06 15:37:59,795 Observation complete.
soxs : [INFO     ] 2024-03-06 15:37:59,798 Writing events to file evt.fits.
soxs : [INFO     ] 2024-03-06 15:37:59,917 Observation complete.

We can use the write_image() function in SOXS to bin the events into an image and write them to a file, restricting the energies between 0.5 and 2.0 keV:

[13]:
soxs.write_image("evt.fits", "img.fits", emin=0.5, emax=2.0, overwrite=True)

Now we can take a quick look at the image:

[14]:
soxs.plot_image(
    "img.fits", stretch="sqrt", cmap="arbre", vmin=0.0, vmax=10.0, width=0.2
)
[14]:
(<Figure size 1000x1000 with 2 Axes>, <WCSAxes: >)
../_images/cookbook_Thermal_Emission_28_1.png