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) = w₀+w₁a(1-a)

where a is the cosmological expansion factor, and w₀ and w₁ 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:

  <parameter>
    <name>cosmologyMethod</name>
    <value>matter-darkEnergy</value>
  </parameter>
  <parameter>
    <name>darkEnergyEquationOfStateW0</name>
    <value>-0.8</value>
  </parameter>
  <parameter>
    <name>darkEnergyEquationOfStateW1</name>
    <value>0.0</value>
  </parameter>
  <parameter>
    <name>criticalOverdensityMethod</name>
    <value>sphericalTopHatDarkEnergy</value>
  </parameter>
  <parameter>
    <name>virialDensityContrastMethod</name>
    <value>sphericalTopHatDarkEnergy</value>
  </parameter>

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 models¹. For the highest mass galaxies there are small differences as the dark energy equation of state is varied.

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

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:

  <parameter>
    <name>outputRotationCurveData</name>
    <value>true</value>
  </parameter>
  <parameter>
    <name>outputRotationCurveRadii</name>
    <value>
      diskRadius:all:baryonic:loaded:2.2
      diskRadius:all:dark:loaded:2.2
      diskRadius:all:dark:unloaded:2.2
    </value>
  </parameter>

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:componentType:massType:loading?:radius

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/plotRotationCurve.pl

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.

Dust Absorption and Emission Spectra with Grasil

Galacticus can compute the spectrum of light emitted by stars in each model galaxy, and the resulting magnitude in any band. But, this spectrum accounts for absorption by dust in only a very simplistic way, and makes no accounting for re-emission of light from dust grains.

But, by coupling Galacticus with the Grasil code it's possible to compute  spectra of galaxies with a detailed model of dust extinction and re-emission, including the effects of a distribution of dust grain types and sizes, fluctuating grain temperatures, emission from polycyclic aromatic hydrocarbons (PAHs), and radio emission from supernova remnants.

We've recently made it much easier to perform these Galacticus+Grasil calculations in v0.9.2 of Galacticus, and added a short tutorial to the Galacticus documentation explaining how to do this.

Briefly, just add:

    <parameter>
      <name>starFormationHistoriesMethod</name>
      <value>metallicitySplit</value>
    </parameter>

to the Galacticus parameter file. This causes Galacticus to store the entire star formation history of each galaxy, as a function of time and metallicity. (The tutorial explains how to tune the details of how star formation histories are stored.) These star formation histories are used by Grasil to compute the emission arising from stellar populations of different ages.

Then, using Galacticus' standard analysis package, simply load the Galacticus::Grasil module and ask for a property such as "grasilFlux850microns", to get the flux (in Janskys) at any wavelength. An example code snippet to do this looks like:

use Galacticus::HDF5;
use Galacticus::Grasil;
.
.
.
my $galacticus;
$galacticus->{'file' } = "galacticus.hdf5";
$galacticus->{'store'} = 0;
$galacticus->{'tree' } = "all";
&HDF5::Get_Parameters        ($galacticus                         );
&HDF5::Select_Output         ($galacticus, 2.0                    );
&HDF5::Get_Datasets_Available($galacticus                         );
&HDF5::Get_Dataset           ($galacticus,['grasilFlux850microns']);

print $galacticus->{'dataSets'}->{'grasilFlux850microns'}."\n";

When this code is run, Grasil will be run (and automatically downloaded if necessary) on each galaxy. Where possible (i.e. if you have many cores available), multiple copies of Grasil will be run simultaneously. The full SED as a function of wavelength and inclination will be stored to the Galacticus model file for future use (which is convenient as Grasil models typically take tens of seconds to run), and the requested flux computed and returned.

Once SEDs have been computed in this way they can be extracted and examined. A simple plotting script (scripts/plotting/plotGrasilSpectrum.pl) is provided which makes plots like this:

which nicely shows the attenuated stellar emission at short wavelengths, PAH emission in the ~10um region, thermal emission from dust at ~100um and synchrotron emission above 1mm.

Enhanced N-body Merger Tree Handling

In addition to generating its own merger trees, Galacticus can read trees from file and process them to populate each tree with galaxies. Typically, these merger trees have been extracted from an N-body simulation.

The original implementation of processing such trees in Galacticus had some limitations - which were necessary to make the trees fit with certain assumptions made in the way in which Galacticus processes trees as it creates and evolves galaxies. This processing modified the structure of trees and so could affect the resulting physical properties of galaxies.

