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_31_rotated.hdf5",
bounding_box=[[-1000.0, 1000], [-1000.0, 1000], [-1000.0, 1000]],
)
ds.add_particle_filter("hot_gas")
yt : [INFO ] 2024-03-06 14:56:08,455 Calculating time from 1.000e+00 to be 4.356e+17 seconds
yt : [INFO ] 2024-03-06 14:56:08,491 Parameters: current_time = 4.355810528213311e+17 s
yt : [INFO ] 2024-03-06 14:56:08,491 Parameters: domain_dimensions = [1 1 1]
yt : [INFO ] 2024-03-06 14:56:08,492 Parameters: domain_left_edge = [-1000. -1000. -1000.]
yt : [INFO ] 2024-03-06 14:56:08,492 Parameters: domain_right_edge = [1000. 1000. 1000.]
yt : [INFO ] 2024-03-06 14:56:08,492 Parameters: cosmological_simulation = True
yt : [INFO ] 2024-03-06 14:56:08,493 Parameters: current_redshift = 2.220446049250313e-16
yt : [INFO ] 2024-03-06 14:56:08,493 Parameters: omega_lambda = 0.6911
yt : [INFO ] 2024-03-06 14:56:08,493 Parameters: omega_matter = 0.3089
yt : [INFO ] 2024-03-06 14:56:08,493 Parameters: omega_radiation = 0.0
yt : [INFO ] 2024-03-06 14:56:08,494 Parameters: hubble_constant = 0.6774
yt : [WARNING ] 2024-03-06 14:56:08,494 A bounding box was explicitly specified, so we are disabling periodicity.
yt : [INFO ] 2024-03-06 14:56:08,536 Allocating for 2.734e+07 particles
- Loading particle index: 0%| | 0/63 [00:00<?, ?it/s]
</pre>
- Loading particle index: 0%| | 0/63 [00:00<?, ?it/s]
end{sphinxVerbatim}
Loading particle index: 0%| | 0/63 [00:00<?, ?it/s]
- Loading particle index: 100%|█████████████████| 63/63 [00:00<00:00, 1988.44it/s]
</pre>
- Loading particle index: 100%|█████████████████| 63/63 [00:00<00:00, 1988.44it/s]
end{sphinxVerbatim}
Loading particle index: 100%|█████████████████| 63/63 [00:00<00:00, 1988.44it/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 ] 2024-03-06 14:56:11,813 kT_min = 0.00431 keV
pyxsim : [INFO ] 2024-03-06 14:56:11,815 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_31_photons", box, redshift, area, exp_time, source_model
)
pyxsim : [INFO ] 2024-03-06 14:56:11,843 Cosmology: h = 0.6774, omega_matter = 0.3089, omega_lambda = 0.6911
pyxsim : [INFO ] 2024-03-06 14:56:11,844 Using emission measure field '('hot_gas', 'emission_measure')'.
pyxsim : [INFO ] 2024-03-06 14:56:11,845 Using temperature field '('hot_gas', 'temperature')'.
pyxsim : [INFO ] 2024-03-06 14:56:11,845 Using nH field '('hot_gas', 'H_nuclei_density')'.
pyxsim : [INFO ] 2024-03-06 15:32:09,329 Finished generating photons.
pyxsim : [INFO ] 2024-03-06 15:32:09,334 Number of photons generated: 9281137
pyxsim : [INFO ] 2024-03-06 15:32:09,335 Number of cells with photons: 1593802
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_31_photons",
"cutout_31_events",
"x",
(30.0, 45.0),
absorb_model="tbabs",
nH=0.02,
)
pyxsim : [INFO ] 2024-03-06 15:32:09,396 Foreground galactic absorption: using the tbabs model and nH = 0.02.
pyxsim : [INFO ] 2024-03-06 15:32:11,226 Detected 1829351 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
)
soxs : [INFO ] 2024-03-06 15:32:11,488 Adding in point-source background.
soxs : [INFO ] 2024-03-06 15:32:18,749 Detecting events from source ptsrc_bkgnd.
soxs : [INFO ] 2024-03-06 15:32:18,749 Applying energy-dependent effective area from lem_300522.arf.
soxs : [INFO ] 2024-03-06 15:32:19,506 Pixeling events.
soxs : [INFO ] 2024-03-06 15:32:19,959 Scattering events with a gaussian-based PSF.
soxs : [INFO ] 2024-03-06 15:32:20,083 3032082 events were detected from the source.
soxs : [INFO ] 2024-03-06 15:32:20,174 Scattering energies with RMF lem_2ev_110422.rmf.
soxs : [INFO ] 2024-03-06 15:32:32,719 Generated 3032082 photons from the point-source background.
soxs : [INFO ] 2024-03-06 15:32:32,720 Adding in astrophysical foreground.
soxs : [INFO ] 2024-03-06 15:32:38,638 Adding in instrumental background.
soxs : [INFO ] 2024-03-06 15:32:39,072 Making 6077801 events from the galactic foreground.
soxs : [INFO ] 2024-03-06 15:32:39,073 Making 3413905 events from the instrumental background.
soxs : [INFO ] 2024-03-06 15:32:39,764 Writing 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,
bkgnd_file="bkgnd_evt_31.fits",
)
soxs : [INFO ] 2024-03-06 15:32:42,441 Making observation of source in evt_31.fits.
soxs : [INFO ] 2024-03-06 15:32:42,782 Detecting events from source cutout_31_events.
soxs : [INFO ] 2024-03-06 15:32:42,783 Applying energy-dependent effective area from lem_300522.arf.
soxs : [INFO ] 2024-03-06 15:32:42,920 Pixeling events.
soxs : [INFO ] 2024-03-06 15:32:42,984 Scattering events with a gaussian-based PSF.
soxs : [INFO ] 2024-03-06 15:32:43,000 355422 events were detected from the source.
soxs : [INFO ] 2024-03-06 15:32:43,018 Scattering energies with RMF lem_2ev_110422.rmf.
soxs : [INFO ] 2024-03-06 15:32:44,405 Adding background events from the file bkgnd_evt_31.fits.
soxs : [INFO ] 2024-03-06 15:32:45,041 Adding 12523788 background events from bkgnd_evt_31.fits.
soxs : [INFO ] 2024-03-06 15:32:45,472 Writing events to file evt_31.fits.
soxs : [INFO ] 2024-03-06 15:32:47,521 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:
[12]:
soxs.write_image("evt_31.fits", "img_31.fits", emin=0.644, emax=0.65, overwrite=True)
Now we can take a quick look at the image:
[13]:
soxs.plot_image("img_31.fits", stretch="log", cmap="arbre", width=0.4, vmin=0.5)
[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", ymin=0.5)
soxs.plot_spectrum("bkgnd_evt_31.pi", xmin=0.5, xmax=0.7, fig=fig, ax=ax, ymin=0.5)
[16]:
(<Figure size 1000x1000 with 1 Axes>,
<Axes: xlabel='Energy (keV)', ylabel='Count Rate (counts/s/keV)'>)