Power Law

The second example shows how to set up a power-law spectral source, as well as how to add two sets of photons together.

This example will also briefly show how to set up a mock dataset “in memory” using yt. For more details on how to do this, check out the yt docs on in-memory datasets.

Load up the necessary modules:

[1]:
%matplotlib inline
import matplotlib
matplotlib.rc("font", size=18, family="serif")
import yt
import numpy as np
import matplotlib.pyplot as plt
from yt.units import mp, keV, kpc
import pyxsim
/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

The cluster we set up will be a simple isothermal \(\beta\)-model system, with a temperature of 4 keV. We’ll set it up on a uniform grid of 2 Mpc and 256 cells on a side. The parameters of the model are:

[2]:
R = 1000. # radius of cluster in kpc
r_c = 100. # scale radius of cluster in kpc
rho_c = 1.673e-26 # scale density in g/cm^3
beta = 1. # beta parameter
kT = 4. # cluster temperature in keV
nx = 256
ddims = (nx,nx,nx)

and we set up the density and temperature arrays:

[3]:
x, y, z = np.mgrid[-R:R:nx*1j,
                   -R:R:nx*1j,
                   -R:R:nx*1j]
r = np.sqrt(x**2+y**2+z**2)
[4]:
dens = np.zeros(ddims)
dens[r <= R] = rho_c*(1.+(r[r <= R]/r_c)**2)**(-1.5*beta)
dens[r > R] = 0.0
temp = (kT*keV).to_value("K","thermal")*np.ones(ddims)

Next, we will take the density and temperature arrays and put them into a dictionary with their units, where we will also set up velocity fields set to zero. Then, we’ll call the yt function load_uniform_grid to set this up as a full-fledged yt dataset.

[5]:
data = {}
data["density"] = (dens, "g/cm**3")
data["temperature"] = (temp, "K")
data["velocity_x"] = (np.zeros(ddims), "cm/s")
data["velocity_y"] = (np.zeros(ddims), "cm/s")
data["velocity_z"] = (np.zeros(ddims), "cm/s")

bbox = np.array([[-0.5, 0.5], [-0.5, 0.5], [-0.5, 0.5]]) # The bounding box of the domain in code units

L = (2.0*R*kpc).to_value("cm")

# We have to set default_species_fields="ionized" because
# we didn't create any species fields above
ds = yt.load_uniform_grid(data, ddims, L, bbox=bbox,
                          default_species_fields="ionized")
yt : [INFO     ] 2021-09-09 08:22:00,217 Parameters: current_time              = 0.0
yt : [INFO     ] 2021-09-09 08:22:00,218 Parameters: domain_dimensions         = [256 256 256]
yt : [INFO     ] 2021-09-09 08:22:00,219 Parameters: domain_left_edge          = [-0.5 -0.5 -0.5]
yt : [INFO     ] 2021-09-09 08:22:00,220 Parameters: domain_right_edge         = [0.5 0.5 0.5]
yt : [INFO     ] 2021-09-09 08:22:00,220 Parameters: cosmological_simulation   = 0

The next thing we have to do is specify a derived field for the normalization of the power-law emission. This could come from a variety of sources, for example, relativistic cosmic-ray electrons. For simplicity, we’re not going to assume a specific model, except that we will only specify that the source of the power law emission is proportional to the gas mass in each cell:

[6]:
norm = yt.YTQuantity(1.0e-19, "photons/s/keV")
def _power_law_emission(field, data):
    return norm*data["cell_mass"]/(1.0*mp)
ds.add_field(("gas","power_law_emission"), function=_power_law_emission, sampling_type='local', units="photons/s/keV")

where we have normalized the field arbitrarily. Note that the emission field for a power-law model is a bit odd in that it is technically a specific luminosity for the cell. This is done primarily for simplicity in designing the underlying algorithm.

Now, let’s set up a sphere to collect photons from:

[7]:
sp = ds.sphere("c", (0.5, "Mpc"))

And set the parameters for the initial photon sample:

[8]:
A = yt.YTQuantity(500., "cm**2")
exp_time = yt.YTQuantity(1.0e5, "s")
redshift = 0.03

Set up two source models, a thermal model and a power-law model:

[9]:
therm_model = pyxsim.ThermalSourceModel("apec", 0.01, 80.0, 10000, Zmet=0.3)
plaw_model = pyxsim.PowerLawSourceModel(1.0, 0.01, 80.0, "power_law_emission", 1.0)
soxs : [INFO     ] 2021-09-09 08:22:00,858 Using APEC version 3.0.9.
soxs : [INFO     ] 2021-09-09 08:22:00,859 Using apec_v3.0.9_coco.fits for generating spectral lines.
soxs : [INFO     ] 2021-09-09 08:22:00,860 Using apec_v3.0.9_line.fits for generating the continuum.

