Thursday, June 5, 2014

Stellar Population Synthesis Including Binary Stars

I'm making available a script which will convert the stellar population spectra from the BPASS project into Galacticus HDF5 format, so that they can be used in Galacticus calculations. Just drop the script into your Galacticus directory, run it, and it will download, unpack, and convert the BPASS data. You'll be left with two files (one for populations including only single stars, and another for populations including binaries) in the data/stellarPopulations/ directory that you can then use in Galacticus when computing galaxy luminosities.

Binaries can have a significant effect on galaxy colors (see the papers on the BPASS page for more details). For example, the following plot shows the distribution of u-z colors for galaxies at z=0 when computed using the FSPS stellar population code, and when using the BPASS code (with single and binary populations).

The galaxies with binary stars are offset bluewards by about 0.2 magnitudes. I wonder what effect this has on stellar mass estimates based on broad-band photometry..........?

Thursday, May 22, 2014

Fitting the Stellar Mass Function at z=0

I've just posted this paper to arXiv (full resolution version here, but be warned, it' about 180Mb), which is the first in a series that attempts to carefully explore constraints on Galacticus from observations of the galaxy population.

In this paper, I use an extremely simple galaxy formation model, and a single constraint - the z=0 stellar mass function of galaxies. What sets this study apart from previous, similar studies, is that I attempt to carefully quantify sources of random and systematic error in both the data and the model, and to account for these when constructing the model likelihood function. For example, I compute the full covariance matrix of the stellar mass function, which turns out to look like this (well, technically this is the correlation matrix):

The result is that a good match to the z=0 stellar mass function can be obtained, along with a moderately good "prediction" (which is really an extrapolation) of the evolution of the stellar mass function to z=1. Other statistics, such as the HI mass function of galaxies at z=0, don't work out as well.

Accounting for all sources of random and systematic error significantly broadens the posterior distributions on model parameters - which is important if we want to avoid ruling out viable models. Ignoring these errors also leads to discernable biases in the parameter posteriors.

This is just step one in constructing a data-driven, robustly constrained galaxy formation model. Work on step two is already underway.......

Wednesday, March 12, 2014

Calculating the Cosmic Background Radiation