We've just finished some improvements to v0.9.2 of Galacticus which remove these assumptions and allow N-body trees to be handled without the need to modify them in any way¹.

The original, simplified trees always looked something like this:

Subhalos (orange circles) were supported, but once a halo became a subhalo, it had to stay a subhalo forever. Also, subhalos could never move between branches of the tree (or between different trees). Furthermore, there was no way for a subhalo to exist unless it had, at some point, a non-subhalo (i.e. isolated halo; green circles) progenitor.

This last point is now fixed, allowing tree configurations such as this:

in which a node begins its life as a subhalo, having never had a progenitor².

The second improvement was to allow "branch jumps", in which a subhalo jumps to another branch as shown below:

This can happen if the subhalo becomes unbound from its host (or, perhaps, was never bound), leaves the host and travels to a new host.

The final improvement was to allow subhalos to become isolated halos once again, as in the following tree:

Obviously, these are highly simplified trees - trees extracted from typical N-body simulations often have ~100,000 nodes and much more complicated structure.

It's worth stating briefly Galacticus' philosophy regarding merger trees read from file. Galacticus assumes that the trees are "correct" - that is, whatever the nodes in the tree do, it should be treated as real, physical behavior³. This means that if you put garbage in, you'll get garbage out! Therefore, the input merger trees should be as realistic as possible - typically this means processing them to remove numerical artifacts as much as is possible. This approach seems to be the only reasonable one - Galacticus doesn't know enough about the simulation, halo finding, or merger tree construction algorithms used to create a particular set of trees to intelligently figure out if the trees are physically reasonable. That's left as a task for whoever builds the trees.

Technically, these improvements involved modifying Galacticus to allow arbitrary numbers of "events" to be attached to each node of the merger tree. An event object simply specifies the time at which the event occurs, a paired node involved in the event, and a function to handle the event when it occurs. This allows the evolution of nodes in merger trees to be "locked" until any dependent events have occurred - even if those events occur in other merger trees. Currently, two types of event are handled. The first type are "branch jumps", in which a subhalo jumps from one host halo to another. The second event type is "subhalo promotion" in which a subhalo stops being a subhalo and becomes an isolated halo.

1 The new code has been tested by running on the entire set of trees from the Millennium Simulation which contain a total of 800 million halos. These contain a lot of edge cases and halos doing things you wouldn't expect - Galacticus now handles all of these cases robustly. If you try running Galacticus on trees from some other simulation and find problems - let us know!

2 Such an "initial subhalo" won't do anything interesting at present - it has no way to accrete gas and so can't form a galaxy. That could change in future though, and, in any case, the presence of the subhalo is important as it may interact gravitationally with other subhalos and galaxies.

3 Some modification of the tree is possible within Galacticus. For example, you can tell Galacticus to prune branches of the tree below a mass threshold (if you think that the trees are not reliable below that mass scale for example).

Galacticus v0.9.1 released

After just over one year of development, v0.9.1 of Galacticus is now released! That means that v0.9.1 is now considered "stable" and will receive only bug fixes (not any new features). It also means that v0.9.0 is now considered obsolete and will no longer be supported.

Most of the work on Galacticus v0.9.1 has been to improve robustness, resulting in a much more stable code. However, there are some new features too:

• Several new "simplified" implementations of galaxy physics have been added. These are designed to allow Galacticus to run simplistic models, perhaps for comparison to other semi-analytic models, or to explore the importance of including more detailed physics. They also typically result in shorter run times (at the cost of reduced realism). These new additions include:
• a galactic structure solver in which radii scale linearly with angular momentum;
• an even simpler galactic structure solver which simply fixes radii of galaxies to the virial radius times the spin parameter (multiplied by a user-specified coefficient);
•  simple "power-law" scaling algorithms for cooling rates, feedback in disks and star formation timescales in disks;
• New algorithms have been added to solve the Press-Schechter excursion set problem, allowing for modelling of non-CDM cosmologies (for further details see the previous blog post);
• The reionization epoch (used to cut off accretion into low mass halos) can now be specified via the optical depth to reionization instead of specifying the redshift directly;
• The Run_Galacticus.pl script has some new features:
• Allows for generic analysis scripts to be run on each completed models;
• Models can be run on a Condor cluster.
• Merger trees constructed within Galacticus can now be exported into either Galacticus' native file format or the IRATE format.
• AGN luminosities can now be computed in any filter, using SEDs which depend on the black hole bolometric luminosity, and including the effects of absorption at X-ray wavelengths;
• The black hole/AGN feedback model has been improved to allow "radio-mode" feedback to occur only when the halo is in the "hot-mode" regime, and to allow energy injected into the hot halo to drive mass out of the halo if the input power exceeds the cooling rate;
• The Wetzel & White (2010) satellite merger timescale algorithm is now available;
• It is now possible to construct lightcones:
• Includes a script which grabs all merger trees in a lightcone from the Millennium Database;
• New output and output filtering options allow for only galaxies within the lightcone to be output and for their positions to be translated to lightcone coordinates.

