Working with Data and Models From the Literature

A very incomplete set of data from from the literature exist in $ARES/input/litdata. Each file, named using the convention <last name of first author><year>.py, is composed of dictionaries containing the information most useful to ARES (at least when first transcribed). To see a complete listing of options, consult the following list:

import ares

ares.util.lit_options

If any of these papers ring a bell, you can check out the contents in the following way:

litdata = ares.util.read_lit('mirocha2017')  # a self-serving example

or, look directly at the source code, which lives in $ARES/input/litdata. Hopefully the contents of these files are fairly self-explanatory!

We’ll cover a few options below that I’ve used often enough to warrant the development of special routines to interface with the data and/or to plot the results nicely.

The high-z galaxy luminosity function

Measured luminosity functions from the following works are included in ARES:

• Bouwens et al. (2015)
• Finkelstein et al. (2015)
• Parsa et al. (2016)
• van der Burg et al. (2010)

Stellar population models

Currently, ARES can handle both the starburst99 original dataset and the BPASS version 1.0 models (both of which are downloaded automatically). You can access the data via,

s99 = ares.util.read_lit('leitherer1999')
bpass = ares.util.read_lit('eldridge2009')

or, to create more useful objects for handling these data,

s99 = ares.sources.SynthesisModel(source_sed='leitherer1999')
bpass = ares.sources.SynthesisModel(source_sed='eldridge2009')

The spectra for these models are stored in the exact same way to facilitate comparison and uniform use throughout ARES. The most important attributes are wavelengths (or energies or frequencies), times, and data (a 2-D array with shape (wavelengths, times)). So, to compare the spectra for continuous star formation in the steady-state limit (ARES assumes continuous star formation by default), you could do:

import matplotlib.pyplot as pl

pl.loglog(s99.wavelengths, s99.data[:,-1])
pl.loglog(bpass.wavelengths, bpass.data[:,-1])

The most common options for these models are: pop_Z, pop_ssp, pop_binaries, pop_imf, and pop_nebular. See Population Parameters for a description of each of these parameters.

Parametric SEDs for galaxies and quasars

So far, there is only one litdata module in this category: the multi-wavelength AGN template described in Sazonov et al. 2004.

Reproducing Models from ARES Papers

If you’re interested in reproducing a model from a paper exactly, you can either (1) contact me directly for the model of interest, or preferably (someday) download it from my website, or (2) re-compute it yourself. In the latter case, you just need to make sure you supply the required parameters. To facilitate this, I store “parameter files” (just dictionaries) in the litdata framework as well. You can access them like any other dataset from the literature, e.g.,

m17 = ares.util.read_lit('mirocha2017')

A few of the models we focused on most get their own dictionary, for example our reference double power law model for the star-formation efficiency is stored in the dpl variable:

sim = ares.simulations.Global21cm(**m17.dpl)
sim.run()
sim.GlobalSignature()  # voila!

Hopefully this results exactly in the solid black curve from Figure 2 of Mirocha, Furlanetto, & Sun (2017), provided you’re using ARES version 0.2. If it doesn’t, please contact me!

Alternatively, you can use the ParameterBundle framework, which also taps into our collection of data from the literature. To access the set of parameters for the “dpl” model, you simply do:

pars = ares.util.ParameterBundle('mirocha2017:dpl')

This tells ARES to retrieve the dpl variable within the mirocha2017 module. See Parameter Bundles for more on these objects.

Mirocha, Furlanetto, & Sun (2017) (mirocha2017)

This model has a few options: dpl, and the extensions floor and steep, as explored in the paper.

Non-standard pre-requisites:

• High resolution optical depth table for X-ray background. To generate one for yourself, navigate to $ARES/input/optical_depth and open the generate_optical_depth_tables.py file. Between lines 35 and 45 there are a block of parameters that set the resolution of the table. Make sure that helium=1, zi=50, zf=5, and Nz=[1e3]. It should only take a few minutes to generate this table.

The following parameters are uncertain and typically treated as free parameters within the ranges denoted by brackets (not all of which are hard limits):