The development version of Galacticus (v0.9.3) gained a nice new feature today - the ability to compute the evolution of cosmic background radiation self-consistently from the population of galaxies formed in the model. This is in preparation for a summer project (part of Caltech's SURF program) looking at the reionization of the Universe. For now, it's a relatively simple treatment - it include only the contribution from stars (no AGN), doesn't consider factors such as escape fractions or reprocessing of star light by dust - but it works! Here's the results of a test model in which the Universe gets reionized at z=10:

To carry out this calculation we added some useful new technical features to Galacticus. It's always been possible to attach "events" to halos in a Galacticus merger tree - when the halo reaches the event, some action is triggered before the halo is allowed to continue evolving. Now, it's possible to attach events to the entire universe of merger trees in Galacticus (and also to individual trees, but that's another story). Such a "universe event" causes all merger trees to halt at the event time. In the case of this background radiation calculation, universe events are used to halt the evolution of merger trees, compute their net emissivity, evolve the background radiation, and then let the trees evolve further until they reach the next universe event. Galacticus handles all of this transparently - it normally tries to evolve one tree at a time and then destroy it in order to keep its memory footprint small, but if a universe event is added it will process all merger trees simultaneously. This uses more memory, but allows for calculations that depend on the properties of the entire population of galaxies across all merger trees.

Look for more results utilizing this background radiation calculation later this year!

Wednesday, February 19, 2014

Postprocessing Stellar Spectra (For Fun and Profit)

It's been quiet around here recently. Not because development of Galacticus has slowed down - quite the opposite, development is proceeding rapidly and v0.9.4 (to be released around September 2014) should have some big new features. To prove that development is still happening, here's a short tutorial on Galacticus' functionality to postprocess stellar spectra (motivated by the recent addition of the Inoue et al. 2014 model of IGM absorption).

Stellar luminosities are computed by convolving a library of simple stellar populations with the star formation history of each galaxy. Galacticus allows the spectra of those simple stellar populations to be postprocessed (after being read from file or internally generated for example) before they are utilized in the convolution integral. This postprocessing can modify the spectra in arbitrary ways that depend on wavelength, redshift, and age of stellar population. Furthermore, Galacticus allows you to chain together stellar spectra postprocessors into a set to allow multiple postprocessings to be applied. Furthermore again(!), you can define an arbitrary number of sets and apply different sets to different luminosities.

Typical uses of stellar spectra postprocessors include accounting for absorption of galaxy light by the intervening IGM, or capturing only the light from recent star formation (erhaps so that additional dust extinction can be applied to the light of recently formed stars). A full list of the available postprocessors can be found in the Galacticus documentation.

If you don't specify a postprocessing set, the "default" set (consisting of the inoue2014 postprocessor) is applied to each luminosity calculation. To specify other postprocessing sets add the following to your parameter file:

   <value>default recent unabsorbded recentUnabsorbed</value>

where one set must be specified for each luminosity specified in the luminosityFilter parameter. Note that set names can be reused in order to apply the same postprocessor set to multiple luminosities.

The chain of postprocessors to apply for each set is then specified as follows:

   <value>inoue2014 recent</value>

In this case we've constructed three new sets, in addition to the default set (which applies just the inoue2014 postprocessor). The recent set applies both the inoue2014 IGM absorption postprocessor, followed by the recent postprocessor to retain only recently emitted light. The unabsorbed set ignores IGM absorption entirely - it does this by using the identity postprocessor which leaves the spectrum unaffected. Finally, the recentUnabsorbed set applies only the recent filter while ignoring IGM absorption.

In this way it is relatively easy to extract multiple different measures of luminosity from a Galacticus model. For example, you could construct four postprocessor sets, each corresponding to one of the four different IGM absorption models (lycSuppress, madau1995, meiksin2006, and inoue2014 - the latter three of which are shown below for a galaxy at z=3) and apply these to the same luminosity filter to assess how luminosity depends on the IGM model used.

Tuesday, August 20, 2013

Galacticus v0.9.2 Released

We're happy to announce the release of Galacticus v0.9.2. As usual, you can get Galacticus direct from our repo at BitBucket, just use:

hg clone -u v0.9.2

to grab v0.9.2, or you can download a source tarball. v0.9.2 becomes the stable version (supported and will receive bug fixes), v0.9.1 becomes deprecated (no longer supported and will not receive bug fixes) and development moves to v0.9.3.

In addition to the usual complement of bug fixes and minor improvements, this release sees a complete re-write of Galacticus' internal representation of galaxies and their halos. While this has no visible consequences for running Galacticus models, it makes the code much cleaner and greatly simplifies the process of adding new components to galaxies. Components are now specified by a domain specific language - essentially a simple description of the components properties. Galacticus' build system then takes care of generating classes and associated functions for the new component and all of its properties. This saves ~30,000 lines of code which would otherwise have to be written by hand!

In addition to this, there are a few other additions and improvements worth mentioning.

New components

We've added a component which tracks the statistics of recent halo mergers, and another which tracks averages of quantities (currently star formation rates) between successive outputs. Additionally, the merging statistics component now also reports the original depth of each halo in the structure formation hierarchy.


We now make use of nested parallelism to speed up various calculations, including the integration of stellar population SEDs, tabulations of stellar recycling rates and metal yields, and numerical solution to the excursion set problem.

Stellar luminosities

An implementation of the Charlot & Fall (2000) model for dust extinction is now available.

The framework for postprocessing of stellar spectra has been made much more flexible, allowing for different chains of postprocessing filters to be applied to each luminosity computed. This allows, for example, calculation of a luminosity after absorption by the IGM, and the same luminosity after absorption by the IGM but including only light emitted in the past 100Myr for example.

New filters have been added, including a Lyman continuum, HST ACS and WFC3, Herschel SPIRE, Spitzer IRAC, UKIRT J & H, and CCAT SWCam and LWCam filters.


We've updated to interface with the latest versions of Cloudy (13.02) and FSPS (2.4).


The survey output module now supports output of angular diameter distances, angular positions, and distance modulus.

The rotation curve and velocity dispersion of galaxies can be output at arbitrary radii (defined as absolute radii, multiples of galaxy scale lengths, or radii enclosing a specified fraction of mass or light).


We've added support for dark energy cosmologies with equation of state w(a) = w0 + w1 a (1-a).

Internal calculations of the density, enclosed mass, etc. at any point in a galaxy now fully account for adiabatic contraction of the dark matter profile. This allows the  "Cole2000" merger remnant size calculation to be fully consistent with what Cole et al. (2000) actually do.

Critical overdensity for halo collapse and halo virial density contrast can now be computed using the Kitayama & Suto (1996) fitting formula. This is less accurate than a numerical solution, but also faster.

A simple cooling radius model which assumes and isothermal density profile and cooling rate proportional to density squared has been added. This allows for an analytic calculation of cooling radius (so can be useful for fast models).

Tree building

We've added a new merger tree constructor, "fullySpecified", in which the full specification of a merger tree (including all properties of components) is read from an XML document. This is extremely useful for setting up idealized test calculations.

Reading of merger trees from file (usually N-body merger trees) now supports the full complexity of possible tree configurations, including handling of halos which never existed as isolated halos, processing of "forests" of connected merger trees, subhalos which jump between branches of the tree, and subhalos which escape their halo to become isolated halos once again.

Grasil interaction

Grasil is now downloaded as needed, and an analysis module allows calculation of the luminosity under any specified filter to be computed from Grasil spectra.

What to expect in Galacticus v0.9.3

Along with many under-the-hood improvements, and a few minor physics and interface improvements (more realistic hot halo profiles, better handling of subhalo merger times in N-body trees, ability to read merger trees in the IRATE format), v0.9.3 will also introduce tools to construct mock images of galaxy samples from Galacticus, and a complete framework for deriving detailed and accurate constraints on Galacticus model parameters from observed datasets.

Tuesday, May 14, 2013

Support for dark energy cosmologies added to Galacticus

Galacticus has recently (v0.9.2 revision 1809:c3577bea7b39) added support for more general dark energy cosmologies. It's always supported dark energy in the form of a cosmological constant, but since we don't know what dark energy is, it's useful to be able to explore more general equations of state.

Specifically, Galacticus can now compute models in which the equation of state of dark energy has the form:

w(a) = w0+w1a(1-a)

where a is the cosmological expansion factor, and w0 and w1 are adjustable parameters ("darkEnergyEquationOfStateW0" and "darkEnergyEquationOfStateW1" in a Galacticus parameter file).

The effects of dark energy are accounted for in computing the expansion history of the universe, the critical overdensity for collapse of spherical perturbations, and the density contrast of virialized dark matter halos (see the nice review by Percival 2005 for details).

To add dark energy to a model, simply add the following to an input file:




The present day density in dark energy is still set by the usual Omega_DE parameter.

The effects of the dark energy equation of state are small, but noticeable. This figure shows the stellar mass function of galaxies at z=0.07 for three different dark energy models1. For the highest mass galaxies there are small differences as the dark energy equation of state is varied.

1 The model isn't a perfect match to the data (from Li & White 2009). I'm working on that......

Tuesday, February 5, 2013

Rotation Curves

With the latest updates to Galacticus it's now possible to output rotation curves for all galaxies in a model, including the effects of adiabatic contraction on the dark matter profile. This will allow more realistic comparison with measures of the dynamical properties of observed galaxies - such as the Tully-Fisher relation - which are often measured at a specific radius (e.g. 2.2 times the disk radius).

To output rotation curve data, add the following to a Galacticus parameter file:



The first option tells Galacticus to output rotation curves. The second parameter specifies which contributions to the should be counted rotation curve and at what radii. These specifiers break down as follows:


radiusType specifies which radius in the galaxy we're referring to. In the above we specified that radii will be given in units of the disk radius. We could have also used, diskHalfMassRadius, spheroidRadius, spheroidHalfMassRadius, darkMatterScaleRadius, virialRadius, or just radius (which implies radii are given in units of Mpc). The actual radius at which the rotation curve will be output is given by the radius value.

componentType specifies which components of the galaxy+halo should be counted. Currently allowed values are all, disk, spheroid, hotHalo, darkHalo, and blackHole.

massType specifies which types of mass should be counted. Currently allowed values are all, dark, baryonic, galactic, gaseous, stellar, and blackHole.

The loading? option should be either loaded or unloaded, and specifies whether the effect of baryonic loading (i.e. adiabatic contraction) should be included in the calculation of the rotation curve.

Using these specifiers you can request rotation curves at any radius, and get the contribution from any subset of galaxy or halo components. In the above example, we requested the rotation curve at 2.2 times the disk radius, and asked for the contribution from all components. In the first entry we requested the baryonic contribution to the rotation curve, including adiabatic contraction. In the second and third entries we requested the dark matter contribution to the rotation curve, first with and then without the effects of adiabatic contraction.

By outputting at many different radii, you can build up a full rotation curve. For example, here's the result for a roughly L* galaxy at z=0, plotted using scripts/plotting/

It shows quite nicely how the rotation curve stays relatively flat out to large radii, and how this results from the interplay of disk and (adiabatically contracted) dark matter contributions.