Cosmology

This module is an implementation of the standard FLRW cosmology with a number of dark energy models including \(\Lambda CDM\), wCDM, and varying dark energy equations of state. The cosmology object models the contributions from dark matter, baryons, curvature, photons, neutrinos, and dark energy.

Basics

In Colossus, the cosmology is set globally, and all functions respect that global cosmology. Colossus does not set a default cosmology, meaning that the user must set a cosmology before using any cosmological functions or any other functions that rely on the Cosmology module. This documentation contains coding examples of the most common operations. Much more extensive code samples can be found in the Tutorials.

Setting and getting cosmologies

First, we import the cosmology module:

from colossus.cosmology import cosmology

Setting a cosmology is almost always achieved with the setCosmology() function, which can be used in multiple ways:

Set one of the pre-defined cosmologies:
```
cosmology.setCosmology('planck18')
```
Set one of the pre-defined cosmologies, but overwrite certain parameters:
```
cosmology.setCosmology('planck18', print_warnings = False)
```
Use an existing cosmology, overwrite parameters, and rename it. This might be better if we want to overwrite important cosmological parameters:
```
cosmology.setCosmology('myCosmo', params = cosmology.cosmologies['planck18'], Om0 = 0.28)
```
Add a new cosmology to the global list of available cosmologies. This has the advantage that the new cosmology can be set from anywhere in the code. Only the main cosmological parameters are mandatory, all other parameters can be left to their default values:
```
params = {'flat': True, 'H0': 67.2, 'Om0': 0.31, 'Ob0': 0.049, 'sigma8': 0.81, 'ns': 0.95}
cosmology.addCosmology('myCosmo', **params)
cosmo = cosmology.setCosmology('myCosmo')
```

Set a new cosmology without adding it to the global list of available cosmologies:

params = {'flat': True, 'H0': 67.2, 'Om0': 0.31, 'Ob0': 0.049, 'sigma8': 0.81, 'ns': 0.95}
cosmo = cosmology.setCosmology('myCosmo', **params)

Set a self-similar cosmology with a power-law power spectrum of a certain slope, and the default settings set in the powerlaw cosmology:
```
cosmo = cosmology.setCosmology('powerlaw_-2.60')
```

Whichever way a cosmology is set, the current cosmology is stored in a global variable and can be obtained at any time:

cosmo = cosmology.getCurrent()

For more extensive examples, please see the Tutorials.

Changing and switching cosmologies

The current cosmology can also be set to an already existing cosmology object, for example when switching between cosmologies:

cosmo1 = cosmology.setCosmology('WMAP9')
cosmo2 = cosmology.setCosmology('planck18')
cosmology.setCurrent(cosmo1)

The user can change the cosmological parameters of an existing cosmology object at run-time, but MUST call the update function checkForChangedCosmology() directly after the changes. This function ensures that the parameters are consistent (e.g., flatness) and that no outdated cached quantities are used:

cosmo = cosmology.setCosmology('WMAP9')
cosmo.Om0 = 0.31
cosmo.checkForChangedCosmology()

Only user-defined cosmological parameters, that is, parameters that can be passed to the constructor of the cosmology object, can be changed in this way. Changing other internal variables of the class can have unintended consequences! If derived parameters (such as Ok0 or Onu0) are changed, those changes will simply be overwritten when checkForChangedCosmology() is called.

Converting to and from Astropy cosmologies

Colossus can easily interface with the cosmology module of the popular Astropy code. Astropy cosmology objects can be converted to Colossus cosmologies with the fromAstropy() function:

import astropy.cosmology

params = dict(H0 = 70, Om0 = 0.27, Ob0 = 0.0457, Tcmb0 = 2.7255, Neff = 3.046)
sigma8 = 0.82
ns = 0.96

astropy_cosmo = astropy.cosmology.FlatLambdaCDM(**params)
colossus_cosmo = cosmology.fromAstropy(astropy_cosmo, sigma8, ns, cosmo_name = 'my_cosmo')

The cosmo_name parameter is not necessary if a name is set in the Astropy cosmology. The sigma8 and ns parameters must be set by the user because the Astropy cosmology does not contain them (because it does not compute power spectrum-related quantities). The conversion supports the LambdaCDM, FlatLambdaCDM, wCDM, FlatwCDM, w0waCDM, and Flatw0waCDM Astropy cosmology classes. Conversely, to convert a Colossus cosmology to Astropy, we simply use the toAstropy() function:

colossus_cosmo = cosmology.setCosmology('WMAP9')
astropy_cosmo = colossus_cosmo.toAstropy()

Naturally, both conversion functions will fail if astropy.cosmology cannot be imported.

Summary of getter and setter functions

`setCosmology`(cosmo_name[, params])	Set a cosmology.
`addCosmology`(cosmo_name[, params])	Add a set of cosmological parameters to the global list.
`setCurrent`(cosmo)	Set the current global cosmology to a cosmology object.
`getCurrent`()	Get the current global cosmology.
`fromAstropy`(astropy_cosmo, sigma8, ns[, ...])	Convert an Astropy cosmology object to Colossus and set it as current cosmology.

Standard cosmologies

The following sets of cosmological parameters can be chosen using the setCosmology() function:

ID	Paper	Location	Explanation
`planck18-only`	Planck Collab. 2018	Table 2	Best-fit, Planck only (column 5)
`planck18`	Planck Collab. 2018	Table 2	Best-fit with BAO (column 6)
`planck15-only`	Planck Collab. 2015	Table 4	Best-fit, Planck only (column 2)
`planck15`	Planck Collab. 2015	Table 4	Best-fit with ext (column 6)
`planck13-only`	Planck Collab. 2013	Table 2	Best-fit, Planck only
`planck13`	Planck Collab. 2013	Table 5	Best-fit with BAO etc.
`WMAP9-only`	Hinshaw et al. 2013	Table 2	Max. likelihood, WMAP only
`WMAP9-ML`	Hinshaw et al. 2013	Table 2	Max. likelihood, with eCMB, BAO and H0
`WMAP9`	Hinshaw et al. 2013	Table 4	Best-fit, with eCMB, BAO and H0
`WMAP7-only`	Komatsu et al. 2011	Table 1	Max. likelihood, WMAP only
`WMAP7-ML`	Komatsu et al. 2011	Table 1	Max. likelihood, with BAO and H0
`WMAP7`	Komatsu et al. 2011	Table 1	Best-fit, with BAO and H0
`WMAP5-only`	Komatsu et al. 2009	Table 1	Max. likelihood, WMAP only
`WMAP5-ML`	Komatsu et al. 2009	Table 1	Max. likelihood, with BAO and SN
`WMAP5`	Komatsu et al. 2009	Table 1	Best-fit, with BAO and SN
`WMAP3-ML`	Spergel et al. 2007	Table 2	Max.likelihood, WMAP only
`WMAP3`	Spergel et al. 2007	Table 5	Best fit, WMAP only
`WMAP1-ML`	Spergel et al. 2003	Table 1/4	Max.likelihood, WMAP only
`WMAP1`	Spergel et al. 2003	Table 7/4	Best fit, WMAP only
`illustris`	Vogelsberger et al. 2014	–	Cosmology of the Illustris simulation
`bolshoi`	Klypin et al. 2011	–	Cosmology of the Bolshoi simulation
`multidark-planck`	Klypin et al. 2016	Table 1	Cosmology of the Multidark-Planck simulations
`millennium`	Springel et al. 2005	–	Cosmology of the Millennium simulation
`EdS`	–	–	Einstein-de Sitter cosmology
`powerlaw`	–	–	Default settings for power-law cosms.

Those cosmologies that refer to particular simulations (such as bolshoi and millennium) are generally set to ignore relativistic species, i.e. photons and neutrinos, because they are not modeled in the simulations. The EdS cosmology refers to an Einstein-de Sitter model, i.e. a flat cosmology with only dark matter and \(\Omega_{\\rm m} = 1\).

Dark energy and curvature

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 (which is otherwise computed from the matter and relativistic contributions):

cosmo = cosmology.setCosmology('planck_curvature', params = cosmology.cosmologies['planck18'], 
                                                                flat = False, Ode0 = 0.85)

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:

cosmo = cosmology.setCosmology('planck_w0wa', params = cosmology.cosmologies['planck18'], 
                                                                de_model = 'w0wa', w0 = -0.8, wa = 0.1)

We can implement more exotic models by supplying an arbitrary function:

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)

Please note that the redshift range into the future is reduced from \(z = -0.995\) (\(a = 200\)) to \(z = 0.9\) (\(a = 10\)) for the w0wa and user dark energy models to avoid numerical issues.

Power spectrum, variance, correlation function

Many calculations in large-scale structure are based on the the power spectrum of matter fluctuations (Cosmology.matterPowerSpectrum()), the variance within a given radius (Cosmology.sigma()), or the 2-point correlation function (Cosmology.correlationFunction()). The latter two quantities represent integrals over the power spectrum. All models for the power spectrum are housed in a separate power_spectrum module, and can be evaluated with the general powerSpectrum() function. However, in general, only functions internal to the Cosmology class should be needed:

k = 10**np.linspace(-2.0, 3.0, 100)
R = 10**np.linspace(-2.0, 2.0, 100)

P = cosmo.matterPowerSpectrum(k, z = 1.0)
variance = cosmo.sigma(R, z = 1.0)
corr = cosmo.correlationFunction(R, z = 1.0)

The following models for the power spectrum are supported, and are listed in the models dictionary. Their ID can be passed as the model parameter to the power spectrum functions.

Model ID	Function	Reference	Comment
`sugiyama95`	`modelSugiyama95()`	Sugiyama 1995	A semi-analytical fitting function
`eisenstein98`	`modelEisenstein98()`	Eisenstein & Hu 1998	A semi-analytical fitting function (default)
`eisenstein98_zb`	`modelEisenstein98ZeroBaryon()`	Eisenstein & Hu 1998	The zero-baryon version, i.e., no BAO
`camb`	`modelCamb()`	Lewis et al. 2000	The CAMB Boltzmann code (installed separately)

By default, Colossus computes the power spectrum using the eisenstein98 fitting function because it is accurate to better than 5%, fast, and does not need externally installed tools. While this model takes numerous physical effects into account, it overestimates the power on small scales (\(k \\sim 100\\ h/{\\rm Mpc}\)) because it neglects baryon pressure. Moreover, the model always assumes three massless neutrino species, regardless of what is chosen in a given cosmology.

One alternative is to use the camb Boltzmann solver (see the CAMB documentation for details). This code (and its python interface) must be installed externally, e.g., using pip. The first evaluation of the power spectrum for a given cosmology can take a few seconds. However, CAMB takes a rich set of physics into account and can be flexibly tuned using the ps_args parameters. One downside is that the range of wavenumbers is limited, which can lead to an underestimation of the variance on small scales. See the modelCamb() function for details.

All functions that rely on the matter power spectrum (such as variance, correlation function, and functions that depend on them) accept a ps_args dictionary through which the user can set the power spectrum model and parameters. When evaluating the power spectrum (or functions that rely on it) with non-default parameters, it is important to keep those parameters consistent across calculations. In the following example, we evaluate the power spectrum and variance of dark matter only based on a CAMB code, using an increased maximum wavenumber:

ps_args = dict(model = 'camb', kmax = 5E3, ps_type = 'cdm')

P = cosmo.matterPowerSpectrum(k, **ps_args)
variance = cosmo.sigma(R, **ps_args)

The additional arguments can also be passed to all functions in the Large-scale structure and Dark matter halos modules that use the power spectrum. See the Tutorials for more code examples.

Derivatives and inverses

Almost all cosmology functions that are interpolated (e.g., age(), luminosityDistance() or sigma()) can be evaluated as an nth derivative. Please note that some functions are interpolated in log space, resulting in a logarithmic derivative, while others are interpolated and differentiated in linear space. Please see the function documentations below for details.

The derivative functions were not systematically tested for accuracy. Their accuracy will depend on how well the function in question is represented by the interpolating spline approximation. In general, the accuracy of the derivatives will be worse that the error quoted on the function itself, and get worse with the order of the derivative.

Furthermore, the inverse of interpolated functions can be evaluated by passing inverse = True. In this case, for a function y(x), x(y) is returned instead. Those functions raise an Exception if the requested value lies outside the range of the interpolating spline.

The inverse and derivative flags can be combined to give the derivative of the inverse, i.e. dx/dy. Once again, please check the function documentation whether that derivative is in linear or logarithmic units.

Performance optimization and accuracy

This module is optimized for fast performance, particularly in computationally intensive functions such as the correlation function. Almost all quantities are, by default, tabulated, stored in files, and re-loaded when the same cosmology is set again (see the storage module for details). For some rare applications (for example, MCMC chains where functions are evaluated few times, but for a large number of cosmologies), the user can turn this behavior off:

cosmo = cosmology.setCosmology('planck18', interpolation = False, persistence = '')

For more details, please see the documentation of the interpolation and persistence parameters. In order to turn off the interpolation temporarily, the user can simply switch the interpolation parameter off:

cosmo.interpolation = False
Pk = cosmo.matterPowerSpectrum(k)
cosmo.interpolation = True

In this example, the power spectrum is evaluated directly without interpolation. The interpolation is fairly accurate (see specific notes in the function documentation), meaning that it is very rarely necessary to use the exact routines.

Module reference

class cosmology.cosmology.Cosmology(name=None, flat=True, Om0=None, Ode0=None, Ob0=None, H0=None, sigma8=None, ns=None, de_model='lambda', w0=None, wa=None, wz_function=None, relspecies=True, Tcmb0=2.7255, Neff=3.046, power_law=False, power_law_n=0.0, interpolation=True, persistence='rw', print_info=False, print_warnings=True)

A cosmology is set via the parameters passed to the constructor. Any parameter whose default value is None must be set by the user. The easiest way to set these parameters is to use the setCosmology() function with one of the pre-defined sets of cosmological parameters listed above.

In addition, the user can choose between different equations of state for dark energy, including an arbitrary \(w(z)\) function.

A few cosmological parameters are free in principle, but are well constrained and have a sub-dominant impact on the computations. For such parameters, default values are pre-set so that the user does not have to choose them manually. This includes the CMB temperature today (Tcmb0 = 2.7255 K, Fixsen 2009) and the effective number of neutrino species (Neff = 3.046, Planck Collaboration 2018). These values are compatible with the most recent observational measurements and can be changed by the user if necessary.

Parameters

