The Metagalactic X-ray Background¶
In previous example we’ve calculated the UV background spectrum at a single redshift. Here, we’ll compute the X-ray background evolution over time, and fast!
import ares
import numpy as np
import matplotlib.pyplot as pl
# Parameters for population of stellar mass BHs
params = \
{
'source_type': 'bh',
'Tmin': 1e4,
'spectrum_type': 'mcd',
'source_mass': 10.,
'approx_xrb': False, # Solve the RTE!
'discrete_xrb': True, # Look in $ARES/input/optical_depth for lookup table
'redshift_bins': 400, # Sample redshift interval with this many points
}
# Initialize RadiationBackground instance
rad = ares.solvers.UniformBackground(**params)
# Compute X-ray flux at all redshifts and photon energies
z, E, fluxes = rad.XrayBackground()
Now, we have the cosmic X-ray background flux at all redshifts, z
, and photon
energies, E
. So, the variable fluxes
above is a 2-D array with dimensions
(len(z), len(E))
. If we wanted to plot the background spectrum at a few
redshifts we could do:
for i in range(0, 400, 100):
pl.loglog(E, fluxes[i], label=r'$z=%.3g$' % z[i])
pl.legend()
pl.xlabel(ares.util.labels['E'])
pl.ylabel(ares.util.labels['flux'])
Computing the Heating/Ionization Rate Density¶
With fluxes in hand, we can compute the heating rate density and/or ionization rate density at each redshift straightforwardly:
heat = np.zeros_like(z)
ioniz = np.zeros_like(z)
for i, redshift in enumerate(z):
heat[i] = rad.igm.HeatingRate(redshift, xray_flux=fluxes[i])
ioniz[i] = rad.igm.IonizationRateIGM(redshift, xray_flux=fluxes[i])
Then, plot the results via:
pl.semilogy(z, heat)
or
pl.semilogy(z, ioniz)
These values could be used as input to cosmological simulations, just beware
that by default the values returned by HeatingRate
and IonizationRateIGM
are rates, not rate coefficients, i.e., they have units of \(\mathrm{erg} \ \mathrm{s}^{-1} \ \mathrm{cm}^{-3}\)
for the heating rate density, and \(\mathrm{s}^{-1}\) for the ionization
rate density.
Tabulating the Optical Depth¶
The above example relied on a pre-existing table of the IGM optical depth over
redshift and photon energy, hence the parameter discrete_xrb
, which tells ares
to go looking in $ARES/input/optical_depth
for lookup tables. This technique
was outlined originally in Appendix C of Haardt & Madau (1996).
The shape of the lookup table is defined by the minimum and maximum redshift
(10 and 40 by default), the number of redshift bins used to sample that
interval, redshift_bins
, the minimum and maximum photon energies (0.2 and
30 keV by default), and the number of photon energies (determined iteratively
from the redshift and energy intervals and the value of redshift_bins
).
To make optical depth tables of your own, see $ARES/examples/generate_optical_depth_tables.py
.
By default, ares generates tables assuming the IGM is fully neutral, but that
is not required. See Section 3 of Mirocha (2014)
for more discussion of this technique.
A more complete example can be found in $ARES/tests/test_cxrb_generator.py
.
Alternative Methods¶
The technique outlined above is the fastest way to integrate the cosmological radiative transfer equation (RTE), but it assumes that we can tabulate the optical depth ahead of time. What if instead we wanted to study the radiation background in a decreasingly opaque IGM? Well, we can solve the RTE at several photon energies in turn:
E = np.logspace(2.5, 4.5, 100)
To determine the background intensity at \(z=10\) due to the same BH population as above, we could do something like:
# Function describing evolution of IGM ionized fraction with respect to redshift
# (fully ionized for all time in this case, meaning IGM is optically thin)
xofz = lambda z: 1.0
# Compute flux at z=10 and each observed energy due to emission from
# sources at 10 <= z <= 20.
F = [rad.AngleAveragedFlux(10., nrg, zf=20., xavg=xofz) for nrg in E]
pl.loglog(E, F)
You’ll notice that computing the background intensity is much slower when we do not pre-compute the IGM optical depth.
Let’s compare this to an IGM with evolving ionized fraction:
# Here's a function describing the ionization evolution for a scenario
# in which reionization is halfway done at z=10 and somewhat extended.
xofz2 = lambda z: ares.util.xHII_tanh(z, zr=10., dz=4.)
# Compute fluxes
F2 = [rad.AngleAveragedFlux(10., nrg, zf=20., xavg=xofz2) for nrg in E]
# Plot results
pl.loglog(E, F2)
# Add some nice axes labels
pl.xlabel(ares.util.labels['E'])
pl.ylabel(ares.util.labels['flux'])
Notice how the plot of F2
has been hardened by neutral absorption in the IGM!