Line Emission

This example shows how to create a line emission model. It uses a galaxy cluster from a Gadget SPH cosmological dataset, and will create a thermal model out of the gas particles and will use the dark matter particles to add line emission to the spectrum, assuming that the emission comes from some decay process of the dark matter.

First, we load the modules:

[1]:
%matplotlib inline
import matplotlib

matplotlib.rc("font", size=18, family="serif")
import yt
import matplotlib.pyplot as plt
from yt.units import mp
import pyxsim

As in the Advanced Thermal Emission example, we’ll set up a hot gas filter:

[2]:
# Note that the units of all numbers in this function are CGS
def hot_gas(pfilter, data):
    pfilter1 = data[pfilter.filtered_type, "density"] < 5e-25
    pfilter2 = data[pfilter.filtered_type, "temperature"] > 3481355.78432401
    pfilter3 = data[pfilter.filtered_type, "temperature"] < 4.5e8
    return (pfilter1) & (pfilter2) & (pfilter3)


yt.add_particle_filter(
    "hot_gas",
    function=hot_gas,
    filtered_type="gas",
    requires=["density", "temperature"],
)

The dataset used here does not have a field for the electron number density, which is required to construct the emission measure field. Because we’ll only be using the hot gas, we can create a ("gas","El_number_density") field which assumes complete ionization (while taking into account the H and He mass fractions vary from particle to particle). This is not strictly true for all of the "gas" type particles, but since we’ll be using the "hot_gas" type it should be sufficiently accurate for our purposes. We’ll define the field here and add it.

[3]:
def _El_number_density(field, data):
    mueinv = data["gas", "H_fraction"] + 0.5 * (1.0 - data["gas", "He_fraction"])
    return data["gas", "density"] * mueinv / (1.0 * mp)


yt.add_field(
    ("gas", "El_number_density"),
    _El_number_density,
    units="cm**-3",
    sampling_type="local",
)

and we load the dataset in yt:

[4]:
ds = yt.load("snapshot_033/snap_033.0.hdf5")
ds.add_particle_filter("hot_gas")
yt : [INFO     ] 2024-03-06 14:55:16,920 Parameters: current_time              = 4.343952725460923e+17 s
yt : [INFO     ] 2024-03-06 14:55:16,920 Parameters: domain_dimensions         = [1 1 1]
yt : [INFO     ] 2024-03-06 14:55:16,921 Parameters: domain_left_edge          = [0. 0. 0.]
yt : [INFO     ] 2024-03-06 14:55:16,921 Parameters: domain_right_edge         = [25. 25. 25.]
yt : [INFO     ] 2024-03-06 14:55:16,921 Parameters: cosmological_simulation   = 1
yt : [INFO     ] 2024-03-06 14:55:16,921 Parameters: current_redshift          = -4.811891664902035e-05
yt : [INFO     ] 2024-03-06 14:55:16,922 Parameters: omega_lambda              = 0.762
yt : [INFO     ] 2024-03-06 14:55:16,922 Parameters: omega_matter              = 0.238
yt : [INFO     ] 2024-03-06 14:55:16,922 Parameters: omega_radiation           = 0.0
yt : [INFO     ] 2024-03-06 14:55:16,922 Parameters: hubble_constant           = 0.73
yt : [INFO     ] 2024-03-06 14:55:16,974 Allocating for 4.194e+06 particles
Loading particle index: 0%| | 0/12 [00:00&lt;?, ?it/s]

</pre>

Loading particle index: 0%| | 0/12 [00:00<?, ?it/s]

end{sphinxVerbatim}

Loading particle index: 0%| | 0/12 [00:00<?, ?it/s]

Loading particle index: 100%|██████████████████| 12/12 [00:00&lt;00:00, 772.99it/s]

</pre>

Loading particle index: 100%|██████████████████| 12/12 [00:00<00:00, 772.99it/s]

end{sphinxVerbatim}

Loading particle index: 100%|██████████████████| 12/12 [00:00<00:00, 772.99it/s]


[4]:
True

To get a sense of what the cluster looks like, let’s take a slice through the density and temperature:

[5]:
slc = yt.SlicePlot(
    ds,
    "z",
    [("gas", "density"), ("gas", "temperature")],
    center="max",
    width=(3.0, "Mpc"),
)
slc.show()
yt : [INFO     ] 2024-03-06 14:55:17,891 max value is 7.22543e-22 at 9.1760425567626953 12.2128591537475586 9.3898868560791016
yt : [INFO     ] 2024-03-06 14:55:17,971 xlim = 8.081095 10.270990
yt : [INFO     ] 2024-03-06 14:55:17,971 ylim = 11.117912 13.307806
yt : [INFO     ] 2024-03-06 14:55:17,973 xlim = 8.081095 10.270990
yt : [INFO     ] 2024-03-06 14:55:17,973 ylim = 11.117912 13.307806
yt : [INFO     ] 2024-03-06 14:55:17,974 Making a fixed resolution buffer of (('gas', 'temperature')) 800 by 800
yt : [INFO     ] 2024-03-06 14:55:18,711 Making a fixed resolution buffer of (('gas', 'density')) 800 by 800