name: str: A name for the cosmology, e.g. WMAP9 or a user-defined name. If a user-defined set of cosmological parameters is used, it is advisable to use a name that does not represent any of the pre-set cosmologies.
flat: bool: If flat, there is no curvature, \(\Omega_{\rm k} = 0\), and the dark energy content of the universe is computed as \(\Omega_{\rm de} = 1 - \Omega_{\rm m} - \Omega_{\gamma} - \Omega_{\nu}\) where \(\Omega_{\rm m}\) is the density of matter (dark matter and baryons) in units of the critical density, \(\Omega_{\gamma}\) is the density of photons, and \(\Omega_{\nu}\) the density of neutrinos. If flat == False, the Ode0 parameter must be passed.
Om0: float: \(\Omega_{\rm m}\), the matter density in units of the critical density at z = 0 (includes all non-relativistic matter, i.e., dark matter and baryons but not neutrinos).
Ode0: float: \(\Omega_{\rm de}\), the dark energy density in units of the critical density at z = 0. This parameter is ignored if flat == True.
Ob0: float: \(\Omega_{\rm b}\), the baryon density in units of the critical density at z = 0.
H0: float: The Hubble constant in km/s/Mpc.
sigma8: float: The normalization of the power spectrum, i.e. the variance when the field is filtered with a top hat filter of radius 8 Mpc/h. See the sigma() function for details on the variance.
ns: float: The tilt of the primordial power spectrum.
de_model: str: An identifier indicating which dark energy equation of state is to be used. The DE equation of state can either be a cosmological constant (de_model = lambda), a constant w (de_model = w0, the w0 parameter must be set), a linear function of the scale factor according to the parameterization of Linder 2003 where \(w(z) = w_0 + w_a (1 - a)\) (de_model = w0wa, the w0 and wa parameters must be set), or a function supplied by the user (de_model = user). In the latter case, the w(z) function must be passed using the wz_function parameter.
w0: float: If de_model == w0, this variable gives the constant dark energy equation of state parameter w. If de_model == w0wa, this variable gives the constant component w (see de_model parameter).
wa: float: If de_model == w0wa, this variable gives the varying component of w, otherwise it is ignored (see de_model parameter).
wz_function: function: If de_model == user, this field must give a function that represents the dark energy equation of state. This function must take z as its only input variable and return w(z).
relspecies: bool: If relspecies == False, all relativistic contributions to the energy density of the universe (such as photons and neutrinos) are ignored. If relspecies == True, their energy densities are computed based on the Tcmb0 and Neff parameters.
Tcmb0: float: The temperature of the CMB at z = 0 in Kelvin.
Neff: float: The effective number of neutrino species.
power_law: bool: Create a self-similar cosmology with a power-law matter power spectrum, \(P(k) = k^{\rm power\_law\_n}\).
power_law_n: float: See power_law.
interpolation: bool: By default, lookup tables are created for certain computationally intensive quantities, cutting down the computation times for future calculations. If interpolation == False, all interpolation is switched off. This can be useful when evaluating quantities for many different cosmologies (where computing the tables takes a prohibitively long time). However, many functions will be much slower if this setting is False, and the derivatives and inverses will not work. Thus, please use interpolation == False only if absolutely necessary.
persistence: str: By default, interpolation tables and other data are stored in a permanent file for each cosmology. This avoids re-computing the tables when the same cosmology is set again. However, if either read or write file access is to be avoided (for example in MCMC chains), the user can set this parameter to any combination of read ('r') and write ('w'), such as 'rw' (read and write, the default), 'r' (read only), 'w' (write only), or '' (no persistence).
print_info: bool: Output information to the console.
print_warnings: bool: Output warnings to the console.

Methods

`Ez`(z)	The Hubble parameter as a function of redshift, in units of \(H_0\).
`Hz`(z)	The Hubble parameter as a function of redshift.
`Ob`(z)	The baryon density of the universe, in units of the critical density.
`Ode`(z)	The dark energy density of the universe, in units of the critical density.
`Ogamma`(z)	The density of photons in the universe, in units of the critical density.
`Ok`(z)	The curvature density of the universe in units of the critical density.
`Om`(z)	The matter density of the universe, in units of the critical density.
`Onu`(z)	The density of neutrinos in the universe, in units of the critical density.
`Or`(z)	The density of relativistic species, in units of the critical density.
`age`(z[, derivative, inverse])	The age of the universe at redshift z.
`angularDiameterDistance`(z[, derivative])	The angular diameter distance to redshift z.
`checkForChangedCosmology`()	Check whether the cosmological parameters have been changed by the user.
`comovingDistance`([z_min, z_max, transverse])	The transverse or line-of-sight comoving distance.
`correlationFunction`(R[, z, derivative, ps_args])	The linear matter-matter correlation function at radius R.
`distanceModulus`(z)	The distance modulus to redshift z in magnitudes.
`filterFunction`(filt, k, R)	The Fourier transform of certain filter functions.
`getName`()	Return the name of this cosmology.
`growthFactor`(z[, derivative, inverse])	The linear growth factor normalized to z = 0, \(D_+(z) / D_+(0)\).
`growthFactorUnnormalized`(z)	The linear growth factor, \(D_+(z)\).
`hubbleTime`(z)	The Hubble time, \(1/H(z)\).
`kpcPerArcsec`(z)	The size (in physical kpc) of an object subtending 1 arcsec at a given redshift.
`lookbackTime`(z[, derivative, inverse])	The lookback time since redshift z.
`luminosityDistance`(z[, derivative, inverse])	The luminosity distance to redshift z.
`matterPowerSpectrum`(k[, z, model, derivative])	The matter power spectrum at a scale k.
`matterPowerSpectrumNorm`(model, **ps_args)	The normalization of the power spectrum to the given sigma8.
`rho_b`(z)	The baryon density of the universe at redshift z.
`rho_c`(z)	The critical density of the universe at redshift z.
`rho_de`(z)	The dark energy density of the universe at redshift z.
`rho_gamma`(z)	The photon density of the universe at redshift z.
`rho_k`(z)	The density of curvature in the universe at redshift z.
`rho_m`(z)	The matter density of the universe at redshift z.
`rho_nu`(z)	The neutrino density of the universe at redshift z.
`rho_r`(z)	The density of relativistic species in the universe at redshift z.
`sigma`(R[, z, j, filt, inverse, derivative, ...])	The rms variance of the linear density field on a scale R, \(\sigma(R)\).
`soundHorizon`()	The sound horizon at recombination.
`toAstropy`()	Create an equivalent Astropy cosmology object.
`wz`(z)	The dark energy equation of state parameter.

getName(): Return the name of this cosmology.

checkForChangedCosmology(): Check whether the cosmological parameters have been changed by the user. If there are changes, all pre-computed quantities (e.g., interpolation tables) are discarded and re-computed if necessary.

toAstropy()

Create an equivalent Astropy cosmology object.

This function throws an error if Astropy cannot be imported. Most standard Colossus cosmologies can be converted, exceptions are self-similar cosmologies with power-law spectra and user-defined dark energy equations of state.

Returns

cosmo_astropy: FLRW: An astropy.cosmology.FLRW class object.

Ez(z)

The Hubble parameter as a function of redshift, in units of \(H_0\).

Parameters

z: array_like: Redshift; can be a number or a numpy array.

Returns

E: array_like: \(H(z) / H_0\); has the same dimensions as z.

See also

Hz: The Hubble parameter as a function of redshift.

Hz(z)

The Hubble parameter as a function of redshift.

Parameters

z: array_like: Redshift; can be a number or a numpy array.

Returns

H: array_like: \(H(z)\) in units of km/s/Mpc; has the same dimensions as z.

See also

Ez: The Hubble parameter as a function of redshift, in units of \(H_0\).

wz(z)

The dark energy equation of state parameter.

The EOS parameter is defined as \(w(z) = P(z) / \rho(z)\). Depending on its chosen functional form (see the de_model parameter to Cosmology()), w(z) can be -1, another constant, a linear function of a, or an arbitrary function chosen by the user.

Parameters

z: array_like: Redshift; can be a number or a numpy array.

Returns

w: array_like: \(w(z)\), has the same dimensions as z.

hubbleTime(z)

The Hubble time, \(1/H(z)\).

Parameters

z: array_like: Redshift; can be a number or a numpy array.

Returns

