Testing the halo model¶
A different approach to compute the non-linear matter power spectrum is provided by the halo model, which became popular at the beginning of the 2000. This code implements all the useful functions for this purpose, using the Sheth-Tormen mass function. Despite it returns a poor description of matter clustering, it can be easily expanded to account for galaxies and/or qualitatively predicting covariance matrices for cosmological observables (for a review see arXiv:0206508 ).
Matter power spectrum¶
Let us start defining a colibri.halo.halo() instance:
import colibri.cosmology as cc
import colibri.halo as hc
# Define a halo instance
H = hc.halo(z = [0., 1., 2.], # 3 redshifts
k = np.logspace(-4., 1., 201), # Scales in h/Mpc
code = 'camb', # Boltzmann code with which to compute linear P(k)
BAO_smearing = False, # Smooth BAO due to non-linearities
cosmology = cc.cosmo()) # ``cosmo`` instance
The halo power spectrum is then computed with the following lines
# Load power spectrum according to the halo model
H.halo_Pk(kwargs_mass_function = {'a': 0.707, 'p': 0.3}, # arguments to pass to the Sheth-Tormen mass function
kwargs_concentration = {'c0': 9., 'b': 0.13}) # arguments to pass to the concentration parameter c0*(M/Mstar)**(-b)
# Store 1-halo, 2-halo and total terms
oneh = H.Pk['matter']['1-halo'] # array of shape (len(z),len(k))
twoh = H.Pk['matter']['2-halo'] # array of shape (len(z),len(k))
tot = H.Pk['matter']['total halo'] # array of shape (len(z),len(k))
The comparison with CAMB Halofit (for our 3 different redshifts) is the following:
Galaxy power spectrum: real space¶
The galaxy power spectrum in real space is computed using the Halo Occupation Distribution (HOD) method. In this pictures, halos are populated of galaxies depending exclusively on their mass. Galaxies are divided into central (max. 1 per halo) and satellites (from 0 to infinity): the 1-halo and 2-halo terms are computed by counting pairs of galaxies in the same and different halos, respectively.
To compute the power spectrum, one must provide the average amount of galaxies per halo.
This is done in a few steps here.
First, we define a colibri.galaxy.galaxy() instance:
import colibri.cosmology as cc
import colibri.galaxy as gc
# Define a halo instance
G = gc.galaxy(z = [0., 1., 2.], # 3 redshifts
k = np.logspace(-4., 2., 201), # Scales in h/Mpc
code = 'camb', # Boltzmann code with which to compute linear P(k)
BAO_smearing = False, # Smooth BAO due to non-linearities
cosmology = cc.cosmo()) # ``cosmo`` instance
Second, we give to the code the ‘ingredients’ for the HOD.
This is done through the function colibri.galaxy.galaxy.load_HOD()
G.load_HOD(kind_satellite = G.power_law,
kwargs_satellite = {'log_Mcut':[13., 12.5, 12.3], 'log_M1': [13., 13.4, 13.], 'alpha': [1., 1.5, 1.]},
kind_central = G.logistic_function,
kwargs_central = {'log_Mmin': [13., 12.4, 11.4], 'sigma_logM': [0.8, 0.5, 0.6]})
The arguments kind_satellite and kind_central are callable functions whose first argument must be the halo mass in units of \(M_\odot/h\).
In this case we have used the functions colibri.galaxy.galaxy.power_law() and colibri.galaxy.galaxy.logistic_function().
The arguments kwargs_satellite and kwargs_central are dictionaries that contain all the remaining arguments to pass to the above functions.
Each value of each key must be a list of same size of the redshifts required, otherwise the code will return an AssertionError.
The last step consists of the actual computation of the galaxy power spectrum, which is done like in the halo case:
G.galaxy_Pk(kwargs_mass_function = {'a': 0.707, 'p': 0.3},
kwargs_concentration = {'c0': 9., 'b': 0.13})
pg1 = G.Pk['galaxies']['real space']['1-halo']
pg2 = G.Pk['galaxies']['real space']['2-halo']
pgt = G.Pk['galaxies']['real space']['total halo']
Note
The last two steps can be combined together.
In fact, when the function colibri.galaxy.galaxy.galaxy_Pk() is called, it will search the quantities self.Ncen and self.Nsat which are generated only after the call of colibri.galaxy.galaxy.load_HOD().
If these are not found, the code will try to compute them on the fly, provided that the arguments of colibri.galaxy.galaxy.load_HOD() are given.
In other words, the last two steps can be gathered in:
G.galaxy_Pk(kind_satellite = G.power_law,
kwargs_satellite = {'log_Mcut':[13., 12.5, 12.3], 'log_M1': [13., 13.4, 13.], 'alpha': [1., 1.5, 1.]},
kind_central = G.logistic_function,
kwargs_central = {'log_Mmin': [13., 12.4, 11.4], 'sigma_logM': [0.8, 0.5, 0.6]},
kwargs_mass_function = {'a': 0.707, 'p': 0.3},
kwargs_concentration = {'c0': 9., 'b': 0.13})
Galaxy power spectrum: redshift space¶
Galaxies in surveys are observed in redshift-space due to the fact that we only measure the recession velocity along the line-of-sight.
The class colibri.RSD.RSD() provides routines able to compute the power spectrum in redshift-space using different configurations of independent variables:
in \((k,\mu)\) space
in multipole space, using the Legendre expansion
in \((k_\parallel,k_\perp)\) space
The RSD instance is called as follows:
import colibri.cosmology as cc
import colibri.RSD as rsd
Z = rsd.RSD(z = [0., 1., 2.], # Redshifts
k = np.geomspace(0.0005, 10., 101), # Scales in h/Mpc
mu = np.linspace(0., 1., 31), # Cosine of angles with LOS
k_par = np.linspace(0.01, 1., 51), # Scales parallel in h/Mpc
k_perp = np.linspace(0.01, 1., 31), # Scales perpendicular in h/Mpc
BAO_smearing = True, # Smooth BAO feature in non-linearities
cosmology = cc.cosmo()) # Cosmology
With our 3 redshifts, we define the following useful quantities
# Galaxy biases
bb = [1.30, 2.60, 3.34] # Galaxy biases
# Growth rates
ff = C.Omega_m_z(zz)**0.55
# Velocity dispersions (in km/s)
sv = [200., 200., 200.]
# HOD functions (like in ``galaxy``)
HOD_central_kind, HOD_satellite_kind = Z.logistic_function, Z.power_law
HOD_central_parameters = {'log_Mmin': [12., 12.5, 12.], 'sigma_logM': [0.8, 0.5, 0.2]}
HOD_satellite_parameters = {'log_Mcut': [13., 12.5, 12.], 'log_M1': [14., 13.2, 13.4], 'alpha': [1., 1.25, 1.5]}
# Kind of damping functions for fingers-of-god effect
FoG_damping = 'Lorentzian' # Kind of damping due to fingers-of-God (choose between 'Lorentzian' and 'Gaussian')
# Base power spectrum is the non-linear matter one (computed with Halofit) (choose among 'linear', 'nonlinear', 'HOD', 'halo model')
RSD_model = 'nonlinear'
At this point the galaxy redshift-space power spectrum can be computed in 3 ways:
in \((k,\mu)\) space
Z.galaxy_RSD_Pk(bias = bb, # Galaxy bias (used only if model = 'HOD' or 'halo model') growth_rate = ff, # Growth rate f = dln(D)/dln(a) velocity_dispersion = sv, # Average velocity dispersion of galaxies in halos model = RSD_model, # Model to compute RSD kwargs_mass_function = {'a': 0.707, 'p': 0.3}, # Parameters to compute halo mass function (used only if model = 'HOD' or 'halo model') kwargs_concentration = {'c0': 9., 'b': 0.13}, # Parameters to compute concentration parameter (used only if model = 'HOD' or 'halo model') fingers_of_god = FoG_damping, # Kind of damping ('Lorentzian' or 'Gaussian', used only if model != 'halo model') kind_central = HOD_central_kind, # Function to compute central galaxies (1st arguments must be mass in Msun/h) kwargs_central = HOD_central_parameters, # Remaining arguments to pass to kind_central kind_satellite = HOD_satellite_kind, # Function to compute satellite galaxies (1st arguments must be mass in Msun/h) kwargs_satellite = HOD_satellite_parameters) # Remaining arguments to pass to kind_satellite
in multipole space, using the Legendre expansion
Z.galaxy_RSD_Pk_multipoles(l = [0,2,4] # Multipoles to compute (monopole, quadrupole, hexadecapole) bias = bb, # Galaxy bias (used only if model = 'HOD' or 'halo model') growth_rate = ff, # Growth rate f = dln(D)/dln(a) velocity_dispersion = sv, # Average velocity dispersion of galaxies in halos model = RSD_model, # Model to compute RSD kwargs_mass_function = {'a': 0.707, 'p': 0.3}, # Parameters to compute halo mass function (used only if model = 'HOD' or 'halo model') kwargs_concentration = {'c0': 9., 'b': 0.13}, # Parameters to compute concentration parameter (used only if model = 'HOD' or 'halo model') fingers_of_god = FoG_damping, # Kind of damping ('Lorentzian' or 'Gaussian', used only if model != 'halo model') kind_central = HOD_central_kind, # Function to compute central galaxies (1st arguments must be mass in Msun/h) kwargs_central = HOD_central_parameters, # Remaining arguments to pass to kind_central kind_satellite = HOD_satellite_kind, # Function to compute satellite galaxies (1st arguments must be mass in Msun/h) kwargs_satellite = HOD_satellite_parameters) # Remaining arguments to pass to kind_satellite
in \((k_\parallel,k_\perp)\) space
Z.galaxy_RSD_Pk_2D(bias = bb, # Galaxy bias (used only if model = 'HOD' or 'halo model') growth_rate = ff, # Growth rate f = dln(D)/dln(a) velocity_dispersion = sv, # Average velocity dispersion of galaxies in halos model = RSD_model, # Model to compute RSD kwargs_mass_function = {'a': 0.707, 'p': 0.3}, # Parameters to compute halo mass function (used only if model = 'HOD' or 'halo model') kwargs_concentration = {'c0': 9., 'b': 0.13}, # Parameters to compute concentration parameter (used only if model = 'HOD' or 'halo model') fingers_of_god = FoG_damping, # Kind of damping ('Lorentzian' or 'Gaussian', used only if model != 'halo model') kind_central = HOD_central_kind, # Function to compute central galaxies (1st arguments must be mass in Msun/h) kwargs_central = HOD_central_parameters, # Remaining arguments to pass to kind_central kind_satellite = HOD_satellite_kind, # Function to compute satellite galaxies (1st arguments must be mass in Msun/h) kwargs_satellite = HOD_satellite_parameters) # Remaining arguments to pass to kind_satellite
The result should be the following:
