Advanced Thermal Emission

In this example, we’ll look at another galaxy cluster, but this time the dataset will have metallicity information for several species in it. In contrast to the thermal emission example, which used a grid-based dataset, the dataset we’ll use here is SPH, taken from the Magneticum suite of simulations. Finally, there are phases of gas in this dataset that will not emit in X-rays, so we also show how to make cuts in phase space and focus only on the X-ray emitting gas. 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]:
%matplotlib inline
import yt
import pyxsim
import soxs
from yt.units import mp
/Users/jzuhone/Source/yt/yt/utilities/logger.py:4: VisibleDeprecationWarning: The configuration file /Users/jzuhone/.config/yt/ytrc is deprecated in favor of /Users/jzuhone/.config/yt/yt.toml. Currently, both are present. Please manually remove the deprecated one to silence this warning.
Deprecated since v4.0.0. This feature will be removed in v4.1.0
  from yt.config import ytcfg

We mentioned above that we wanted to make phase spaces cuts on the gas. Because this is an SPH dataset, the best way to do that is with a “particle filter” in yt.

[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 Magneticum 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")

As mentioned above, a number of elements are tracked in the SPH particles in the dataset (10 to be precise, along with a trace field for the remaining, unspecified metals). We will deal with these elements later, since we want to use the specific mass fractions of these elements to determine their emission line strengths in the mock observation. However, we also need to specify the metallicity for the rest of the (non-hydrogen) elements. The best way to do this is to sum the masses for the metals over every particle and divide by the mass of the particle to get the metallicity for that particle. That will be assumed to be the metallicity for all non-tracked metals in this pyXSIM run. The field in the Magneticum dataset to do this is called ("Gas", "ElevenMetalMasses"), which has a shape of (number of SPH particles, 11). We’ll define this field here and add it to the dataset specifically later:

[4]:
def _metallicity(field, data):
    # We index the array starting with 1 here because the first element is
    # helium (thus not a metal)
    return data["Gas","ElevenMetalMasses"][:,1:].sum(axis=1)/data["Gas","Mass"]

Next, we load the dataset with yt:

[5]:
ds = yt.load('MagneticumCluster/snap_132', long_ids=True, field_spec='magneticum_box2_hr')
yt : [INFO     ] 2021-09-09 11:28:28,337 Calculating time from 9.081e-01 to be 3.929e+17 seconds
yt : [INFO     ] 2021-09-09 11:28:28,339 Assuming length units are in kpc/h (comoving)
yt : [INFO     ] 2021-09-09 11:28:28,429 Parameters: current_time              = 3.9286591383929274e+17 s
yt : [INFO     ] 2021-09-09 11:28:28,430 Parameters: domain_dimensions         = [1 1 1]
yt : [INFO     ] 2021-09-09 11:28:28,431 Parameters: domain_left_edge          = [0. 0. 0.]
yt : [INFO     ] 2021-09-09 11:28:28,431 Parameters: domain_right_edge         = [352000. 352000. 352000.]
yt : [INFO     ] 2021-09-09 11:28:28,433 Parameters: cosmological_simulation   = 1
yt : [INFO     ] 2021-09-09 11:28:28,433 Parameters: current_redshift          = 0.10114286171886899
yt : [INFO     ] 2021-09-09 11:28:28,434 Parameters: omega_lambda              = 0.728
yt : [INFO     ] 2021-09-09 11:28:28,434 Parameters: omega_matter              = 0.272
yt : [INFO     ] 2021-09-09 11:28:28,435 Parameters: omega_radiation           = 0.0
yt : [INFO     ] 2021-09-09 11:28:28,435 Parameters: hubble_constant           = 0.704

and now we add the derived fields and the "hot_gas" particle filter to this dataset. Note that for the derived fields to be picked up by the filter, they must be specified first:

[6]:
ds.add_field(("gas","metallicity"), _metallicity, units="",
             sampling_type="local")
ds.add_particle_filter("hot_gas")
yt : [INFO     ] 2021-09-09 11:28:28,453 Allocating for 3.718e+06 particles
Loading particle index: 100%|██████████| 6/6 [00:00<00:00, 9131.29it/s]
Loading particle index: 100%|██████████| 6/6 [00:00<00:00, 32598.22it/s]
yt : [WARNING  ] 2021-09-09 11:28:34,860 The Derived Field ('hot_gas', 'particle_cylindrical_velocity_z') is deprecated as of yt v4.0.0 and will be removed in yt v4.1.0. Use ('hot_gas', 'particle_velocity_cylindrical_z') instead.
yt : [WARNING  ] 2021-09-09 11:28:35,733 The Derived Field ('hot_gas', 'particle_spherical_velocity_theta') is deprecated as of yt v4.0.0 and will be removed in yt v4.1.0. Use ('hot_gas', 'particle_velocity_spherical_theta') instead.
yt : [WARNING  ] 2021-09-09 11:28:36,700 The Derived Field ('hot_gas', 'particle_spherical_position_theta') is deprecated as of yt v4.0.0 and will be removed in yt v4.1.0. Use ('hot_gas', 'particle_position_spherical_theta') instead.
yt : [WARNING  ] 2021-09-09 11:28:37,331 The Derived Field ('hot_gas', 'particle_cylindrical_velocity_theta') is deprecated as of yt v4.0.0 and will be removed in yt v4.1.0. Use ('hot_gas', 'particle_velocity_cylindrical_theta') instead.
yt : [WARNING  ] 2021-09-09 11:28:38,270 The Derived Field ('hot_gas', 'particle_spherical_position_radius') is deprecated as of yt v4.0.0 and will be removed in yt v4.1.0. Use ('hot_gas', 'particle_position_spherical_radius') instead.
yt : [WARNING  ] 2021-09-09 11:28:38,409 The Derived Field ('hot_gas', 'particle_spherical_velocity_phi') is deprecated as of yt v4.0.0 and will be removed in yt v4.1.0. Use ('hot_gas', 'particle_velocity_spherical_phi') instead.
yt : [WARNING  ] 2021-09-09 11:28:39,705 The Derived Field ('hot_gas', 'particle_spherical_position_phi') is deprecated as of yt v4.0.0 and will be removed in yt v4.1.0. Use ('hot_gas', 'particle_position_spherical_phi') instead.
yt : [WARNING  ] 2021-09-09 11:28:39,786 The Derived Field ('hot_gas', 'particle_spherical_velocity_radius') is deprecated as of yt v4.0.0 and will be removed in yt v4.1.0. Use ('hot_gas', 'particle_velocity_spherical_radius') instead.
[6]:
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).

