Creating New Source Populations

If you’re interested in using ARES to evolve radiation backgrounds at high-\(z\) but it doesn’t have the source populations you’re interested in, it shouldn’t be too hard to add your own. This section provides a brief tutorial on how to do that.

Before diving in, however, you may want to consider whether the current machinery has all the functionality you need. For example, though the main set of ARES source populations were designed to model galaxies, and so their parameters refer to galaxy properties, fundamentally these populations are all connected to dark matter halo populations. So, if that’s true of your source-population-of-interest as well, you may be able to trick ARES into doing what you want.

The main source populations that you might manipulate are:

  • GalaxyAggregate

    Assumes star formation rate density (SFRD) proportional to rate at which matter collapses into halos. Radiation backgrounds proportional to SFRD. Key parameters include pop_fstar, pop_rad_yield, pop_fesc, which are all effectively population-averaged quantities, hence the “aggregate” designation.

  • GalaxyCohort

    Assumes 1:1 galaxy:halo mapping, and assumes star formation rate (SFR) in galaxies is proportional to the mass accretion rate onto their parent halo. Luminosity of galaxies is proportional to their SFR, so radiation backgrounds are generated by performing an integral over the luminosity function of galaxies. Key parameters include pop_fstar (now a parameterizable function; see The ParameterizedQuantity Framework), pop_rad_yield, etc. Unlike GalaxyAggregate, galaxy properties now vary as a function of mass and redshift. But, because all galaxies formed at the same time with the same mass are identical, they still evolve as cohorts, i.e., galaxies don’t break mass-rank ordering.

If neither of these will accommodate your goals, continue onward! Perhaps someday your new source population will be listed in the list above.

Step One: A Simple Source Population

Let’s create a new source population with the following properties:

  • Objects form at a constant rate (per co-moving volume element).

  • Objects have a flat spectrum.

This isn’t a terribly interesting source population, but it’s one whose behavior we can easily test.

To get started, make a new file, let’s call it DummyPopulation.py, which we’ll eventually add to the aptly-named populations sub-module in ARES, i.e., in $ARES/ares/populations. For now, you can put this file anywhere.

# A few important imports
import numpy as np
from ares.util import ParameterFile
from ares.populations.Halo import HaloPopulation

class DummyPopulation(HaloPopulation):
    def __init__(self, **kwargs):
        HaloPopulation.__init__(self, **kwargs)

A few things to notice off the bat. First, we inherited HaloPopulation (which itself inherits a generic Population object that holds all sorts of meta-data), and called its constructor. This creates a ParameterFile instance and stores it in the attribute pf, which is a dictionary that contains all the parameters relevant to this source population. If you haven’t supplied any parameters, this dictionary will be filled with a bunch of defaults.

Note

Once you move DummyPopulation.py into the ARES source code directory ($ARES/ares/populations), you’ll need to convert to relative imports for ARES. For example, from ares.util import ParameterFile should instead be from ..util import ParameterFile.

Inheritance of the HaloPopulation gives us immediate access to the halo mass function through the halos attribute.

To do some testing, copy-paste the above snippet into a Python terminal and then initialize an instance of your class:

pop = DummyPopulation()

print(pop.halos.tab_z.shape, pop.halos.tab_z.min(), pop.halos.tab_z.max())
print(pop.halos.tab_M.shape, pop.halos.tab_M.min(), pop.halos.tab_M.max())
print(pop.halos.tab_dndm.shape)

So, the halo mass function is stored as a lookup table of shape (1145, 1200), with dimensions corresponding to (redshift, halo mass).

At this stage, your source population isn’t terribly useful – you could have just initialized a HaloPopulation and been done! However, the thing that sets apart source populations from halo populations is the radiation they emit. In order to evolve the cosmic radiation backgrounds generated by sources, the solvers in ARES simply loop over all population objects and access their Emissivity method, which means our DummyPopulation object needs such a method if it’s to be any use to us.

