Colossus tutorial: Cosmology Basics¶
Welcome to the Colossus cosmology tutorial.
from __future__ import print_function
import matplotlib.pyplot as plt
import copy
%matplotlib inline
Setting, changing, and adjusting cosmologies¶
Let's begin by importing the Colossus cosmology module:
from colossus.cosmology import cosmology
If this call results in an error, please see the chapter Installation in the documentation. Let's set a cosmology, e.g. the most recent Planck cosmology. This function returns a cosmology object:
cosmo = cosmology.setCosmology('planck18')
Colossus comes with more than 20 hard-coded cosmologies:
for k in cosmology.cosmologies:
print(k)
If we want to retrieve details about a given cosmology, we can just print the cosmology object:
print(cosmo)
In Colossus, the cosmology is set globally. This means that we can obtain the current cosmology object from anywhere in the code:
cosmo = cosmology.getCurrent()
print(cosmo.name)
However, cosmologies are objects, meaning that we can keep multiple cosmologies around. This can be useful when switching back and forth between different cosmologies, for example:
cosmo2 = cosmology.setCosmology('WMAP9')
cosmology.setCurrent(cosmo2)
print(cosmology.getCurrent().name)
cosmology.setCurrent(cosmo)
print(cosmology.getCurrent().name)
Of course, you can also set a user-defined cosmology. In this case, arguments to the constructor of the Cosmology object are passed as a dictionary. Note that there are many more possible arguments than shown here, as described in the documentation.
my_cosmo = {'flat': True, 'H0': 72.0, 'Om0': 0.25, 'Ob0': 0.043, 'sigma8': 0.8, 'ns': 0.97}
cosmo = cosmology.setCosmology('my_cosmo', **my_cosmo)
print(cosmo)
By default, relativistic species (neutrinos and photons) are included, meaning that they contribute to the energy density of the universe and are subtracted from the remaining components. Thus, when setting the flat cosmology above with $\Omega_{\rm m,0}=0.25$ and relspecies=True, the density of dark energy is not 0.75 but rather 0.75 minus the density of relativistic species. In practice, the relativistic contribution is small except at very high redshifts and makes a negligible contribution to most calculations.
It is possible to change parameters of a cosmology after it has been created, but we need to alert Colossus of the change because it needs to clear the cache:
cosmo.Om0 = 0.27
cosmo.checkForChangedCosmology()
print(cosmo)
All the default parameter sets above represent flat $\Lambda CDM$ cosmologies, i.e. model dark energy as a cosmological constant and contain no curvature. To add curvature, the default for flatness must be overwritten, and the dark energy content of the universe must be set explicitly (for flat cosmologies, it is computed from the matter and relativistic contributions). We can either create an updated parameter dictionary and pass it as params, or we can additionally use keyword arguments that overwrite the content of params:
cosmo = cosmology.setCosmology('planck_curv', params = cosmology.cosmologies['planck18'], flat = False, Ode0 = 0.75)
print(cosmo.Ok0)
Multiple models for the dark energy equation of state parameter $w(z)$ are implemented, namely a cosmological constant ($w = -1$), a constant $w$, a linearly varying $w(z) = w_0 + w_a (1 - a)$, and arbitrary user-supplied functions for $w(z)$. To set, for example, a linearly varying EOS, we change the de_model parameter. We can evaluate $w(z)$ using the wz() function:
cosmo = cosmology.setCosmology('planck_w0wa', params = cosmology.cosmologies['planck18'],
de_model = 'w0wa', w0 = -0.8, wa = 0.1)
print(cosmo.wz(0.5))
We can implement more exotic models by supplying an arbitrary function for $w(z)$:
def wz_func(z):
return -1.0 + 0.1 * z
cosmo = cosmology.setCosmology('planck_wz', params = cosmology.cosmologies['planck18'],
de_model = 'user', wz_function = wz_func)
print(cosmo.wz(0.5))
Note that for the cases of the $w_0$-$w_a$ and user-defined cosmologies, the range of future redshifts is restricted to $z=-0.9$ by default because the dark energy component's exponential growth can otherwise lead to overflow errors.
Standard calculations: Densities, distances, times¶
Let's perform some simple calculations such as cosmological densities, distances, and times. As before, we begin by setting a particular cosmology:
cosmo = cosmology.setCosmology('planck18')
Using the object returned by this function (or cosmology.getCurrent()), we can evaluate standard cosmological quantities. For example, let's print the age of the universe today. The cosmology functions always take redshift rather than scale factor as an argument:
print(cosmo.age(0.0))
Note that the returned age is in Gyr. The colossus cosmology functions adhere to the following unit system:
- Length: comoving Mpc/h
- Mass: Msun
- Wavenumber: comoving h/Mpc
- Time: Gyr
- Density: physical Msun * h^2 / kpc^3
Virtually all Colossus functions accept numpy arrays as input and return arrays of the same dimensions:
import numpy as np
z = np.array([0.0, 1.0, 2.0, 3.0])
cosmo.angularDiameterDistance(z)
The colossus cosmology module works to pretty extreme redshifts, namely between $z = -0.995$ and $200$ (though some functions may have stricter limits, for example $z>=0$). For example, let's plot the relative density contribution of the various components in the Planck cosmology:
zp1 = 10**np.arange(-1.0, 2.0, 0.05)
z = zp1 - 1.0
O_b = cosmo.Ob(z)
O_dm = cosmo.Om(z) - O_b
O_de = cosmo.Ode(z)
O_gamma = cosmo.Ogamma(z)
O_nu = cosmo.Onu(z)
plt.figure()
plt.loglog()
plt.xlabel('z+1')
plt.ylabel('Fraction of total density')
plt.plot(zp1, O_dm, '-', label = 'Dark matter')
plt.plot(zp1, O_b, '-', label = 'Baryons')
plt.plot(zp1, O_de, '-', label = 'Dark Energy')
plt.plot(zp1, O_gamma, '-', label = 'Photos')
plt.plot(zp1, O_nu, '-', label = 'Neutrinos')
plt.legend(loc = 4);
There is more to the functions mentioned above: they can also give us their own inverse, e.g. $z(t)$ instead of $t(z)$:
z = np.array([0.0, 1.0, 2.0, 3.0])
t = cosmo.age(z)
z2 = cosmo.age(t, inverse = True)
print(z2)
The slight errors in the recovered redshifts give a hint as to how colossus works internally: the majority of functions rely on interpolation for performance. The interpolating splines can be inverted numerically. If, for some reason, you really need the exact calculation, you can set interpolation = False in the constructor, or at a later time.
Once again using interpolating splines, we can also evaluate the derivatives of many functions, e.g. $dt/dz$ or $d^2t/dz^2$:
cosmo.age(z, derivative = 1)
cosmo.age(z, derivative = 2)
We can even combine the inverse and derivative arguments to evaluate the derivative of the inverse, e.g. $dz/dt$:
t = cosmo.age(z)
cosmo.age(t, inverse = True, derivative = 1)
A word of warning: while the documentation quotes a particular accuracy for each function, the accuracy of the derivatives and inverses is not guaranteed. There are many more functions related to densities, distances, and times. Please consult the documentation of the cosmology module for an exhaustive list.