Development now moves to v0.9.2 of Galacticus which has an entirely re-written internal class structure for nodes and galaxies. This makes it very easy to add new components and extend existing ones, and does away with the need to write significant amounts of boilerplate for each new component. Other plans for v0.9.2 include doing away with several of the assumptions which are forced on N-body merger trees to allow handling of fly-by encounters, further optimizations, and improvements to the physics in several areas.

What if dark matter isn't just plain vanilla CDM?

The cold dark matter (CDM) hypothesis has been enormously successful in describing a wide variety of cosmological observations - from the Cosmic Microwave Background to the large scale distribution of galaxies and the properties of galaxies themselves. But, there are a few areas where things aren't so clear cut. These cases are all somewhat ambiguous - because they all involve the complex physics of galaxies and baryons which we don't understand well enough to make really strong statements. But, it's intriguing enough that there's been a lot of interest in alternative types of dark matter recently.

As a result of this, we've just pushed substantial new functionality to Galacticus v0.9.1 (coinciding with a new paper, which appears on arXiv today) which allows it to follow the formation of dark matter structure in non-CDM universes. It's currently set up to work specifically for warm dark matter (WDM) cases, but the underlying algorithms can handle any type of dark matter (for those in the know, I'll add the caveat that it's any type of dark matter whose physics can be described by modifications to the linear theory power spectrum and the barrier for collapse in excursion set theory), due to a neat numerical algorithm figured out by my former student Arya Farahi.

The results look good. Here's a comparison of the dark matter halo mass function with an N-body simulation of warm dark matter by Schneider et al. (2012):

The yellow/red circles were measured from Schneider et al.'s simulation (at low masses the simulation becomes unreliable - shown by the tiny circles - so ignore it in that region). The solid blue line is the result from Galacticus - in excellent agreement with the N-body result and showing the expected cut-off at low masses (well, not that low, Schneider et al. use a relatively light, 0.25keV, WDM particle which isn't consistent with current constraints - but here we're only concerned with comparing the two techniques).

Schneider et al. only include one effect of WDM - the suppression of small scale power in the linear theory power spectrum. A secondary effect, due to the velocity dispersion of the WDM particles, is to make it more difficult for small overdense regions of the Universe to collapse. To make a fair comparison, I ignored that effect in Galacticus also. But, if I add it back in, the result is the green dot-dashed line - it has a very strong effect on the number of low mass halos around the cut-off. (If you want to know what the pink dotted line is all about, check out the preprint!)

The modular nature of Galacticus makes it relatively easy to add in new algorithms such as this - all of its galaxy formation physics will continue to work happily with merger trees built from WDM.

Galacticus Statistics Page

Sometimes procrastination can have interesting results... I was avoiding diving back into a debugging session earlier in the week and so finally finished putting up a simple web page showing a few interesting statistics on the Galacticus project.

The page currently shows three charts (these are annotated timelines powered by Google - easy to use and a nice way to show this information, even if they are Flash based....).

The first chart shows Galacticus revision number (i.e. the latest entry in its Bazaar archive as a function of time). Turns out that we average about 1.1 commits per day!

The second chart shows source lines of code (measured with the excellent SLOCCount, augmented by my own count of embedded directives for the build system). There's been about 20k of new lines added in the past year, brining in a lot of new functionality. I'm hoping to see this curve drop significantly in the next few months as I merge in some very substantial changes to the way in which Galacticus stores and manipulates the objects that it uses to represent galaxies.

Finally, the third chart shows the total instruction fetch count (measured by callgrind) for a simple Galacticus model that gets run automatically every day to benchmark the code. I've only been saving this data for the past few weeks, but already there's a noticeable decline due to a few recent optimizations. This is another line that I'm hoping to see drop significantly - there are a few patches waiting to merge in that should improve performance (particularly when running under OpenMP).