More Advanced Thermal Emission and Internal Absorption¶
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. We’ll also show how to use a column density map to apply internal absorption to the photons.
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
import numpy as np
from astropy.io import fits
from astropy.convolution import Gaussian2DKernel
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[pfilter.filtered_type, "star_formation_rate"] == 0.0
pfilter3 = data[pfilter.filtered_type, "density"] < 3e-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_31_rotated.hdf5",
bounding_box=[[-1000.0, 1000], [-1000.0, 1000], [-1000.0, 1000]],
)
ds.add_particle_filter("hot_gas")
yt : [INFO ] 2026-04-13 11:58:15,164 Calculating time from 1.000e+00 to be 4.356e+17 seconds
yt : [INFO ] 2026-04-13 11:58:15,190 Parameters: current_time = 4.355810528213311e+17 s
yt : [INFO ] 2026-04-13 11:58:15,191 Parameters: domain_dimensions = [1 1 1]
yt : [INFO ] 2026-04-13 11:58:15,191 Parameters: domain_left_edge = [-1000. -1000. -1000.]
yt : [INFO ] 2026-04-13 11:58:15,191 Parameters: domain_right_edge = [1000. 1000. 1000.]
yt : [INFO ] 2026-04-13 11:58:15,191 Parameters: cosmological_simulation = True
yt : [INFO ] 2026-04-13 11:58:15,191 Parameters: current_redshift = 2.220446049250313e-16
yt : [INFO ] 2026-04-13 11:58:15,192 Parameters: omega_lambda = 0.6911
yt : [INFO ] 2026-04-13 11:58:15,192 Parameters: omega_matter = 0.3089
yt : [INFO ] 2026-04-13 11:58:15,192 Parameters: omega_radiation = 0.0
yt : [INFO ] 2026-04-13 11:58:15,192 Parameters: hubble_constant = 0.6774
yt : [WARNING ] 2026-04-13 11:58:15,192 A bounding box was explicitly specified, so we are disabling periodicity.
yt : [INFO ] 2026-04-13 11:58:15,216 Allocating for 2.734e+07 particles
Loading particle index: 100%|███████████████████████████████████████████████████| 63/63 [00:00<00:00, 3163.88it/s]
[3]:
True
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.
[4]:
source_model = pyxsim.PionSourceModel(
0.3,
1.4,
4000,
("hot_gas", "metallicity"),
binscale="log",
temperature_field=("hot_gas", "temperature"),
emission_measure_field=("hot_gas", "emission_measure"),
nh_field=("hot_gas", "H_nuclei_density"),
)
pyxsim : [INFO ] 2026-04-13 11:58:16,891 kT_min = 0.025 keV
pyxsim : [INFO ] 2026-04-13 11:58:16,892 kT_max = 64 keV
As in other examples, we choose big numbers for the collecting area and exposure time, and a redshift:
[5]:
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:
[6]:
c = ds.arr([0.0, 0.0, 0.0], "code_length")
width = ds.quan(0.5, "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:
[7]:
n_photons, n_cells = pyxsim.make_photons(
"cutout_31_photons", box, redshift, area, exp_time, source_model
)
pyxsim : [INFO ] 2026-04-13 11:58:16,904 Cosmology: h = 0.6774, omega_matter = 0.3089, omega_lambda = 0.6911
pyxsim : [INFO ] 2026-04-13 11:58:16,905 Using emission measure field '('hot_gas', 'emission_measure')'.
pyxsim : [INFO ] 2026-04-13 11:58:16,905 Using temperature field '('hot_gas', 'temperature')'.
pyxsim : [INFO ] 2026-04-13 11:58:16,905 Using nH field '('hot_gas', 'H_nuclei_density')'.
pyxsim : [INFO ] 2026-04-13 12:00:37,082 Finished generating photons.
pyxsim : [INFO ] 2026-04-13 12:00:37,082 Number of photons generated: 1296906
pyxsim : [INFO ] 2026-04-13 12:00:37,083 Number of cells with photons: 523954
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}\):
[8]:
theta = -40.0 * np.pi / 180.0
N = [0.0, 0.0, 1.0]
L = [np.cos(theta), 0.0, np.sin(theta)]
[9]:
n_events = pyxsim.project_photons(
"cutout_31_photons",
"cutout_31_events",
L,
(30.0, 45.0),
absorb_model="tbabs",
nH=0.02,
north_vector=N,
prng=24,
)
pyxsim : [INFO ] 2026-04-13 12:00:37,089 Foreground galactic absorption: using the tbabs model and nH = 0.02.
pyxsim : [INFO ] 2026-04-13 12:00:37,467 Detected 1026503 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 observe it with the 2eV LEM model for 1 Ms. First, we’ll create a background file that we’ll use for the background:
[10]:
soxs.make_background_file(
"bkgnd_evt_31.fits", (1000.0, "ks"), "lem_2eV", [30.0, 45.0], overwrite=True, prng=23,
)
soxs : [INFO ] 2026-04-13 12:00:37,646 Adding in point-source background.
soxs : [INFO ] 2026-04-13 12:00:52,480 Simulating events from 1 sources using instrument lem_2eV for 1000 ks.
soxs : [INFO ] 2026-04-13 12:00:55,616 Scattering energies with RMF lem_2ev_110422.rmf.
soxs : [INFO ] 2026-04-13 12:01:05,007 Detected 3510666 events in total.
soxs : [INFO ] 2026-04-13 12:01:05,029 Generated 3510666 photons from the point-source background.
soxs : [INFO ] 2026-04-13 12:01:05,029 Adding in astrophysical foreground.
soxs : [INFO ] 2026-04-13 12:01:10,380 Adding in instrumental background.
soxs : [INFO ] 2026-04-13 12:01:10,640 Making 6736536 events from the galactic foreground.
soxs : [INFO ] 2026-04-13 12:01:10,640 Making 3414303 events from the instrumental background.
soxs : [INFO ] 2026-04-13 12:01:12,217 Writing background events to file bkgnd_evt_31.fits.
Now we simulate the source itself, adding in the background:
[11]:
soxs.instrument_simulator(
"cutout_31_events.h5",
"evt_31.fits",
(1000.0, "ks"),
"lem_2eV",
[30.0, 45.0],
overwrite=True,
subpixel_res=True,
bkgnd_file="bkgnd_evt_31.fits",
prng=25,
)
soxs : [INFO ] 2026-04-13 12:01:13,089 Simulating events from 1 sources using instrument lem_2eV for 1000 ks.
soxs : [INFO ] 2026-04-13 12:01:13,243 Scattering energies with RMF lem_2ev_110422.rmf.
soxs : [INFO ] 2026-04-13 12:01:13,836 Detected 295971 events in total.
soxs : [INFO ] 2026-04-13 12:01:13,840 Adding background events from the file bkgnd_evt_31.fits.
soxs : [INFO ] 2026-04-13 12:01:15,282 Adding 13661505 background events from bkgnd_evt_31.fits.
soxs : [INFO ] 2026-04-13 12:01:16,443 Observation complete.
soxs : [INFO ] 2026-04-13 12:01:16,465 Writing events to file evt_31.fits.
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:
[12]:
soxs.write_image("evt_31.fits", "img_31.fits", emin=0.644, emax=0.65, reblock=0.5, overwrite=True)
Now we can take a quick look at the image:
[13]:
soxs.plot_image(
"img_31.fits", stretch="log", cmap="inferno", width=0.4, vmin=0.1, smoothing_kernel=Gaussian2DKernel(0.75)
)
[13]:
(<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:
[14]:
soxs.filter_events(
"evt_31.fits",
"evt_31_filter.fits",
overwrite=True,
region='fk5\ncircle(30.0000000,45.0000000,540.000")',
)
soxs.filter_events(
"bkgnd_evt_31.fits",
"bkgnd_evt_31_filter.fits",
overwrite=True,
region='fk5\ncircle(30.0000000,45.0000000,540.000")',
)
Now bin up spectra for these new event files:
[15]:
soxs.write_spectrum("evt_31_filter.fits", "evt_31.pi", overwrite=True)
soxs.write_spectrum("bkgnd_evt_31_filter.fits", "bkgnd_evt_31.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.
[16]:
fig, ax, _ = soxs.plot_spectrum(
"evt_31.pi", xmin=0.5, xmax=0.7, xscale="linear", yscale="log", color="C0", ymin=0.5,
plot_steps=True
)
soxs.plot_spectrum(
"bkgnd_evt_31.pi",
xmin=0.5,
xmax=0.7,
fig=fig,
ax=ax,
yscale="log",
color="C1",
ymin=0.4,
plot_steps=True,
)
[16]:
(<Figure size 1000x1000 with 1 Axes>,
<Axes: xlabel='Energy (keV)', ylabel='Count Rate (counts/s/keV)'>,
array([0.1 , 0.10025, 0.1005 , ..., 3.0995 , 3.09975, 3.1 ],
shape=(12001,)))
Including Absorption by Neutral Gas Internal to the Source¶
We can also use this dataset to explore the effects of taking into account absorption of the photons by the neutral gas in the source itself. To do this, we need to create a column density map of the neutral hydrogen in the source as a function of depth, and then use that as an input to the project_photons function. The column density map is created using the make_column_density_map() function in pyXSIM. We will use the same normal and north vectors as before, and we will save the column
density map to a file. The field we use for the column density is the number density of neutral hydrogen, which is given by the "H_p0_number_density" field in the dataset by default, but an alternative field may be used.
[17]:
width = (0.5, "Mpc") # width of the column density map
depth = (0.5, "Mpc") # depth of the column density map
nx = 128 # number of pixels in the x and y directions (width)
nz = 128 # number of pixels in the z direction (depth)
column_file = "nH.h5" # name of the file to save the column density map to
pyxsim.make_column_density_map(
ds,
L,
c,
width,
depth,
nx,
nz,
column_file,
field=("gas", "H_p0_number_density"),
north_vector=N,
)
Having created this map of column density, we can now use it in the project_photons() function to create a new set of events that includes the effects of absorption by the neutral gas in the source itself. We will use the same foreground absorption model and Galactic column density as before, but now we will also include the internal absorption using the column density map we just created:
[18]:
n_events_absorb = pyxsim.project_photons(
"cutout_31_photons",
"cutout_31_nH_events",
L,
(30.0, 45.0),
absorb_model="tbabs",
nH=0.02,
north_vector=N,
column_file=column_file,
prng=24,
)
pyxsim : [INFO ] 2026-04-13 12:05:14,795 Foreground galactic absorption: using the tbabs model and nH = 0.02.
pyxsim : [INFO ] 2026-04-13 12:05:14,795 Internal absorption: using the tbabs model and column density map nH.h5.
pyxsim : [INFO ] 2026-04-13 12:05:18,593 Detected 896182 events.
We notice that ~20% of the emitted photons are absorbed by the neutral gas in the source itself:
[19]:
print((n_events_absorb-n_events)/n_events)
-0.12695627777025492
Now we can simulate the observation of this new set of events with internal absorption using the instrument simulator in SOXS, using the same parameters as before:
[20]:
soxs.instrument_simulator(
"cutout_31_nH_events.h5",
"evt_nH_31.fits",
(1000.0, "ks"),
"lem_2eV",
[30.0, 45.0],
overwrite=True,
subpixel_res=True,
bkgnd_file="bkgnd_evt_31.fits",
prng=25,
)
soxs : [INFO ] 2026-04-13 12:05:18,771 Simulating events from 1 sources using instrument lem_2eV for 1000 ks.
soxs : [INFO ] 2026-04-13 12:05:18,867 Scattering energies with RMF lem_2ev_110422.rmf.
soxs : [INFO ] 2026-04-13 12:05:19,388 Detected 263881 events in total.
soxs : [INFO ] 2026-04-13 12:05:19,392 Adding background events from the file bkgnd_evt_31.fits.
soxs : [INFO ] 2026-04-13 12:05:20,898 Adding 13661505 background events from bkgnd_evt_31.fits.
soxs : [INFO ] 2026-04-13 12:05:22,043 Observation complete.
soxs : [INFO ] 2026-04-13 12:05:22,060 Writing events to file evt_nH_31.fits.
We can write an image of this new set of events with internal absorption, and then compare it to the image we made before without internal absorption. We will write the new image to a file, and then read in both images and create a difference image:
[21]:
soxs.write_image("evt_nH_31.fits", "img_nH_31.fits", emin=0.644, emax=0.65, reblock=0.5, overwrite=True)
[22]:
with fits.open("img_31.fits") as i1:
img1 = i1[0].data
header = i1[0].header
with fits.open("img_nH_31.fits") as i2:
img2 = i2[0].data
fits.ImageHDU(img1-img2, header=header).writeto("img_diff.fits", overwrite=True)
[23]:
soxs.plot_image("img_diff.fits", stretch="log", cmap="inferno", width=0.4, vmin=0.1)
[23]:
(<Figure size 1000x1000 with 2 Axes>, <WCSAxes: >)
