Saturday, January 19, 2013

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.

Thursday, January 10, 2013

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 way1.

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 progenitor2.

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 behavior3. 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).