Note

Generally, this emissivity is tied to the SFRD, which (by convention) we simply call SFRD, though ARES will only go looking for the SFRD if requested by the user.

There are a few rules for the Emissivity function.

  • It must take redshift, z, as its first positional argument.

  • It must have three additional keyword arguments, E, Emin, and Emax, all set to None by default.

  • It should return the volume-averaged emissivity of the source population in units of \(\mathrm{erg} \ \mathrm{s}^{-1} \ \mathrm{cm}^{-3}\) if E==None. If E!=None, then the user has requested a specific emissivity, and so the units should be \(\mathrm{erg} \ \mathrm{s}^{-1} \ \mathrm{cm}^{-3} \ \mathrm{eV}^{-1}\). In other words, in the latter case if you integrate the resulting emissivity over an array of photon energies (in eV) from pop_Emin to pop_Emax, you should recover the luminosity density of the population. In both cases, the \(\mathrm{cm}^3\) volume should be a co-moving volume!

So, here’s a start:

def Emissivity(self, z, E=None, Emin=None, Emax=None):
    """
    An informative docstring.
    """

    # Check to see if these sources are on, i.e., is the supplied redshift
    # in the range (pop_zdead, pop_zform).
    on = self.on(z)
    if not np.any(on):
        return z * on

    # Compute the emissivity
    if E is not None:
        # Compute a specific emissivity and return
        raise NotImplemented('do something!')
    else:
        # Compute a band-integrated emissivity and return
        raise NotImplemented('do something else!')

No matter what, the Emissivity method needs to know about the spectral energy distribution (SED) of the source population. The SED is defined by the user via the pop_sed parameter, which can be a variety of things (see Models for Radiation Emitted by Galaxies). This information is parsed by the Population class and stored in an attribute src, which we have access to already because we inherited Population!

The attribute src is an instance of yet another class that represents the SED of an object (these are defined in the sources sub-module of ARES). It’s most important method is Spectrum, which is a function of photon energy (in eV), and returns a normalized version of the SED. See params_sources` for a more detailed description of how this is done.

So, in completing your Emissivity function, you’ll need to use the src.Spectrum function to properly account for the SED of your sources.

Step Two: Integrating with ARES

In the previous section we created a stand-alone (but useless) source population. In order for it to work within ARES, we need to make a few changes in other areas of the code.

First, we must decide how the user is to indicate that this population is of interest, which means setting the pop_sfr_model parameter. Let’s set things up so that passing in pop_sfr_model='dummy' will trigger the creation of our DummyPopulation source. To make sure ARES knows about this, we need to navigate to ares.populations.GalaxyPopulation.py.

Note

Again, since the sources are generally assumed to be galaxies, whose luminosity is mostly from star formation, the main parameter is pop_sfr_model and the class that initializes a “generic” source population is GalaxyPopulation. Apologies if your model does not invoke star formation or galaxies!

Near the bottom of ares.populations.GalaxyPopulation.py, there is a series of if/else statements that are checking the value of pop_sfr_model, and initializing the appropriate class depending on its value. For example, if pop_sfr_model=='sfe-func' we initialize a GalaxyCohort, if ``pop_sfr_model=='fcoll' we initialize a GalaxyAggregate, and so on. Within this if/else block, you need only add

elif model in ['dummy']:
    return DummyPopulation(**kwargs)

You’ll of course need to be sure to add an import statement for DummyPopulation at the top of the file.

Finally, if you haven’t already, move DummyPopulation.py into the ares.populations module, and make sure to convert to relative imports (as noted in previous section).

Step Three: Testing the new population

First, let’s make sure we can initialize an instance of the new source population through ARES:

import ares

pop = ares.populations.GalaxyPopulation(pop_sfr_model='dummy')

and verify that its routines behave as expected.

Now, to verify that the population works within an ARES simulation, let’s compare the results of two calculations: one standard calculation, and the same calculation with this new source population.

Under construction!