More Advanced Thermal Emission¶
In this example, we’ll look at the emission from a disk galaxy from the Illustris TNG simulations. This dataset has metallicity information for several species in it. We’ll make a cut in phase space like we did in the previous example. The dataset we want to use for this example is available for download here.
First, import our necessary modules:
[1]:
import yt
import pyxsim
import soxs
We will make phase space cuts on the gas cells using density, temperature, and star formation rate:
[2]:
# Note that the units of all numbers in this function are CGS
# define hot gas filter
def hot_gas(pfilter, data):
pfilter1 = data[pfilter.filtered_type, "temperature"] > 3.0e5
pfilter2 = data["PartType0", "StarFormationRate"] == 0.0
pfilter3 = data[pfilter.filtered_type, "density"] < 5e-25
return pfilter1 & pfilter2 & pfilter3
yt.add_particle_filter(
"hot_gas",
function=hot_gas,
filtered_type="gas",
requires=["temperature", "density"],
)
Next, we load
the dataset with yt, and add the "hot_gas"
filter to the dataset:
[3]:
ds = yt.load(
"cutout_37.hdf5", bounding_box=[[-1000.0, 1000], [-1000.0, 1000], [-1000.0, 1000]]
)
ds.add_particle_filter("hot_gas")
yt : [INFO ] 2023-01-10 14:13:02,289 Calculating time from 1.000e+00 to be 4.356e+17 seconds
yt : [INFO ] 2023-01-10 14:13:02,326 Parameters: current_time = 4.355810528213311e+17 s
yt : [INFO ] 2023-01-10 14:13:02,327 Parameters: domain_dimensions = [1 1 1]
yt : [INFO ] 2023-01-10 14:13:02,327 Parameters: domain_left_edge = [-1000. -1000. -1000.]
yt : [INFO ] 2023-01-10 14:13:02,327 Parameters: domain_right_edge = [1000. 1000. 1000.]
yt : [INFO ] 2023-01-10 14:13:02,327 Parameters: cosmological_simulation = 1
yt : [INFO ] 2023-01-10 14:13:02,327 Parameters: current_redshift = 2.220446049250313e-16
yt : [INFO ] 2023-01-10 14:13:02,328 Parameters: omega_lambda = 0.6911
yt : [INFO ] 2023-01-10 14:13:02,329 Parameters: omega_matter = 0.3089
yt : [INFO ] 2023-01-10 14:13:02,329 Parameters: omega_radiation = 0.0
yt : [INFO ] 2023-01-10 14:13:02,329 Parameters: hubble_constant = 0.6774
yt : [WARNING ] 2023-01-10 14:13:02,329 A bounding box was explicitly specified, so we are disabling periodicity.
yt : [INFO ] 2023-01-10 14:13:02,339 Allocating for 2.761e+06 particles
Initializing coarse index : 100%|███████████████| 11/11 [00:00<00:00, 62.34it/s]
yt : [INFO ] 2023-01-10 14:13:02,528 Updating index_order2 from 2 to 2
Initializing refined index: 100%|███████████████| 11/11 [00:00<00:00, 12.83it/s]
[3]:
True
We also need to tell pyXSIM which elements have fields in the dataset that should be used. To do this we create a var_elem
dictionary of (key, value) pairs corresponding to the element name and the yt field name (assuming the "hot_gas"
type).
[4]:
# metal fields to use
metals = [
"C_fraction",
"N_fraction",
"O_fraction",
"Ne_fraction",
"Mg_fraction",
"Si_fraction",
"Fe_fraction",
]
var_elem = {elem.split("_")[0]: ("hot_gas", elem) for elem in metals}
Now that we have everything we need, we’ll set up the IGMSourceModel
, which is based on Cloudy and includes resonant scattering off of the CXB (see here for more details). Because we created a hot gas filter, we will use the "hot_gas"
field type for the emission measure, temperature, and metallicity fields.
[5]:
source_model = pyxsim.IGMSourceModel(
0.1,
4.0,
5000,
("hot_gas", "metallicity"),
binscale="log",
resonant_scattering=True,
temperature_field=("hot_gas", "temperature"),
emission_measure_field=("hot_gas", "emission_measure"),
nh_field=("hot_gas", "H_nuclei_density"),
var_elem=var_elem,
)
pyxsim : [INFO ] 2023-01-10 14:13:03,849 kT_min = 0.00431 keV
pyxsim : [INFO ] 2023-01-10 14:13:03,849 kT_max = 64 keV
As in other examples, we choose big numbers for the collecting area and exposure time, and a redshift:
[6]:
exp_time = (1.0, "Ms") # exposure time
area = (5000.0, "cm**2") # collecting area
redshift = 0.01
Next, we’ll create a box object to serve as a source for the photons. The dataset consists of only the galaxy at a specific location, which we use below, and pick a width of 1 Mpc:
[7]:
c = ds.arr([0.0, 0.0, 0.0], "code_length")
width = ds.quan(1.0, "Mpc")
le = c - 0.5 * width
re = c + 0.5 * width
box = ds.box(le, re)
So, that’s everything–let’s create the photons! We use the make_photons
function for this:
[8]:
n_photons, n_cells = pyxsim.make_photons(
"cutout_37_photons", box, redshift, area, exp_time, source_model
)
pyxsim : [INFO ] 2023-01-10 14:13:03,861 Cosmology: h = 0.6774, omega_matter = 0.3089, omega_lambda = 0.6911
pyxsim : [INFO ] 2023-01-10 14:13:03,862 Using emission measure field '('hot_gas', 'emission_measure')'.
pyxsim : [INFO ] 2023-01-10 14:13:03,862 Using temperature field '('hot_gas', 'temperature')'.
pyxsim : [INFO ] 2023-01-10 14:13:03,862 Using nH field '('hot_gas', 'H_nuclei_density')'.
pyxsim : [INFO ] 2023-01-10 14:32:33,004 Finished generating photons.
pyxsim : [INFO ] 2023-01-10 14:32:33,012 Number of photons generated: 1859222
pyxsim : [INFO ] 2023-01-10 14:32:33,013 Number of cells with photons: 563639
And now we create events using the project_photons
function. Let’s project along the "z"
axis. We’ll use the "tbabs"
foreground absorption model this time, with a neutral hydrogen column of \(N_H = 2 \times 10^{20}~{\rm cm}^{-2}\):
[9]:
n_events = pyxsim.project_photons(
"cutout_37_photons",
"cutout_37_events",
"x",
(30.0, 45.0),
absorb_model="tbabs",
nH=0.02,
)
pyxsim : [INFO ] 2023-01-10 14:32:33,088 Foreground galactic absorption: using the tbabs model and nH = 0.02.
pyxsim : [INFO ] 2023-01-10 14:32:33,514 Detected 374368 events.
Now that we have a set of “events” on the sky, we can read them in and write them to a SIMPUT file:
[10]:
events = pyxsim.EventList("cutout_37_events.h5")
events.write_to_simput("cutout_37", overwrite=True)
soxs : [WARNING ] 2023-01-10 14:32:33,556 Overwriting cutout_37_phlist.fits.
soxs : [INFO ] 2023-01-10 14:32:33,576 Writing source 'cutout_37' to cutout_37_phlist.fits.
We can then use this SIMPUT file as an input to the instrument simulator in SOXS. We’ll observe it with the 2eV LEM model for 1 Ms. First, we’ll create a background file that we’ll use for the background:
[11]:
soxs.make_background_file(
"bkgnd_evt_37.fits", (1000.0, "ks"), "lem_2eV", [30.0, 45.0], overwrite=True
)
soxs : [INFO ] 2023-01-10 14:32:33,733 Adding in point-source background.
soxs : [INFO ] 2023-01-10 14:32:42,064 Detecting events from source ptsrc_bkgnd.
soxs : [INFO ] 2023-01-10 14:32:42,064 Applying energy-dependent effective area from lem_300522.arf.
soxs : [INFO ] 2023-01-10 14:32:42,961 Pixeling events.
soxs : [INFO ] 2023-01-10 14:32:43,491 Scattering events with a gaussian-based PSF.
soxs : [INFO ] 2023-01-10 14:32:43,716 3525153 events were detected from the source.
soxs : [INFO ] 2023-01-10 14:32:43,832 Scattering energies with RMF lem_2ev_110422.rmf.
soxs : [INFO ] 2023-01-10 14:33:11,589 Generated 3525153 photons from the point-source background.
soxs : [INFO ] 2023-01-10 14:33:11,589 Adding in astrophysical foreground.
soxs : [INFO ] 2023-01-10 14:33:17,413 Adding in instrumental background.
soxs : [INFO ] 2023-01-10 14:33:17,770 Making 6081885 events from the galactic foreground.
soxs : [INFO ] 2023-01-10 14:33:17,772 Making 3412450 events from the instrumental background.
soxs : [INFO ] 2023-01-10 14:33:18,570 Writing events to file bkgnd_evt_37.fits.
Now we simulate the source itself, adding in the background:
[12]:
soxs.instrument_simulator(
"cutout_37_simput.fits",
"evt_37.fits",
(1000.0, "ks"),
"lem_2eV",
[30.0, 45.0],
overwrite=True,
bkgnd_file="bkgnd_evt_37.fits",
)
soxs : [INFO ] 2023-01-10 14:33:20,827 Making observation of source in evt_37.fits.
soxs : [INFO ] 2023-01-10 14:33:20,958 Detecting events from source cutout_37.
soxs : [INFO ] 2023-01-10 14:33:20,960 Applying energy-dependent effective area from lem_300522.arf.
soxs : [INFO ] 2023-01-10 14:33:20,990 Pixeling events.
soxs : [INFO ] 2023-01-10 14:33:21,005 Scattering events with a gaussian-based PSF.
soxs : [INFO ] 2023-01-10 14:33:21,011 79988 events were detected from the source.
soxs : [INFO ] 2023-01-10 14:33:21,015 Scattering energies with RMF lem_2ev_110422.rmf.
soxs : [INFO ] 2023-01-10 14:33:21,948 Adding background events from the file bkgnd_evt_37.fits.
soxs : [INFO ] 2023-01-10 14:33:22,797 Adding 13019488 background events from bkgnd_evt_37.fits.
soxs : [INFO ] 2023-01-10 14:33:23,224 Writing events to file evt_37.fits.
soxs : [INFO ] 2023-01-10 14:33:26,554 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.644 and 0.65 keV, which focuses on the redshifted OVIII line:
[13]:
soxs.write_image("evt_37.fits", "img_37.fits", emin=0.644, emax=0.65, overwrite=True)
Now we can take a quick look at the image:
[14]:
soxs.plot_image("img_37.fits", stretch="log", cmap="arbre", width=0.4, vmin=0.5)
[14]:
(<Figure size 1000x1000 with 2 Axes>, <WCSAxes: >)