tH: float: \(1/H\) in units of Gyr; has the same dimensions as z.

See also

lookbackTime: The lookback time since z.
age: The age of the universe at redshift z.

lookbackTime(z, derivative=0, inverse=False)

The lookback time since redshift z.

The lookback time corresponds to the difference between the age of the universe at redshift z and today.

Parameters

z: array_like: Redshift, where \(-0.995 < z < 200\); can be a number or a numpy array.
derivative: int: If greater than 0, evaluate the nth derivative, \(d^nt/dz^n\).
inverse: bool: If True, evaluate \(z(t)\) instead of \(t(z)\). In this case, the z field must contain the time(s) in Gyr.

Returns

t: array_like: The lookback time (or its derivative) since z in units of Gigayears; has the same dimensions as z.

See also

hubbleTime: The Hubble time, \(1/H_0\).
age: The age of the universe at redshift z.

age(z, derivative=0, inverse=False)

The age of the universe at redshift z.

Parameters

z: array_like: Redshift, where \(-0.995 < z < 200\); can be a number or a numpy array.
derivative: int: If greater than 0, evaluate the nth derivative, \(d^nt/dz^n\).
inverse: bool: If True, evaluate \(z(t)\) instead of \(t(z)\). In this case, the z field must contain the time(s) in Gyr.

Returns

t: array_like: The age of the universe (or its derivative) at redshift z in Gigayears; has the same dimensions as z.

See also

hubbleTime: The Hubble time, \(1/H_0\).
lookbackTime: The lookback time since z.

comovingDistance(z_min=0.0, z_max=0.0, transverse=True)

The transverse or line-of-sight comoving distance.

This function returns the comoving distance between two points. Depending on the chosen geometry, the output can have two different meanings. If transverse = False, the line-of-sight distance is returned,

\[d_{\rm com,los}(z) = \frac{c}{H_0} \int_{0}^{z} \frac{1}{E(z)} {\rm d}z .\]

However, if transverse = False, the function returns the comoving distance between two points separated by an angle of one radian at z_max (if z_min is zero). This quantity depends on the spatial curvature of the universe,