[7]:
var_elem = {
    elem: ("hot_gas", f"{elem}_fraction")
    for elem in
    ["He", "C", "Ca", "O", "N", "Ne", "Mg", "S", "Si", "Fe"]
}

Now that we have everything we need, we’ll set up the ThermalSourceModel. Because we created a hot gas filter, we will use the "hot_gas" field type for the emission measure, temperature, and metallicity fields.

[8]:
source_model = pyxsim.ThermalSourceModel(
    "apec", 0.1, 10.0, 1000, ("hot_gas","metallicity"),
    temperature_field=("hot_gas","temperature"),
    emission_measure_field=("hot_gas", "emission_measure"),
    var_elem=var_elem
)
soxs : [INFO     ] 2021-09-09 11:28:40,124 Using APEC version 3.0.9.
soxs : [INFO     ] 2021-09-09 11:28:40,127 Using apec_v3.0.9_coco.fits for generating spectral lines.
soxs : [INFO     ] 2021-09-09 11:28:40,127 Using apec_v3.0.9_line.fits for generating the continuum.

As before, we choose big numbers for the collecting area and exposure time, but the redshift should be taken from the cluster itself, since this dataset has a redshift:

[9]:
exp_time = (300., "ks") # exposure time
area = (1000.0, "cm**2") # collecting area
redshift = ds.current_redshift

Next, we’ll create a box object to serve as a source for the photons. The dataset consists of only the galaxy cluster at a specific location, which we use below, and pick a width of 3 Mpc:

[10]:
c = ds.arr([310306.53, 340613.47, 265758.47], "code_length")
width = ds.quan(3.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:

[11]:
n_photons, n_cells = pyxsim.make_photons("snap_132_photons", box, redshift, area, exp_time, source_model)
pyxsim : [INFO     ] 2021-09-09 11:28:40,417 Cosmology: h = 0.704 100*km/(Mpc*s), omega_matter = 0.272, omega_lambda = 0.728
pyxsim : [INFO     ] 2021-09-09 11:28:40,418 Using emission measure field '('hot_gas', 'emission_measure')'.
pyxsim : [INFO     ] 2021-09-09 11:28:40,419 Using temperature field '('hot_gas', 'temperature')'.
pyxsim : [INFO     ] 2021-09-09 11:29:55,186 Finished generating photons.
pyxsim : [INFO     ] 2021-09-09 11:29:55,187 Number of photons generated: 15942426
pyxsim : [INFO     ] 2021-09-09 11:29:55,187 Number of cells with photons: 628704

And now we create events using the project_photons function. Let’s pick an off-axis normal vector, and a north_vector to decide which way is “up.” We’ll use the "wabs" foreground absorption model this time, with a neutral hydrogen column of \(N_H = 10^{20}~{\rm cm}^{-2}\):

[12]:
L = [0.1, 0.2, -0.3] # normal vector
N = [0.0, 1.0, 0.0] # north vector
n_events = pyxsim.project_photons("snap_132_photons", "snap_132_events", L, (45.,30.),
                                  absorb_model="wabs", nH=0.01, north_vector=N)
pyxsim : [INFO     ] 2021-09-09 11:29:55,201 Foreground galactic absorption: using the wabs model and nH = 0.01.
Projecting photons from cells/particles: 100%|██████████| 628704/628704 [00:04<00:00, 131313.54it/s]
pyxsim : [INFO     ] 2021-09-09 11:30:00,087 Detected 10798451 events.

Now that we have a set of “events” on the sky, we can read them in and write them to a SIMPUT file:

[13]:
events = pyxsim.EventList("snap_132_events.h5")
events.write_to_simput("snap_132", overwrite=True)
soxs : [WARNING  ] 2021-09-09 11:30:00,675 Overwriting snap_132_phlist.fits.
soxs : [INFO     ] 2021-09-09 11:30:00,971 Writing source to snap_132_phlist.fits.

We can then use this SIMPUT file 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:

[14]:
soxs.instrument_simulator("snap_132_simput.fits", "evt.fits", (100.0, "ks"), "chandra_acisi_cy0",
                          [45., 30.], overwrite=True)
soxs : [INFO     ] 2021-09-09 11:30:02,403 Making observation of source in evt.fits.
soxs : [INFO     ] 2021-09-09 11:30:02,736 Detecting events from source snap_132
soxs : [INFO     ] 2021-09-09 11:30:02,737 Applying energy-dependent effective area from acisi_aimpt_cy0.arf.
soxs : [INFO     ] 2021-09-09 11:30:03,539 753529 events detected.
soxs : [INFO     ] 2021-09-09 11:30:03,598 Pixeling events.
soxs : [INFO     ] 2021-09-09 11:30:03,807 Scattering events with a multi_image-based PSF.
soxs : [INFO     ] 2021-09-09 11:30:04,443 61063 events were rejected because they do not fall on any CCD.
soxs : [INFO     ] 2021-09-09 11:30:04,503 Scattering energies with RMF acisi_aimpt_cy0.rmf.
soxs : [INFO     ] 2021-09-09 11:30:05,313 Adding background events.
soxs : [INFO     ] 2021-09-09 11:30:05,402 Adding in point-source background.
soxs : [INFO     ] 2021-09-09 11:30:06,103 Detecting events from source ptsrc_bkgnd
soxs : [INFO     ] 2021-09-09 11:30:06,104 Applying energy-dependent effective area from acisi_aimpt_cy0.arf.
soxs : [INFO     ] 2021-09-09 11:30:06,110 13663 events detected.
soxs : [INFO     ] 2021-09-09 11:30:06,114 Pixeling events.
soxs : [INFO     ] 2021-09-09 11:30:06,121 Scattering events with a multi_image-based PSF.
soxs : [INFO     ] 2021-09-09 11:30:06,215 3533 events were rejected because they do not fall on any CCD.
soxs : [INFO     ] 2021-09-09 11:30:06,217 Scattering energies with RMF acisi_aimpt_cy0.rmf.
soxs : [INFO     ] 2021-09-09 11:30:06,511 Generated 10130 photons from the point-source background.
soxs : [INFO     ] 2021-09-09 11:30:06,512 Adding in astrophysical foreground.
soxs : [INFO     ] 2021-09-09 11:30:06,519 Making 4143 events from the astrophysical foreground.
soxs : [INFO     ] 2021-09-09 11:30:06,521 Scattering energies with RMF acisi_aimpt_cy0.rmf.
soxs : [INFO     ] 2021-09-09 11:30:06,632 Adding in instrumental background.
soxs : [INFO     ] 2021-09-09 11:30:06,720 Making 122215 events from the instrumental background.
soxs : [INFO     ] 2021-09-09 11:30:06,767 Writing events to file evt.fits.
soxs : [INFO     ] 2021-09-09 11:30:07,187 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:

[15]:
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:

[16]:
soxs.plot_image("img.fits", stretch='sqrt', cmap='arbre', vmin=0.0, vmax=6.0, width=0.2)
[16]:
(<Figure size 720x720 with 2 Axes>,
 <astropy.visualization.wcsaxes.core.WCSAxes at 0x180b4ed30>)
../_images/cookbook_Advanced_Thermal_Emission_32_1.png