Now set up a sphere centered on the maximum density in the dataset:

[6]:
sp = ds.sphere("max", (1.0, "Mpc"))
yt : [INFO     ] 2024-03-06 14:55:20,106 max value is 7.22543e-22 at 9.1760425567626953 12.2128591537475586 9.3898868560791016

and create a thermal source model:

[7]:
thermal_model = pyxsim.CIESourceModel("apec", 0.2, 11.0, 10000, ("gas", "metallicity"))
pyxsim : [INFO     ] 2024-03-06 14:55:20,346 kT_min = 0.025 keV
pyxsim : [INFO     ] 2024-03-06 14:55:20,347 kT_max = 64 keV

Now we’ll set up a line emission field for the dark matter particles. We won’t try to replicate a specific model, but will simply assume that the emission is proportional to the dark matter mass. Note that this field is a particle field.

[8]:
norm = yt.YTQuantity(100.0, "g**-1*s**-1")


def _dm_emission(field, data):
    return data["PartType1", "particle_mass"] * norm


ds.add_field(
    ("PartType1", "dm_emission"),
    function=_dm_emission,
    sampling_type="particle",
    units="photons/s",
    force_override=True,
)

Now we can set up the LineSourceModel object. The first argument is the line center energy in keV, and the second is the field we just created, that sets up the line amplitude. There is another parameter, sigma, for adding in broadening of the line, but in this case we’ll rely on the velocities of the dark matter particles themselves to produce the line broadening.

[9]:
line_model = pyxsim.LineSourceModel(3.5, ("PartType1", "dm_emission"))

Now set up the parameters for generating the photons:

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

and actually generate the photons:

[11]:
ntp, ntc = pyxsim.make_photons(
    "therm_photons.h5", sp, redshift, area, exp_time, thermal_model
)
nlp, nlc = pyxsim.make_photons(
    "line_photons.h5", sp, redshift, area, exp_time, line_model
)
pyxsim : [INFO     ] 2024-03-06 14:55:20,365 Cosmology: h = 0.73, omega_matter = 0.238, omega_lambda = 0.762
pyxsim : [INFO     ] 2024-03-06 14:55:20,365 Using emission measure field '('gas', 'emission_measure')'.
pyxsim : [INFO     ] 2024-03-06 14:55:20,365 Using temperature field '('gas', 'temperature')'.
pyxsim : [INFO     ] 2024-03-06 14:56:00,557 Finished generating photons.
pyxsim : [INFO     ] 2024-03-06 14:56:00,557 Number of photons generated: 1532702
pyxsim : [INFO     ] 2024-03-06 14:56:00,558 Number of cells with photons: 34507
pyxsim : [INFO     ] 2024-03-06 14:56:00,560 Cosmology: h = 0.73, omega_matter = 0.238, omega_lambda = 0.762
---------------------------------------------------------------------------
YTFieldNotParseable                       Traceback (most recent call last)
Cell In[11], line 4
      1 ntp, ntc = pyxsim.make_photons(
      2     "therm_photons.h5", sp, redshift, area, exp_time, thermal_model
      3 )
----> 4 nlp, nlc = pyxsim.make_photons(
      5     "line_photons.h5", sp, redshift, area, exp_time, line_model
      6 )

File ~/Source/pyxsim/pyxsim/photon_list.py:422, in make_photons(photon_prefix, data_source, redshift, area, exp_time, source_model, point_sources, parameters, center, dist, cosmology, velocity_fields, bulk_velocity, observer, fields_to_keep)
    420 for i, ax in enumerate("xyz"):
    421     pos = chunk[p_fields[i]][idxs].to_value("kpc")
--> 422     dx = chunk[w_field][idxs].to_value("kpc")
    423     # Fix photon coordinates for regions crossing a periodic boundary
    424     if ds.periodicity[i]:

File ~/Source/yt/yt/data_objects/data_containers.py:229, in YTDataContainer.__getitem__(self, key)
    225 def __getitem__(self, key):
    226     """
    227     Returns a single field.  Will add if necessary.
    228     """
--> 229     f = self._determine_fields([key])[0]
    230     if f not in self.field_data and key not in self.field_data:
    231         if f in self._container_fields:

File ~/Source/yt/yt/data_objects/data_containers.py:1472, in YTDataContainer._determine_fields(self, fields)
   1469     explicit_fields.append(field)
   1470     continue
-> 1472 finfo = self.ds._get_field_info(field)
   1473 ftype, fname = finfo.name
   1474 # really ugly check to ensure that this field really does exist somewhere,
   1475 # in some naming convention, before returning it as a possible field type

File ~/Source/yt/yt/data_objects/static_output.py:957, in Dataset._get_field_info(self, field)
    952 def _get_field_info(
    953     self,
    954     field: Union[FieldKey, ImplicitFieldKey, DerivedField],
    955     /,
    956 ) -> DerivedField:
--> 957     field_info, candidates = self._get_field_info_helper(field)
    959     if field_info.name[1] in ("px", "py", "pz", "pdx", "pdy", "pdz"):
    960         # escape early as a bandaid solution to
    961         # https://github.com/yt-project/yt/issues/3381
    962         return field_info

File ~/Source/yt/yt/data_objects/static_output.py:1029, in Dataset._get_field_info_helper(self, field)
   1027     ftype, fname = field.name
   1028 else:
-> 1029     raise YTFieldNotParseable(field)
   1031 if ftype == "unknown":
   1032     candidates: list[FieldKey] = [
   1033         (ft, fn) for ft, fn in self.field_info if fn == fname
   1034     ]

YTFieldNotParseable: Cannot identify field None

Next, project the photons for the total set and the line set by itself:

[12]:
nte = pyxsim.project_photons(
    "therm_photons.h5",
    "therm_events.h5",
    "y",
    (30.0, 45.0),
    absorb_model="wabs",
    nH=0.02,
)
nle = pyxsim.project_photons(
    "line_photons.h5", "line_events.h5", "y", (30.0, 45.0), absorb_model="wabs", nH=0.02
)
pyxsim : [INFO     ] 2024-03-06 14:56:00,781 Foreground galactic absorption: using the wabs model and nH = 0.02.
pyxsim : [INFO     ] 2024-03-06 14:56:01,125 Detected 1165454 events.
pyxsim : [INFO     ] 2024-03-06 14:56:01,126 Foreground galactic absorption: using the wabs model and nH = 0.02.
pyxsim : [INFO     ] 2024-03-06 14:56:01,151 Detected 0 events.

Write the raw, unconvolved spectra to disk:

[13]:
et = pyxsim.EventList("therm_events.h5")
el = pyxsim.EventList("line_events.h5")
et.write_spectrum("therm_spec.fits", 0.1, 10.0, 5000, overwrite=True)
el.write_spectrum("line_spec.fits", 0.1, 10.0, 5000, overwrite=True)
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Cell In[13], line 4
      2 el = pyxsim.EventList("line_events.h5")
      3 et.write_spectrum("therm_spec.fits", 0.1, 10.0, 5000, overwrite=True)
----> 4 el.write_spectrum("line_spec.fits", 0.1, 10.0, 5000, overwrite=True)

File ~/Source/pyxsim/pyxsim/event_list.py:350, in EventList.write_spectrum(self, specfile, emin, emax, nchan, overwrite)
    347 col2 = fits.Column(name="ENERGY", format="1D", array=emid.astype("float64"))
    348 col3 = fits.Column(name="COUNTS", format="1J", array=spec.astype("int32"))
    349 col4 = fits.Column(
--> 350     name="COUNT_RATE", format="1D", array=spec / self.parameters["exp_time"]
    351 )
    353 coldefs = fits.ColDefs([col1, col2, col3, col4])
    355 tbhdu = fits.BinTableHDU.from_columns(coldefs)

KeyError: 'exp_time'

Now let’s plot up both spectra. We see that we have a thermal spectrum with the addition of a line at 3.5 keV (in real life such a line would not be so prominent, but it makes the example easier to see):

[14]:
import astropy.io.fits as pyfits

f1 = pyfits.open("therm_spec.fits")
f2 = pyfits.open("line_spec.fits")
fig = plt.figure(figsize=(9, 7))
ax = fig.add_subplot(111)
ax.loglog(
    f1["SPECTRUM"].data["ENERGY"],
    f1["SPECTRUM"].data["COUNTS"] + f2["SPECTRUM"].data["COUNTS"],
)
ax.loglog(f2["SPECTRUM"].data["ENERGY"], f2["SPECTRUM"].data["COUNTS"])
ax.set_xlim(0.2, 10)
ax.set_ylim(1, 3.0e4)
ax.set_xlabel("E (keV)")
ax.set_ylabel("counts/bin")
---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
Cell In[14], line 4
      1 import astropy.io.fits as pyfits
      3 f1 = pyfits.open("therm_spec.fits")
----> 4 f2 = pyfits.open("line_spec.fits")
      5 fig = plt.figure(figsize=(9, 7))
      6 ax = fig.add_subplot(111)

File ~/mambaforge/envs/py312/lib/python3.12/site-packages/astropy/io/fits/hdu/hdulist.py:222, in fitsopen(name, mode, memmap, save_backup, cache, lazy_load_hdus, ignore_missing_simple, use_fsspec, fsspec_kwargs, decompress_in_memory, **kwargs)
    219 if not name:
    220     raise ValueError(f"Empty filename: {name!r}")
--> 222 return HDUList.fromfile(
    223     name,
    224     mode,
    225     memmap,
    226     save_backup,
    227     cache,
    228     lazy_load_hdus,
    229     ignore_missing_simple,
    230     use_fsspec=use_fsspec,
    231     fsspec_kwargs=fsspec_kwargs,
    232     decompress_in_memory=decompress_in_memory,
    233     **kwargs,
    234 )

File ~/mambaforge/envs/py312/lib/python3.12/site-packages/astropy/io/fits/hdu/hdulist.py:486, in HDUList.fromfile(cls, fileobj, mode, memmap, save_backup, cache, lazy_load_hdus, ignore_missing_simple, **kwargs)
    467 @classmethod
    468 def fromfile(
    469     cls,
   (...)
    477     **kwargs,
    478 ):
    479     """
    480     Creates an `HDUList` instance from a file-like object.
    481
   (...)
    484     documentation for details of the parameters accepted by this method).
    485     """
--> 486     return cls._readfrom(
    487         fileobj=fileobj,
    488         mode=mode,
    489         memmap=memmap,
    490         save_backup=save_backup,
    491         cache=cache,
    492         ignore_missing_simple=ignore_missing_simple,
    493         lazy_load_hdus=lazy_load_hdus,
    494         **kwargs,
    495     )

File ~/mambaforge/envs/py312/lib/python3.12/site-packages/astropy/io/fits/hdu/hdulist.py:1157, in HDUList._readfrom(cls, fileobj, data, mode, memmap, cache, lazy_load_hdus, ignore_missing_simple, use_fsspec, fsspec_kwargs, decompress_in_memory, **kwargs)
   1154 if fileobj is not None:
   1155     if not isinstance(fileobj, _File):
   1156         # instantiate a FITS file object (ffo)
-> 1157         fileobj = _File(
   1158             fileobj,
   1159             mode=mode,
   1160             memmap=memmap,
   1161             cache=cache,
   1162             use_fsspec=use_fsspec,
   1163             fsspec_kwargs=fsspec_kwargs,
   1164             decompress_in_memory=decompress_in_memory,
   1165         )
   1166     # The Astropy mode is determined by the _File initializer if the
   1167     # supplied mode was None
   1168     mode = fileobj.mode

File ~/mambaforge/envs/py312/lib/python3.12/site-packages/astropy/io/fits/file.py:218, in _File.__init__(self, fileobj, mode, memmap, overwrite, cache, use_fsspec, fsspec_kwargs, decompress_in_memory)
    216     self._open_fileobj(fileobj, mode, overwrite)
    217 elif isinstance(fileobj, (str, bytes)):
--> 218     self._open_filename(fileobj, mode, overwrite)
    219 else:
    220     self._open_filelike(fileobj, mode, overwrite)

File ~/mambaforge/envs/py312/lib/python3.12/site-packages/astropy/io/fits/file.py:641, in _File._open_filename(self, filename, mode, overwrite)
    638 ext = os.path.splitext(self.name)[1]
    640 if not self._try_read_compressed(self.name, magic, mode, ext=ext):
--> 641     self._file = open(self.name, IO_FITS_MODES[mode])
    642     self.close_on_error = True
    644 # Make certain we're back at the beginning of the file
    645 # BZ2File does not support seek when the file is open for writing, but
    646 # when opening a file for write, bz2.BZ2File always truncates anyway.

FileNotFoundError: [Errno 2] No such file or directory: 'line_spec.fits'

Let’s zoom into the region surrounding the line, seeing that it has some broadening due to the random velocities of the dark matter particles:

[15]:
ax.set_xscale("linear")
ax.set_xlim(3, 3.7)
ax.set_ylim(1.0, 3.0e2)
fig
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[15], line 1
----> 1 ax.set_xscale("linear")
      2 ax.set_xlim(3, 3.7)
      3 ax.set_ylim(1.0, 3.0e2)

NameError: name 'ax' is not defined