astropy/astropy — pull request #365 by aconley (2012-08-28 23:52:47+00:00)
Title: Implementation of relativistic species (photons, neutrinos) into cosmology classes.

As promised (threatened?) the next step in my master plan is to include photons and neutrinos into the cosmology classes. This is a branch that does just that. It isn't fully tested yet, but I wanted to push it out for discussion because it isn't totally clear what the best way to include photons and neutrinos. The basic problem is just that neutrinos are messy, and the simplest way of doing this (just specify omega_rad) doesn't allow for non-zero neutrino masses at all, even if they are ultra-relativistic.

So the current branch takes a more complicated, but general approach by adding three parameters (Tcmb0, Neff, and m_nu) that allow an exact treatment of massive nuetrinos (under certain assumptions). They all have sensible defaults so that most users need never know they exist, and will just quietly benefit from slightly more accurate calculations.

Before I went any further, I wanted to see if this parameterization made sense -- further discussion below about the alternatives.

The neutrinos treatment is based on Komatsu et al. 2011 ApJS 180, 330 (the WMAP7 cosmology paper),
section 3.3.

* Add parameters Tcmb0, Neff, and m_nu to all classes, with defaults.
* Add nu_relative_density fucntion to compute rho_nu / rho_gamma
* Update all efunc/inv_efunc calculations to include relativistic species
* Because specifying a flat cosmology is slightly tricker when there is
a small non-zero Omega_rel, added a few Flat cosmologies to make things
easier on the user (FlatLambdaCDM, FlatwCDM).
* Updated tests.
* Updated docs.

The new code is not tested in detail -- but the old tests have been
modified to work (since they all assumed Orel = 0.0).

The neutrino treatement is exact but assumes that all neutrino species
have the same mass. For some high masses they aren't properly relativistic,
but in practice they will be so it seems convenient to include them
in Omega_relativistic (note that I don't call it Omega_radiation for this
reason).

It isn't totally clear what the best way to include relativistic species is.
There are a number of possibilities.

* Ignore neutrino masses, assume the normal number of neutrinos, and
work purely with Omega_rad. This is the simplest, but is rather limiting
because it makes totally vanilla assumptions, and lots of people are
interested in massive neutrinos -- which absolutely cannot be implemented
using the pure Omega_rad approach, even if they are very relativistic.
Also, the best fit these days tends to favor N_eff ~ 4, which is rather
interesting, and this parameterization doesn't allow that.
* Allow for neutrinos with masses, but use Omega parameters as much
as possible. This requires: Omega_photon, Omega_nu (the usual
Omega_nu h^2 = \sum m_nu / 94 eV), N_eff (the effective number
of nuetrino species), N_nu (the actual number of neutrino species).
Note that N_eff != N_nu (and they can't be easily calculated from each
other, although they don't differ by much; if you are curious the vanilla model
has N_eff = 3.04 and N_nu=3). Also note that this, the usual way Omega_nu is written,
is very confusing because it virtually -never- actually corresponds to the contribution of neutrinos
to the total energy density-- that is, it looks like an Omega_nu, but actually
is nothing of the sort... Even for low mass neutrinos, Omega_relativistic != Omega_photon + Omega_nu
if you write it this way, which is really confusing. But the advantage is that people tend to expect
Omega_nu to be defined that way, even though doing so forces you to introduce extra parameters
because it's a rather silly definition (you need -both- N_eff and N_nu, whereas more physical
treatments only need N_eff).
* Try to use observable parameters as much as possible; These are Tcmb0, m_nu, and N_eff.
The advantage is that the paramaterization has 1) fewer parameters than the previous one because it
relates better to the actual physics, and 2) it is more directly related to observables. However, the
parameterization may be a little less familiar.

The latter is used in this branch right now because it is simpler, but other people may not feel that way.