Now, generate the photons for each source model:

[10]:
ntp, ntc = pyxsim.make_photons("therm_photons", sp, redshift, A, exp_time, therm_model)
npp, npc = pyxsim.make_photons("plaw_photons", sp, redshift, A, exp_time, plaw_model)
pyxsim : [INFO     ] 2021-09-09 08:22:00,983 Cosmology: h = 0.71 100*km/(Mpc*s), omega_matter = 0.27, omega_lambda = 0.73
pyxsim : [INFO     ] 2021-09-09 08:22:00,984 Using emission measure field '('gas', 'emission_measure')'.
pyxsim : [INFO     ] 2021-09-09 08:22:00,985 Using temperature field '('gas', 'temperature')'.
pyxsim : [INFO     ] 2021-09-09 08:27:12,840 Finished generating photons.
pyxsim : [INFO     ] 2021-09-09 08:27:12,841 Number of photons generated: 3436048
pyxsim : [INFO     ] 2021-09-09 08:27:12,842 Number of cells with photons: 316584
pyxsim : [INFO     ] 2021-09-09 08:27:12,906 Cosmology: h = 0.71 100*km/(Mpc*s), omega_matter = 0.27, omega_lambda = 0.73
pyxsim : [INFO     ] 2021-09-09 08:27:14,568 Finished generating photons.
pyxsim : [INFO     ] 2021-09-09 08:27:14,568 Number of photons generated: 111835
pyxsim : [INFO     ] 2021-09-09 08:27:14,569 Number of cells with photons: 93256

Now, we want to project the photons along a line of sight. We’ll specify the "wabs" model for foreground galactic absorption. We’ll create events from the total set of photons as well as the power-law only set, to see the difference between the two.

[11]:
nte = pyxsim.project_photons("therm_photons", "therm_events", "x", (30.0, 45.0), absorb_model="wabs", nH=0.02)
npe = pyxsim.project_photons("plaw_photons", "plaw_events", "x", (30.0, 45.0), absorb_model="wabs", nH=0.02)
pyxsim : [INFO     ] 2021-09-09 08:27:14,584 Foreground galactic absorption: using the wabs model and nH = 0.02.
Projecting photons from cells/particles: 100%|██████████| 316584/316584 [00:00<00:00, 595334.50it/s]
pyxsim : [INFO     ] 2021-09-09 08:27:15,141 Detected 819923 events.
pyxsim : [INFO     ] 2021-09-09 08:27:15,144 Foreground galactic absorption: using the wabs model and nH = 0.02.
Projecting photons from cells/particles: 100%|██████████| 93256/93256 [00:00<00:00, 2768553.55it/s]
pyxsim : [INFO     ] 2021-09-09 08:27:15,196 Detected 68729 events.

Finally, create energy spectra for both sets of events. We won’t bother convolving with instrument responses here, because we just want to see what the spectra look like.

[12]:
et = pyxsim.EventList("therm_events.h5")
ep = pyxsim.EventList("plaw_events.h5")
et.write_spectrum("therm_spec.fits", 0.1, 80.0, 8000, overwrite=True)
ep.write_spectrum("plaw_spec.fits", 0.1, 80.0, 8000, overwrite=True)

To visualize the spectra, we’ll load them up using AstroPy’s FITS I/O and use Matplotlib to plot the spectra:

[13]:
import astropy.io.fits as pyfits
f1 = pyfits.open("therm_spec.fits")
f2 = pyfits.open("plaw_spec.fits")
plt.figure(figsize=(9,7))
plt.loglog(f2["SPECTRUM"].data["ENERGY"], f2["SPECTRUM"].data["COUNTS"], label="Power-Law")
plt.loglog(f2["SPECTRUM"].data["ENERGY"], f1["SPECTRUM"].data["COUNTS"]+f2["SPECTRUM"].data["COUNTS"],
           label="Power-Law+Thermal")
plt.xlim(0.1, 50)
plt.ylim(1, 2.0e4)
plt.legend()
plt.xlabel("E (keV)")
plt.ylabel("counts/bin")
[13]:
Text(0, 0.5, 'counts/bin')
../_images/cookbook_Power_Law_25_1.png

As you can see, the orange line shows a sum of a thermal and power-law spectrum, the latter most prominent at high energies. The blue line shows the power-law spectrum. Both spectra are absorbed at low energies, thanks to the Galactic foreground absorption.