Now we will make spectra to look at. First, filter the events of both the combined source and background files and the background-only files within 0.15 degree of the center:
[15]:
soxs.filter_events(
"evt_37.fits",
"evt_37_filter.fits",
overwrite=True,
region='fk5\ncircle(30.0000000,45.0000000,540.000")',
)
soxs.filter_events(
"bkgnd_evt_37.fits",
"bkgnd_evt_37_filter.fits",
overwrite=True,
region='fk5\ncircle(30.0000000,45.0000000,540.000")',
)
Now bin up spectra for these new event files:
[16]:
soxs.write_spectrum("evt_37_filter.fits", "evt_37.pi", overwrite=True)
soxs.write_spectrum("bkgnd_evt_37_filter.fits", "bkgnd_evt_37.pi", overwrite=True)
Finally, we can plot the spectra. Below, the total spectrum is in blue and the background/foreground spectrum is in orange. The lines from the emission of the distant galaxy are redshifted away from the foreground Milky Way lines.
[17]:
fig, ax = soxs.plot_spectrum("evt_37.pi", xmin=0.5, xmax=0.7, xscale="linear", ymin=0.5)
soxs.plot_spectrum("bkgnd_evt_37.pi", xmin=0.5, xmax=0.7, fig=fig, ax=ax, ymin=0.5)
[17]:
(<Figure size 1000x1000 with 1 Axes>,
<AxesSubplot: xlabel='Energy (keV)', ylabel='Count Rate (counts/s/keV)'>)