• pop_Z{0}, \([0.001, 0.04]\)
• pop_Tmin{0}, \([300, \sim \mathrm{few} \times 10^5]\) (pop_Tmin{1} is tied to this value by default).
• pop_fesc{0}, in general can lie in range \([0, 1]\), but for consistency with observational constraints on \(\tau_e\) (from, e.g., Planck), it’s probably best to limit values to \(\lesssim 0.2\).
• pop_fesc_LW{0}, \([0, 1]\)
• pop_rad_yield{1}, \([10^{38}, 10^{42}]\) \(2.6 \times 10^{39}\) by default
• pop_logN{1}, \([18, 23]\), \(-\infty\) by default.

Note

Changes in the metallicity (pop_Z{0}) in general affect the luminosity function (LF). However, by default, the normalization of the star formation efficiency will automatically be adjusted to guarantee that the LF does not change upon changes to pop_Z{0}. Set the pop_calib_lum{0} parameter to None to remove this behavior.

To re-make the right-hand panel of Figure 1 from the paper, you could do something like:

import ares

dpl = ares.util.ParameterBundle('mirocha2017:dpl')

ax = None
for model in ['floor', 'dpl', 'steep']:
    pars = dpl + ares.util.ParameterBundle('mirocha2017:%s' % model)
    pars.update()
    sim = ares.simulations.Global21cm(**pars)
    sim.run()
    ax = sim.GlobalSignature(ax=ax)

For more thorough parameter space explorations, you might want to consider using the ModelGrid (Simple Parameter Study: 2-D Model Grid) or ModelSample (Monte-Carlo Sampling Higher Dimensional Spaces) machinery. If you’d like to do some forecasting or fitting with these models, check out Fitting the Global 21-cm Signal and Fitting Galaxy Luminosity Functions.

Note

Notice that the floor and steep options are defined relative to the dpl model, i.e., they only contain the parameters that are different from the dpl model, which is why we updated the parameter dictionary rather than creating a new one just with the steep or floor parameters.

Furlanetto et al. (2017) (furlanetto2017)

The main options in this model are whether to use momentum-driven or energy-driven feedback, what are accessible separately via, e.g.,

E = ares.util.ParameterBundle('furlanetto2017:energy')
p = ares.util.ParameterBundle('furlanetto2017:momentum')
fshock = ares.util.ParameterBundle('furlanetto2017:fshock')

The only difference is the assumed slope of the star formation efficiency in low-mass halos, which is defined in the parameter 'pq_func_par2[0]', i.e., the third parameter (par2) of the first parameterized quantity ([0]), in addition to a power-law index that describes the rate of redshift evolution, pq_func_par2[1].

To do a quick comparison, you could simply do:

import ares

ls = ['-', '--']
for i, model in enumerate([E, p]):
    pars = model + fshock
    pop = ares.populations.GalaxyPopulation(**pars)
    M = np.logspace(7, 13)
    pl.loglog(M, pop.fstar(z=6, Mh=M), ls=ls[i])

Creating your own

As with parameter bundles, you can write your own litdata modules without modifying the ARES source code. Just create a new .py file and stick it in one of the following places (searched in this order!):

• Your current working directory.
• $HOME/.ares
• $ARES/input/litdata

For example, if I created the following file (junk_lf.py; which you’ll notice resembles the other LF litdata modules) in my current directory:

import numpy as np

redshifts = [4, 5]
wavelength = 1600.
units = {'phi': 1}  # i.e., not abundances not recorded as log10 values

data = {}
data['lf'] = \
{
 4: {
     'M': [-23, -22, -21, -20],
     'phi': list(np.random.rand(4) * 1e-4),
     'err': [tuple(np.random.rand(2) * 1e-7) for i in range(4)]
    },
 5: {
     'M': [-23, -22, -21, -20],
      'phi': list(np.random.rand(4) * 1e-4),
      'err': [tuple(np.random.rand(2) * 1e-7) for i in range(4)],
    }
}

then the built-in plotting routines will automatically find it. For example, you could compare this completely made-up LF with the rest

obslf = ares.analysis.GalaxyPopulation()

ax = obslf.Plot(z=4, sources='junk_lf')
ax = obslf.Plot(z=4, sources='all', round_z=0.2, ax=ax)