\[\begin{split}d_{\rm com,trans}(z) = \left\{ \begin{array}{ll} \frac{c/H_0}{\sqrt{\Omega_{\rm k,0}}} \sinh \left(\frac{\sqrt{\Omega_{\rm k,0}}}{c/H_0} d_{\rm com,los} \right) & \forall \, \Omega_{\rm k,0} > 0 \\ d_{\rm com,los} & \forall \, \Omega_{\rm k,0} = 0 \\ \frac{c/H_0}{\sqrt{-\Omega_{\rm k,0}}} \sin \left(\frac{\sqrt{-\Omega_{\rm k,0}}}{c/H_0} d_{\rm com,los} \right) & \forall \, \Omega_{\rm k,0} < 0 \\ \end{array} \right. \end{split}\]

In Colossus, this distance is referred to as the “transverse comoving distance” (e.g., Hogg 1999), but a number of other terms are used in the literature, e.g., “comoving angular diameter distance” (Dodelson 2003), “comoving coordinate distance” (Mo et al. 2010),or “angular size distance” (Peebles 1993). The latter is not to be confused with the angular diameter distance.

Either z_min or z_max can be a numpy array; in those cases, the same z_min / z_max is applied to all values of the other. If both are numpy arrays, they need to have the same dimensions, and the comoving distance returned corresponds to a series of different z_min and z_max values.

This function does not use interpolation (unlike the other distance functions) because it accepts both z_min and z_max parameters which would necessitate a 2D interpolation. Thus, for fast evaluation, the luminosityDistance() and angularDiameterDistance() functions should be used.

Parameters

zmin: array_like: Redshift; can be a number or a numpy array.
zmax: array_like: Redshift; can be a number or a numpy array.
transverse: bool: Whether to return the transverse of line-of-sight comoving distance. The two are the same in flat cosmologies.

Returns

d: array_like: The comoving distance in Mpc/h; has the same dimensions as zmin and/or zmax.

See also

luminosityDistance: The luminosity distance to redshift z.
angularDiameterDistance: The angular diameter distance to redshift z.

luminosityDistance(z, derivative=0, inverse=False)

The luminosity distance to redshift z.

Parameters

z: array_like: Redshift, where \(-0.995 < z < 200\); can be a number or a numpy array.
derivative: int: If greater than 0, evaluate the nth derivative, \(d^nD/dz^n\).
inverse: bool: If True, evaluate \(z(D)\) instead of \(D(z)\). In this case, the z field must contain the luminosity distance in Mpc/h.

Returns

d: array_like: The luminosity distance (or its derivative) in Mpc/h; has the same dimensions as z.

See also

comovingDistance: The comoving distance between redshift \(z_{\rm min}\) and \(z_{\rm max}\).
angularDiameterDistance: The angular diameter distance to redshift z.

angularDiameterDistance(z, derivative=0)

The angular diameter distance to redshift z.

The angular diameter distance is the line-of-sight distance that allows us to compute the transverse size of an object as if in a non-expanding Universe. In other words, the angular diameter distance is the transverse distance that, at redshift z, corresponds to an angle of one radian. Note that the inverse is not available for this function because it is not strictly increasing or decreasing with redshift, making its inverse multi-valued.

Parameters

z: array_like: Redshift, where \(-0.995 < z < 200\); can be a number or a numpy array.
derivative: int: If greater than 0, evaluate the nth derivative, \(d^nD/dz^n\).

Returns

d: array_like: The angular diameter distance (or its derivative) in Mpc/h; has the same dimensions as z.

See also

comovingDistance: The comoving distance between redshift \(z_{\rm min}\) and \(z_{\rm max}\).
luminosityDistance: The luminosity distance to redshift z.

kpcPerArcsec(z)

The size (in physical kpc) of an object subtending 1 arcsec at a given redshift.

Parameters

z: array_like: Redshift; can be a number or a numpy array.

Returns

size: array_like: The physical size in kpc (not kpc/h); has the same dimensions as z.

distanceModulus(z)

The distance modulus to redshift z in magnitudes.

Parameters

z: array_like: Redshift; can be a number or a numpy array.

Returns

mu: array_like: The distance modulus in magnitudes; has the same dimensions as z.

soundHorizon()

The sound horizon at recombination.

This function returns the sound horizon in Mpc/h, according to Eisenstein & Hu 1998, Equation 26. This fitting function is accurate to 2% where \(\Omega_{\rm b} h^2 > 0.0125\) and \(0.025 < \Omega_{\rm m} h^2 < 0.5\).

Returns

s: float: The sound horizon at recombination in Mpc/h.

rho_c(z)

The critical density of the universe at redshift z.

Parameters

z: array_like: Redshift; can be a number or a numpy array.

Returns

rho_critical: array_like: The critical density in units of physical \(M_{\odot} h^2 / {\rm kpc}^3\); has the same dimensions as z.

rho_m(z)

The matter density of the universe at redshift z.

Parameters

z: array_like: Redshift; can be a number or a numpy array.

Returns

rho_matter: array_like: The matter density in units of physical \(M_{\odot} h^2 / {\rm kpc}^3\); has the same dimensions as z.

See also

Om: The matter density of the universe, in units of the critical density.

rho_b(z)

The baryon density of the universe at redshift z.

Parameters

z: array_like: Redshift; can be a number or a numpy array.

Returns

rho_baryon: array_like: The baryon density in units of physical \(M_{\odot} h^2 / {\rm kpc}^3\); has the same dimensions as z.

See also

Ob: The baryon density of the universe, in units of the critical density.

rho_de(z)

The dark energy density of the universe at redshift z.

Parameters

z: array_like: Redshift; can be a number or a numpy array.

Returns

rho_de: float: The dark energy density in units of physical \(M_{\odot} h^2 / {\rm kpc}^3\); has the same dimensions as z.

See also

Ode: The dark energy density of the universe, in units of the critical density.

rho_gamma(z)

The photon density of the universe at redshift z.

If relspecies == False, this function returns 0.

Parameters

z: array_like: Redshift; can be a number or a numpy array.

Returns

rho_gamma: array_like: The photon density in units of physical \(M_{\odot} h^2 / {\rm kpc}^3\); has the same dimensions as z.

See also

Ogamma: The density of photons in the universe, in units of the critical density.

rho_nu(z)

The neutrino density of the universe at redshift z.

If relspecies == False, this function returns 0.

Parameters

z: array_like: Redshift; can be a number or a numpy array.

Returns

rho_nu: array_like: The neutrino density in units of physical \(M_{\odot} h^2 / {\rm kpc}^3\); has the same dimensions as z.

See also

Onu: The density of neutrinos in the universe, in units of the critical density.

rho_r(z)

The density of relativistic species in the universe at redshift z.

This density is the sum of the photon and neutrino densities. If relspecies == False, this function returns 0.

Parameters

z: array_like: Redshift; can be a number or a numpy array.

Returns

rho_relativistic: array_like: The density of relativistic species in units of physical \(M_{\odot} h^2 / {\rm kpc}^3\); has the same dimensions as z.

See also

Or: The density of relativistic species in the universe, in units of the critical density.

rho_k(z)

The density of curvature in the universe at redshift z.

While the meaning of the “density” of curvature is not easy to interpret, this term enters the Friedmann equations (for non-flat universes) just like all the other densities.

Parameters

z: array_like: Redshift; can be a number or a numpy array.

Returns

rho_curvature: array_like: The equivalent density of curvature in units of physical \(M_{\odot} h^2 / {\rm kpc}^3\); has the same dimensions as z.

Om(z)

The matter density of the universe, in units of the critical density.

Parameters

z: array_like: Redshift; can be a number or a numpy array.

Returns

Omega_matter: array_like: Has the same dimensions as z.

See also

rho_m: The matter density of the universe at redshift z.

Ob(z)

The baryon density of the universe, in units of the critical density.

Parameters

z: array_like: Redshift; can be a number or a numpy array.

Returns

Omega_baryon: array_like: Has the same dimensions as z.

See also

rho_b: The baryon density of the universe at redshift z.

Ode(z)

The dark energy density of the universe, in units of the critical density.

Parameters

z: array_like: Redshift; can be a number or a numpy array.

Returns

Omega_de: array_like: Has the same dimensions as z.

See also

rho_de: The dark energy density of the universe at redshift z.

Ogamma(z)

The density of photons in the universe, in units of the critical density.

Parameters

z: array_like: Redshift; can be a number or a numpy array.

Returns

Omega_gamma: array_like: Has the same dimensions as z.

See also

rho_gamma: The photon density of the universe at redshift z.

Onu(z)

The density of neutrinos in the universe, in units of the critical density.

Parameters

z: array_like: Redshift; can be a number or a numpy array.

Returns

Omega_nu: array_like: Has the same dimensions as z.

See also

rho_nu: The neutrino density of the universe at redshift z.

Or(z)

The density of relativistic species, in units of the critical density.

This function returns the sum of the densities of photons and neutrinos.

Parameters

z: array_like: Redshift; can be a number or a numpy array.

Returns

Omega_relativistic: array_like: Has the same dimensions as z.

See also

rho_r: The density of relativistic species in the universe at redshift z.

Ok(z)

The curvature density of the universe in units of the critical density.

In a flat universe, \(\Omega_{\rm k} = 0\).

Parameters

z: array_like: Redshift; can be a number or a numpy array.

Returns

Omega_curvature: array_like: Has the same dimensions as z.

growthFactorUnnormalized(z)

The linear growth factor, \(D_+(z)\).

The growth factor describes the linear evolution of over- and underdensities in the dark matter density field. There are three regimes:

In the matter-radiation regime, we use an approximate analytical formula (Equation 5 in Gnedin et al. 2011. If relativistic species are ignored, \(D_+(z) \propto a\).
In the matter-dominated regime, \(D_+(z) \propto a\).
In the matter-dark energy regime, we evaluate \(D_+(z)\) through integration as defined in Eisenstein & Hu 1999, Equation 8 (see also Heath 1977) for LCDM cosmologies. For cosmologies where \(w(z) \neq -1\), this expression is not valid and we instead solve the ordinary differential equation for the evolution of the growth factor (Equation 11 in Linder & Jenkins 2003).

At the transition between the integral and analytic approximation regimes, the two expressions do not quite match up, with differences of the order <1E-3. in order to avoid a discontinuity, we introduce a transition regime where the two quantities are linearly interpolated.

The normalization is such that the growth factor approaches \(D_+(a) = a\) in the matter-dominated regime. There are other normalizations of the growth factor (e.g., Percival 2005, Equation 15), but since we almost always care about the growth factor normalized to z = 0, the normalization does not matter too much (see the growthFactor() function).

Parameters

z: array_like: Redshift, where \(-0.995 < z\); the high end of z is only limited by the validity of the analytical approximation mentioned above. Can be a number or a numpy array.

Returns

D: array_like: The linear growth factor; has the same dimensions as z.

Warning

This function directly evaluates the growth factor by integration or analytical approximation. In most cases, the growthFactor() function should be used since it interpolates and is thus much faster.

See also

growthFactor: The linear growth factor normalized to z = 0, \(D_+(z) / D_+(0)\).

growthFactor(z, derivative=0, inverse=False)

The linear growth factor normalized to z = 0, \(D_+(z) / D_+(0)\).

The growth factor describes the linear evolution of over- and underdensities in the dark matter density field. This function is sped up through interpolation which barely degrades its accuracy, but if you wish to evaluate the exact integral or compute the growth factor for very high redshifts (z > 200), please use the growthFactorUnnormalized() function.

Parameters

z: array_like: Redshift, where \(-0.995 < z < 200\); can be a number or a numpy array.
derivative: int: If greater than 0, evaluate the nth derivative, \(d^nD_+/dz^n\).
inverse: bool: If True, evaluate \(z(D_+)\) instead of \(D_+(z)\). In this case, the z field must contain the normalized growth factor.

Returns

D: array_like: The linear growth factor (or its derivative); has the same dimensions as z.

See also

growthFactorUnnormalized: The linear growth factor, \(D_+(z)\).

matterPowerSpectrum(k, z=0.0, model='eisenstein98', derivative=False, **ps_args)

The matter power spectrum at a scale k.

This general function can return power spectra computed by a variety of models, tables, or Boltzmann codes, as described in the powerSpectrum() function. By default, an interpolation table is generated once and future calls to this function will return interpolated values.

By default, the power spectrum is computed using a model for the transfer function (see the powerSpectrum() function). The default Eisenstein & Hu 1998 approximation is accurate to about 5%, and the interpolation introduces errors significantly smaller than that.

Alternatively, the user can supply the path to a file with a tabulated power spectrum by passing a path parameter in the keyword arguments. The file must contain two columns, \(\log_{10}(k)\) and \(\log_{10}(P)\) where k and P(k) are in the same units as in this function. This table is interpolated with a third-order spline. Note that the tabulated spectrum is normalized to the value if \(\sigma_8\) set in the cosmology. If a power spectrum is to be read from a file, the persistence parameter must allow for reading (though not necessarily writing) of files.

Parameters

k: array_like: The wavenumber k (in comoving h/Mpc), where \(10^{-20} < k < 10^{20}\); can be a number or a numpy array. If a user-supplied table is used, the limits of that table apply.
z: float: The redshift at which the power spectrum is evaluated, zero by default. If non-zero, the power spectrum is scaled with the linear growth factor, \(P(k, z) = P(k, 0) D_{+}^2(z)\).
model: str: A model for the power spectrum (see the cosmology.power_spectrum module). If a tabulated power spectrum is used (see path parameter), this name must still be passed. Internally, the power spectrum is saved using this name, so the name must not overlap with any other models.
derivative: bool: If False, return P(k). If True, return \(d \log(P) / d \log(k)\).
ps_args: kwargs: Arguments passed to the powerSpectrum() function. One parameter can be path, which can contain a tabulated power spectrum, where the two columns are \(\log_{10}(k)\) (in comoving h/Mpc) and \(\log_{10}(P)\) (in \(({\rm Mpc}/h)^3\)).

Returns

P: array_like: The matter power spectrum; has the same dimensions as k and units of \(({\rm Mpc}/h)^3\). Alternatively, the dimensionless logarithmic derivative if derivative == True.

Warning

If a user-supplied power spectrum table is used, integrals over the power spectrum such as the variance and correlation function are integrated only within the limits of the given power spectrum. By default, Boltzmann codes return a relatively small range in wavenumber. Please increase this range if necessary, and check that the computed quantities are converged.

matterPowerSpectrumNorm(model, **ps_args)

The normalization of the power spectrum to the given sigma8.

Internally, Colossus normalizes all power spectra such that they integrate to the value of sigma8 set in the cosmology object, regardless of whether the power spectrum was generated using a fitting function, table, or Boltzmann code. In some cases, these normalizations may be useful, for example when comparing the original amplitude of the power spectra of components such as dark matter and baryons.

Parameters

model: str: A model for the power spectrum (see the cosmology.power_spectrum module).
ps_args: kwargs: Arguments passed to the powerSpectrum() function.

Returns

norm: float: The number that, when multiplied with the power spectrum as output by the powerSpectrum() function, normalizes the spectrum to the sigma8 set in the cosmology. Note that this number can take on arbitrary values that depend on the power spectrum models and/or user-defined parameters.

filterFunction(filt, k, R)

The Fourier transform of certain filter functions.

The main use of the filter function is in computing the variance, please see the documentation of the sigma() function for details. This function is dimensionless, the input units are k in comoving h/Mpc and R in comoving Mpc/h. Available filters are tophat,

\[\tilde{W}_{\rm tophat} = \frac{3}{(kR)^3} \left[ \sin(kR) - kR \times \cos(kR) \right] \,,\]

a gaussian filter,

\[\tilde{W}_{\rm gaussian} = \exp \left[ \frac{-(kR)^2}{2} \right] \,,\]

and a sharp-k filter,

\[\tilde{W}_{\rm sharp-k} = \Theta(1 - kR) \,, \]

where \(\Theta\) is the Heaviside step function.

Parameters

filt: str: Either tophat (a top-hat filter in real space), sharp-k (a top-hat filter in Fourier space), or gaussian (a Gaussian in both real and Fourier space).
k: float: A wavenumber k (in comoving h/Mpc).
R: float: A radius R (in comoving Mpc/h).

Returns

filter: float: The value of the filter function.

sigma(R, z=0.0, j=0, filt='tophat', inverse=False, derivative=False, kmin=None, kmax=None, ps_args={'model': 'eisenstein98', 'path': None})

The rms variance of the linear density field on a scale R, \(\sigma(R)\).

The variance and its higher moments are defined as the integral

\[\sigma^2(R,z) = \frac{1}{2 \pi^2} \int_0^{\infty} k^2 k^{2j} P(k,z) |\tilde{W}(kR)|^2 dk\]

where \(\tilde{W}(kR)\) is the Fourier transform of the filterFunction(), and \(P(k,z) = D_+^2(z)P(k,0)\) is the matterPowerSpectrum(). See the documentation of filterFunction() for possible filters.

By default, the power spectrum is computed using the transfer function approximation of Eisenstein & Hu 1998 (see the cosmology.power_spectrum module). With this approximation, the variance is accurate to about 2% or better (see the Colossus code paper for details). Using a tabulated power spectrum can make this computation more accurate, but please note that the limits of the corresponding table are used for the integration.

Higher moments of the variance (such as \(\sigma_1\), \(\sigma_2\) etc) can be computed by setting j > 0 (see Bardeen et al. 1986). Furthermore, the logarithmic derivative \(d \log(\sigma) / d \log(R)\) can be evaluated by setting derivative == True.

Parameters

R: array_like: The radius of the filter in comoving Mpc/h, where \(10^{-12} < R < 10^3\); can be a number or a numpy array.
z: float: Redshift; for z > 0, \(\sigma(R)\) is multiplied by the linear growth factor.
j: integer: The order of the integral. j = 0 corresponds to the variance, j = 1 to the same integral with an extra \(k^2\) term etc; see Bardeen et al. 1986 for mathematical details.
filt: str: Either tophat, sharp-k or gaussian (see filterFunction()). Higher moments (j > 0) can only be computed for the gaussian filter.
inverse: bool: If True, compute \(R(\sigma)\) rather than \(\sigma(R)\).
derivative: bool: If True, return the logarithmic derivative, \(d \log(\sigma) / d \log(R)\), or its inverse, \(d \log(R) / d \log(\sigma)\) if inverse == True.
kmin: float: The lower limit of the variance integral in k-space. If None, the limit is determined automatically (it should be zero in principle, but will be set to a very small number depending on the type of power spectrum). Setting kmin can be useful when considering finite simulation volumes where the largest scales (smallest k-modes) are not taken into account. If kmin is not None and a top-hat filter is used, the variance will oscillate at certain radii. Thus, the number of points in the interpolation table is automatically increased, but the solution is nevertheless unreliable close to the cutoff scale.
kmax: float: The upper limit of the variance integral in k-space. If None, the limit is determined automatically (it should be infinity in principle, but will be set to a large number where the power spectrum has fallen off sufficiently that the integral is converged). If kmax is not None and a top-hat filter is used, the variance will oscillate at certain radii. Thus, the number of points in the interpolation table is automatically increased, but the solution is nevertheless unreliable close to the cutoff scale.
ps_args: dict: Arguments passed to the matterPowerSpectrum() function.

Returns

sigma: array_like: The rms variance; has the same dimensions as R. If inverse and/or derivative are True, the inverse, derivative, or derivative of the inverse are returned. If j > 0, those refer to higher moments.

See also

matterPowerSpectrum: The matter power spectrum at a scale k.

correlationFunction(R, z=0.0, derivative=False, ps_args={'model': 'eisenstein98', 'path': None})

The linear matter-matter correlation function at radius R.

The linear correlation function is defined as

\[\xi(R,z) = \frac{1}{2 \pi^2} \int_0^\infty k^2 P(k,z) \frac{\sin(kR)}{kR} dk\]

where P(k) is the matterPowerSpectrum(). By default, the power spectrum is computed using the transfer function approximation of Eisenstein & Hu 1998 (see the cosmology.power_spectrum module). With this approximation, the correlation function is accurate to ~5% over the range \(10^{-2} < R < 200\) (see the Colossus code paper for details). Using a tabulated power spectrum can make this computation more accurate, but please note that the limits of the corresponding table are used for the integration.

Parameters

R: array_like: The radius in comoving Mpc/h; can be a number or a numpy array.
z: float: Redshift; if non-zero, the correlation function is scaled with the linear growth factor, \(\xi(R, z) = \xi(R, 0) D_{+}^2(z)\).
derivative: bool: If derivative == True, the linear derivative \(d \xi / d R\) is returned.
ps_args: dict: Arguments passed to the matterPowerSpectrum() function.

Returns

xi: array_like: The correlation function, or its derivative; has the same dimensions as R.

See also

matterPowerSpectrum: The matter power spectrum at a scale k.

cosmology.cosmology.setCosmology(cosmo_name, params=None, **kwargs)

Set a cosmology.

This function provides a convenient way to create a cosmology object without setting the parameters of the Cosmology class manually. See the Basic Usage section for examples. Whichever way the cosmology is set, the global variable is updated so that the getCurrent() function returns the set cosmology.

Parameters

cosmo_name: str: The name of the cosmology. Can be the name of a pre-set cosmology, or another name in which case the params dictionary needs to be provided.
params: dictionary: The parameters of the constructor of the Cosmology class. Not necessary if cosmo_name is the name of a pre-set cosmology. Alternatively, keyword arguments can be set directly.
kwargs: kwargs: A set of keyword arguments that are passed to the constructor of the cosmology class. Equivalent to setting a dictionary in params, but these arguments overwrite params if a parameter is present in both.

Returns

cosmo: Cosmology: The created cosmology object.

cosmology.cosmology.addCosmology(cosmo_name, params={}, **kwargs)

Add a set of cosmological parameters to the global list.

After this function is executed, the new cosmology can be set using setCosmology() from anywhere in the code.

Parameters

cosmo_name: str: The name of the cosmology.
params: dictionary: A set of parameters for the constructor of the Cosmology class.
kwargs: kwargs: A set of keyword arguments that are passed to the constructor of the cosmology class. Equivalent to setting a dictionary in params, but these arguments overwrite params if a parameter is present in both.

cosmology.cosmology.setCurrent(cosmo)

Set the current global cosmology to a cosmology object.

Unlike setCosmology(), this function does not create a new cosmology object, but allows the user to set a cosmology object to be the current cosmology. This can be useful when switching between cosmologies, since many routines use the getCurrent() routine to obtain the current cosmology.

Parameters

cosmo: Cosmology: The cosmology object to be set as the global current cosmology.

cosmology.cosmology.getCurrent()

Get the current global cosmology.

This function should be used whenever access to the cosmology is needed. By using the globally set cosmology, there is no need to pass cosmology objects around the code. If no cosmology is set, this function raises an Exception that reminds the user to set a cosmology.

Returns

cosmo: Cosmology: The current globally set cosmology.

cosmology.cosmology.fromAstropy(astropy_cosmo, sigma8, ns, cosmo_name=None, **kwargs)

Convert an Astropy cosmology object to Colossus and set it as current cosmology.

This function throws an error if Astropy cannot be imported. The user needs to supply sigma8 and ns because Astropy does not include power spectrum calculations in its cosmology module. There are a few dark energy model implemented in Astropy that are not implemented in Colossus, namely “wpwaCDM” and “w0wzCDM”. The corresponding classes cannot be converted with this function.

Parameters

astropy_cosmo: FLRW: Astropy cosmology object of type FLRW or its derivative classes.
sigma8: float: The normalization of the power spectrum, i.e. the variance when the field is filtered with a top hat filter of radius 8 Mpc/h. This parameter is not part of an Astropy cosmology but must be set in Colossus.
ns: float: The tilt of the primordial power spectrum. This parameter is not part of an Astropy cosmology but must be set in Colossus.
cosmo_name: str: The name of the cosmology. Can be None if a name is supplied in the Astropy cosmology, which is optional in Astropy.
kwargs: dictionary: Additional parameters that are passed to the Colossus cosmology.

Returns

cosmo: Cosmology: The cosmology object derived from Astropy, which is also set globally.

Power spectrum

This module implements models for the matter power spectrum, which can be evaluated using the matterPowerSpectrum() function. This module is automatically imported with the cosmology module.

Module contents

`PowerSpectrumModel`([output, allowed_types, ...])	Characteristics of power spectrum models.
`models`	Dictionary containing a list of models.
`powerSpectrum`(k, model, cosmo[, output])	The power spectrum (or transfer function) as a function of wavenumber.
`modelSugiyama95`(k, h, Om0, Ob0, Tcmb0)	The transfer function according to Sugiyama 1995.
`modelEisenstein98`(k, h, Om0, Ob0, Tcmb0)	The transfer function according to Eisenstein & Hu 1998.
`modelEisenstein98ZeroBaryon`(k, h, Om0, Ob0, ...)	The zero-baryon transfer function according to Eisenstein & Hu 1998.
`modelCamb`(k, cosmo[, ps_type, kmax])	The power spectrum as computed by the CAMB Boltzmann solver.

Module reference

cosmology.power_spectrum.CAMB_KMIN = 0.0001: The minimum wavenumber for which P(k) can be evaluated by CAMB.

cosmology.power_spectrum.CAMB_KMAX = 1000.0: The default maximum wavenumber for which P(k) can be evaluated by CAMB. The user can set a different upper limit, but that may increase the runtime significantly.

class cosmology.power_spectrum.PowerSpectrumModel(output=None, allowed_types=None, long_name='')

Characteristics of power spectrum models.

The models dictionary contains one item of this class for each available model.

Parameters

output: str: Indicates the quantity a given model outputs, which can be the transfer function (tf) or the power spectrum (ps).
allowed_types: array_like: List of strings that indicate different types of power spectra that can be returned. By default, the total PS is returned, but some models (i.e., Boltzmann codes) can also compute the spectra of components such as CDM only.

cosmology.power_spectrum.models = {'camb': <cosmology.power_spectrum.PowerSpectrumModel object>, 'eisenstein98': <cosmology.power_spectrum.PowerSpectrumModel object>, 'eisenstein98_zb': <cosmology.power_spectrum.PowerSpectrumModel object>, 'sugiyama95': <cosmology.power_spectrum.PowerSpectrumModel object>}

Dictionary containing a list of models.

An ordered dictionary containing one PowerSpectrumModel entry for each model.

cosmology.power_spectrum.powerSpectrum(k, model, cosmo, output='ps', **kwargs)

The power spectrum (or transfer function) as a function of wavenumber.

The transfer function transforms the spectrum of primordial fluctuations into the linear power spectrum of the matter density fluctuations. The primordial power spectrum is usually described as a power law, leading to a power spectrum of

\[P(k) = T(k)^2 k^{n_s}\]

where P(k) is the matter power spectrum, T(k) is the transfer function, and \(n_s\) is the tilt of the primordial power spectrum. See the Cosmology class for further details on the cosmological parameters.

Parameters

k: array_like: The wavenumber k (in comoving h/Mpc); can be a number or a numpy array.
model: str: The power spectrum model (see table above).
cosmo: cosmology: A Cosmology object.
output: str: Indicates which quantity should be returned, namely the transfer function (tf) or the power spectrum (ps).
kwargs: kwargs: Keyword arguments that are passed to the function evaluating the given model.

Returns

Pk: array_like: The power spectrum; has the same dimensions as k.

cosmology.power_spectrum.powerSpectrumModelName(model, **ps_args)

A unique internal name for the given power spectrum model (and parameters).

By default, the internal name used for power spectrum models is just the name of the model itself, but there are certain parameters that alter the spectrum so much that Colossus keeps track of separate spectra. In particular, the user can request the spectra of components such as dark matter and baryons (see powerSpectrum()), which are folded into the unique model name.

Parameters

model: str: The power spectrum model (see table above).
ps_args: kwargs: Keyword arguments that are passed to the function evaluating the given model. These arguments need to be consistent with those passed when the model is evaluated.

Returns

name: str: A unique name for this power spectrum model.

cosmology.power_spectrum.powerSpectrumLimits(model, **ps_args)

The lower and upper wavenumbers between which a model can be evaluated.

This function returns (None, None) for fitting functions that do not have a k-limit, and two floats for tabulated power spectra, Boltzmann codes, or other models that do have defined limits.

Parameters

model: str: The power spectrum model (see table above).
ps_args: kwargs: Keyword arguments that are passed to the function evaluating the given model. These arguments need to be consistent with those passed when the model is evaluated.

Returns

kmin: float: The lowest wavenumber (in comoving h/Mpc) where the model can be evaluated.
kmax: float: The highest wavenumber (in comoving h/Mpc) where the model can be evaluated.

cosmology.power_spectrum.modelCamb(k, cosmo, ps_type='tot', kmax=1000.0, **kwargs)

The power spectrum as computed by the CAMB Boltzmann solver.

This function translates a Colossus Cosmology object into parameters for the CAMB code and computes the power spectrum. See the CAMB documentation for information on possible keyword arguments and details about the calculations. We deliberately turn off the reionization component of the power spectrum, since that is usually not desired for large-scale structure and halo calculations.

In general, we leave as many parameters as possible to their default values in order to take advantage of future optimizations in the CAMB code. This means, on the other hand, that the results depend slightly on the code version. This function was tested with versions up to CAMB 1.3.5.

Important note: Colossus automatically keeps track of multiple versions of the power spectrum if different ps_type and kmax parameters are passed, but NOT for different keyword arguments (see kwargs below).

Parameters

k: array_like: The wavenumber k (in comoving h/Mpc); can be a number or a numpy array, but the input must be larger than 1E-4 (a fixed limit in CAMB) and smaller than CAMB_KMAX() unless a larger upper limit is passed as a kmax keyword argument.
cosmo: cosmology: A Cosmology object.
ps_type: str: CAMB can evaluate the tot (total) power spectrum (the default) or that of components, e.g. cdm (only dark matter) or baryon (only baryons). Other options can be passed by setting var1 = 'delta_XXX in the keyword arguments, but such a choice is not taken into account when internally naming the power spectrum model and will thus likely lead to inconsistent behavior. Conversely, Colossus can handle multiple power spectrum types within the same cosmology if ps_type is set.
kmax: float: The maximum wavenumber for which the CAMB calculation will be set up. This parameter is separate from the k array because kmax is fixed after CAMB has been initialized when this function is called for the first time. Note that increasing kmax significantly increases the runtime. On the other hand, it is up to the user to make sure that an insufficient kmax does not lead to inaccurate results. For example, when evaluating the variance sigma() for radii close to kmax, the variance will be underestimated due to the missing high-k modes in the power spectrum.
kwargs: kwargs: Arguments that are passed to the set_params() function in CAMB (see documentation). Note that Colossus does not keep track of different versions of the power spectra created with different kwargs. If multiple spectra are to be computed, it is easiest to create multiple cosmology objects.

Returns

P: array_like: The power spectrum in units of \(({\rm Mpc}/h)^3\); has the same dimensions as k.

cosmology.power_spectrum.transferFunction(k, h, Om0, Ob0, Tcmb0, model='eisenstein98')

The transfer function (deprecated).

This function is deprecated, powerSpectrum() should be used instead.

\[P(k) = T(k)^2 k^{n_s}\]

Parameters

k: array_like: The wavenumber k (in comoving h/Mpc); can be a number or a numpy array.
h: float: The Hubble constant in units of 100 km/s/Mpc.
Om0: float: \(\Omega_{\rm m}\), the matter density in units of the critical density at z = 0.
Ob0: float: \(\Omega_{\rm b}\), the baryon density in units of the critical density at z = 0.
Tcmb0: float: The temperature of the CMB at z = 0 in Kelvin.

Returns

Tk: array_like: The transfer function; has the same dimensions as k.

cosmology.power_spectrum.modelSugiyama95(k, h, Om0, Ob0, Tcmb0)

The transfer function according to Sugiyama 1995.

This function computes the Sugiyama 1995 approximation to the transfer function at a scale k, which is based on the Bardeen et al. 1986 formulation. Note that this approximation is not as accurate as the eisenstein98 model, with deviations of about 10-20% in the power spectrum, variance, and correlation function.

Parameters

k: array_like: The wavenumber k (in comoving h/Mpc); can be a number or a numpy array.
h: float: The Hubble constant in units of 100 km/s/Mpc.
Om0: float: \(\Omega_{\rm m}\), the matter density in units of the critical density at z = 0.
Ob0: float: \(\Omega_{\rm b}\), the baryon density in units of the critical density at z = 0.
Tcmb0: float: The temperature of the CMB at z = 0 in Kelvin.

Returns

Tk: array_like: The transfer function; has the same dimensions as k.

cosmology.power_spectrum.modelEisenstein98(k, h, Om0, Ob0, Tcmb0)

The transfer function according to Eisenstein & Hu 1998.

This function computes the Eisenstein & Hu 1998 approximation to the transfer function at a scale k. The code was adapted from Matt Becker’s cosmocalc code.

This function was tested against numerical calculations based on the CAMB code (Lewis et al. 2000) and found to be accurate to 5% or better up to k of about 100 h/Mpc (see the Colossus code paper for details).

Parameters

k: array_like: The wavenumber k (in comoving h/Mpc); can be a number or a numpy array.
h: float: The Hubble constant in units of 100 km/s/Mpc.
Om0: float: \(\Omega_{\rm m}\), the matter density in units of the critical density at z = 0.
Ob0: float: \(\Omega_{\rm b}\), the baryon density in units of the critical density at z = 0.
Tcmb0: float: The temperature of the CMB at z = 0 in Kelvin.

Returns

Tk: array_like: The transfer function; has the same dimensions as k.

See also

modelEisenstein98ZeroBaryon: The zero-baryon transfer function according to Eisenstein & Hu 1998.

cosmology.power_spectrum.modelEisenstein98ZeroBaryon(k, h, Om0, Ob0, Tcmb0)

The zero-baryon transfer function according to Eisenstein & Hu 1998.

This fitting function is significantly simpler than the full modelEisenstein98() version, and still approximates numerical calculations from a Boltzmann code to better than 10%, and almost as accurate when computing the variance or correlation function (see the Colossus code paper for details).

If Ob > 0, the assumptions of zero baryons is obviously inconsistent. However, the function executes without a warning because it is most commonly intended as a simplification of the power spectrum (e.g., to avoid the BAO peaks), not as an actual model for a Universe without baryons.

Parameters

k: array_like: The wavenumber k (in comoving h/Mpc); can be a number or a numpy array.
h: float: The Hubble constant in units of 100 km/s/Mpc.
Om0: float: \(\Omega_{\rm m}\), the matter density in units of the critical density at z = 0.
Ob0: float: \(\Omega_{\rm b}\), the baryon density in units of the critical density at z = 0.
Tcmb0: float: The temperature of the CMB at z = 0 in Kelvin.

Returns

Tk: array_like: The transfer function; has the same dimensions as k.

See also

modelEisenstein98: The transfer function according to Eisenstein & Hu 1998.

Cosmology

Previous topic

Units

Next topic

Large-scale structure

Cosmology

Basics

Setting and getting cosmologies

Changing and switching cosmologies

Converting to and from Astropy cosmologies

Summary of getter and setter functions

Standard cosmologies

Dark energy and curvature

Power spectrum, variance, correlation function

Derivatives and inverses

Performance optimization and accuracy

Module reference

Power spectrum

Module contents

Module reference

Table of Contents

Previous topic

Next topic