Self-Consistent Dust Reddening

In Mirocha, Mason, & Stark (2020), we describe a simple extension to the now-standard galaxy model in ARES that generates ultraviolet colours self-consistently, rather than invoking IRX-\(\beta\) relations to perform a dust correction. Here, we describe how to work with these new models.

Preliminaries

Before getting started, two lookup tables are required that don’t ship by default with ARES (via the remote.py script; see Installation):

• A new halo mass function lookup table that employs the Tinker et al. 2010 results, rather than Sheth-Tormen (used in most earlier work with ARES).

• A lookup table of halo mass assembly histories.

To create these yourself, you’ll need the hmf code, which is installable via pip.

This should only take a few minutes even in serial. First, navigate to $ARES/input/hmf, where you should see a few .hdf5 files and Python scripts. Open the file generate_hmf_tables.py and make the following adjustments to the parameter dictionary:

"hmf_model": 'Tinker10',

"hmf_tmin": 30,
"hmf_tmax": 2e3,
"hmf_dt": 1,

"cosmology_id": 'paul',
"cosmology_name": 'user',
"sigma_8": 0.8159,
'primordial_index': 0.9652,
'omega_m_0': 0.315579,
'omega_b_0': 0.0491,
'hubble_0': 0.6726,
'omega_l_0': 1. - 0.315579,

The new HMF table will use constant 1 Myr time-steps, rather than the default redshift steps, and employ a cosmology designed to remain consistent with another project (led by a collaborator whose name you can probably guess)!

Once you’ve run the generate_hmf_tables.py script, you should have a new file, hmf_Tinker10_user_paul_logM_1400_4-18_t_1971_30-2000.hdf5, sitting inside $ARES/input/hmf. Now, we’re almost done. Simply execute:

python generate_halo_histories.py hmf_Tinker10_user_paul_logM_1400_4-18_t_1971_30-2000.hdf5

The additional resulting file, hgh_Tinker10_user_paul_logM_1400_4-18_t_1971_30-2000_xM_10_0.10.hdf5, will be found automatically in subsequent calculations.

Example

With the required lookup tables now in hand, we can start in the usual way:

import ares
import numpy as np
import matplotlib.pyplot as pl

and read-in the best-fit parameters via

pars = ares.util.ParameterBundle('mirocha2020:univ')

Now, we create a GalaxyPopulation object,

pop = ares.populations.GalaxyPopulation(**pars)

which is an instance of the GalaxyEnsemble class, rather than the GalaxyCohort class, the latter of which is described in example_galaxy_pop. This offers a more general approach, including more complex star formation histories and proper treatment of spectral synthesis.

We can plot the UVLF in the usual way,

# First, some data
obslf = ares.analysis.GalaxyPopulation()
ax = obslf.PlotLF(z=6, round_z=0.2)

# Now, the predicted/calibrated UVLF
mags = np.arange(-25, -5, 0.1)
phi = pop.LuminosityFunction(6, mags)

ax.semilogy(mags, phi)

The main difference between these models and the simpler models (from, e.g., the mirocha2017 parameter bundle) is access to UV colours. The following high-level routine will make a plot like Figure 4 from the paper:

axes = obslf.PlotColors(pop, fig=2)

This will take order \(\simeq 10\) seconds, as modeling UV colours requires synthesizing a reasonbly high resolution (\(\Delta \lambda = 20 \unicode{x212B}\) by default) spectrum for each galaxy in the model, so that there are multiple points within photometric windows. You can adjust the keyword arguments z_uvlf and z_beta to see different redshifts, while sources will control the datasets being plotted.

Note

Computing colors from fits in the Calzetti et al. 1994 windows requires higher resolution, given that these windows are each only \(\Delta \lambda \sim 10 \unicode{x212B}\) wide. Adjust the dlam keyword argument to increase the wavelength-sampling used prior to photometric UV slope estimation.

To access the magnitudes and colours more directly, you can do something like

mags = pop.Magnitude(6., presets='hst')
beta = pop.Beta(6., presets='hst')

fig3, ax3 = pl.subplots(1, 1, num=3)
ax3.scatter(mags, beta, alpha=0.1, color='b', edgecolors='none')

which computes the UV slope \(\beta\) using the Hubble filters appropriate for the input redshift (see Table A1 in the paper).

Note

Other presets include 'jwst', 'jwst-w', 'jwst-m', and 'calzetti1994'.

You can bin points in the \(M_{\text{UV}}-\beta\) plane, if you prefer it, via the return_binned keyword argument,

mags = np.arange(-25, -10, 0.5) # bin centers
beta, beta_s = pop.Beta(6., presets='hst', Mbins=mags, return_binned=True,
        return_scatter=True)

# Plot scatter in each MUV bin as errorbars
ax3.errorbar(mags, beta, yerr=beta_s.T, color='b', marker='s', fmt='o')

Recall that each galaxy in the model actually represents an “average” galaxy in some halo mass bin, i.e., there is not a 1:1 correspondence between galaxies and elements in mags and beta above, which is why we generally weight by halo abundance in each bin. The default mass-bins have \(\Delta \log_{10} M_h = 0.01\) wide, within which ARES inject pop_thin_hist halos and down-weights their abundance accordingly.

Lastly, to extract the raw'' galaxy population properties, like SFR, stellar mass, etc., you can use the ``get_field method, e.g.,

Ms = pop.get_field(6., 'Ms') # stellar mass
Md = pop.get_field(6., 'Md') # dust mass
Sd = pop.get_field(6., 'Sd') # dust surface density
# etc.

To see what’s available, check out

pop.histories.keys()

From these simple commands, most plots and analyses from the paper can be reproduced in relatively short order.

Notes on performance

The compute time for these models is determined largely by three factors:

• The number of halos (really halo bins) we evolve forward in time. The default mirocha2020:univ models use \(\Delta \log_{10} M_h = 0.01\) mass-bins, each with 10 halos (pop_thin_hist=10). For a larger sample, e.g., when lots of scatter is being employed, larger values of pop_thin_hist may be warranted.

• The number of observed redshifts at which \(\beta\) is synthesize, since new spectra must be generated.

• The number of wavelengths used to sample the intrinsic UV spectrum of galaxies. When computing \(\beta\) via photometry (Hubble or JWST), \(\Delta \lambda = 20 \unicode{x212B}\) will generally suffice. However,

For the models in Mirocha, Mason, & Stark (2020), with \(\Delta \log_{10} M_h = 0.01 (\)), pop_thin_hist=10, calibrating at two redshifts for \(\beta\) (4 and 6), with \(\Delta \lambda = 20 \unicode{x212B}\), each model takes \(\sim 10-20\) seconds, depending on your machine.

Note

Generating the UVLF is essentially free compared to computing \(\beta\).

For more information on what’s happening under the hood, e.g., with regards to spectral synthesis and noisy star-formation histories, see More General Star Formation Histories and Spectral Synthesis.