Angular spectra in Limber’s approximation

For a usage example of this library, see Power spectra in Limber’s approximation

class colibri.limber.limber(z_limits=(0.0, 5.0), cosmology=<colibri.cosmology.cosmo object>)

The class colibri.limber.limber() contains all the functions useful to compute the angular power spectra and correlation functions in the flat sky and Limber’s approximation. It also contains routines to compute widely-used window functions as well as a routine to add custom ones. At initialization it takes as inputs a redshift range for integration and a colibri.cosmology.cosmo() instance.

Parameters:

cosmology (cosmo instance, default = cosmology.cosmo()) – Fixes the cosmological parameters. If not declared, the default values are chosen (see colibri.cosmology.cosmo() documentation).
z_limits (2-uple or list/array of length 2, default = (0., 5.)) – Lower and upper limit of integration along the line of sight. Both numbers must be non-negative and the first number must be smaller than the second. If the lower limit is set to 0, it will be enhanced by 1e-10 to avoid divergences at the origin of the lightcone.

Warning

All the power spectra are computed in the Limber approximation and the window function are assumed to be dependent only on redshift and not on scales (see e.g. colibri.limber.limber.load_lensing_window_functions()). Typically the scale dependence of the window functions can be factorized out (e.g. ISW effect, different orders of cosmological perturbation theory…) and in this code it can be added to the power spectrum (see colibri.limber.limber.load_power_spectra()).

constant_distribution(zmin, zmax, step=0.005)

Example function for the distribution of source galaxy. Here we use a constant distribution of sources.

Parameters:

z (array) – Redshifts.
zmin (float) – Lower edge of the bin.
zmax (float) – Upper edge of the bin.
step (float, default = 0.005) – width of the cutoff (better to avoid a sharp one for numerical reasons, better set it to be at least 0.001)

Returns:

array

euclid_distribution(z, zmin, zmax, a=2.0, b=1.5, z_med=0.9, step=0.005)

Example function for the distribution of source galaxy. This distribution in particular is expected to be used in the Euclid mission:

\[n(z) \propto z^a \ \exp{\left[-\left(\frac{z}{z_{med}/\sqrt 2}\right)^b\right]}\]

This distribution will eventually be normalized such that its integral on all redshifts is 1.

Parameters:

z (array) – Redshifts.
zmin (float) – Lower edge of the bin (a small width of 0.005 will be applied for convergence reasons).
zmax (float) – Upper edge of the bin (a small width of 0.005 will be applied for convergence reasons).
a (float, default = 2.0) – Parameter of the distribution.
b (float, default = 1.5) – Parameter of the distribution.
z_med (float, default = 0.9) – Median redshift of the distribution.
step (float, default = 0.005) – width of the cutoff (better to avoid a sharp one for numerical reasons, better set it to be at least 0.001)

Returns:

array

euclid_distribution_with_photo_error(z, zmin, zmax, a=2.0, b=1.5, z_med=0.9, f_out=0.1, c_b=1.0, z_b=0.0, sigma_b=0.05, c_o=1.0, z_o=0.1, sigma_o=0.05)

Example function for the distribution of source galaxy. This distribution in particular is expected to be used in the Euclid mission. Here also the effect of photometric errors is included.

\[n^{(i)}(z) \propto \int_{z_i^-}^{z_i^+} dy \ z^a \ \exp{\left[-\left(\frac{z}{z_{med}/\sqrt 2}\right)^b\right]} \ p_\mathrm{ph}(y|z)\]

where

\[ \begin{align}\begin{aligned}p_\mathrm{ph}(y|z) = \frac{1-f_\mathrm{out}}{\sqrt{2\pi}\sigma_b(1+z)} \ \exp\left[-\frac{1}{2} \left(\frac{z-c_b y -z_b}{\sigma_b(1+z)}\right)^2\right] +\\+ \frac{f_\mathrm{out}}{\sqrt{2\pi}\sigma_o(1+z)} \ \exp\left[-\frac{1}{2} \left(\frac{z-c_o y -z_o}{\sigma_o(1+z)}\right)^2\right]\end{aligned}\end{align} \]

Parameters:

z (array) – Redshifts.
zmin (float) – Lower edge of the bin.
zmax (float) – Upper edge of the bin.
a (float, default = 1.5) – Parameter of the distribution.
b (float, default = 1.5) – Parameter of the distribution.
z_med (float, default = 0.9) – Median redshift of the distribution.
f_out (float, default = 0.1) – Fraction of outliers
c_b (float, default = 1.0) – Parameter of the Gaussian (normalization) representing the uncertainty on the photometric error for in-liers.
z_b (float, default = 0.0) – Parameter of the Gaussian (scale-redshift) representing the uncertainty on the photometric error for in-liers.
sigma_b (float, default = 0.05) – Parameter of the Gaussian (width) representing the uncertainty on the photometric error for in-liers.
c_o (float, default = 1.0) – Parameter of the Gaussian (normalization) representing the uncertainty on the photometric error for out-liers.
z_o (float, default = 0.1) – Parameter of the Gaussian (scale-redshift) representing the uncertainty on the photometric error for out-liers.
sigma_o (float, default = 0.05) – Parameter of the Gaussian (width) representing the uncertainty on the photometric error for out-liers.

Returns:

array

gaussian_distribution(z, mean, sigma)

Example function for the distribution of source galaxy. Here we use a Gaussian galaxy distribution

Parameters:

z (array) – Redshifts.
mean (float) – Mean redshift of the distribution.
sigma (float) – Width of the Gaussian

Returns:

array

geometric_factor_f_K(z, z0=0.0)

Geometric factor (distance) between two given redshifts z and z0. It assumes neutrinos as matter, which is a good approximation at low redshifts. In fact, this latter assumption introduces a bias of less than 0.02% at \(z<10\) for even the lowest neutrino masses allowed by particle physics.

Parameters:

z (array) – Redshifts.
z0 (float, default = 0) – Pivot redshift.

Returns:

array

limber_angular_correlation_functions(theta, l, Cl, order)

This function computes the angular correlation function from an angular power spectrum. The equation is as follows

\[\xi^{(ij)}_{XY}(\theta) = \int_0^\infty \frac{d\ell}{2\pi} \ \ell \ C_{XY}^{(ij)}(\ell) \ J_{\nu} (\ell\theta),\]

where \(\nu\) is the order of the transform and it changes from observable to observable.

Warning

For example, for shear \(\nu=0\) or \(\nu=4\), for galaxy clustering \(\nu=0\) and for galaxy-galaxy lensing \(\nu=2\).

Parameters:

theta (array) – Angles (in \(\mathrm{arcmin}\)) where to compute the shear correlation functions
l (array) – Multipoles at which the spectrum is computed
Cl (3D array) – 3D array, where first and second dimensions are the bins and the third is the multipoles, i.e. Cl[bin i, bin j, multipole l]. The last dimension has to have the same length as l.
order (float) – Order of Hankel transform.

Returns:

3D array containing xi[bin i, bin j, angle theta]

limber_angular_power_spectra(l, windows=None)

This function computes the angular power spectra (using the Limber’s and the flat-sky approximations) for the window function specified. Given two redshift bins i and j the equation is

\[C^{(ij)}(\ell) = \int_0^\infty dz \ \frac{c}{H(z)} \ \frac{W^{(i)}(z) W^{(j)}(z)}{f_K^2[\chi(z)]} \ P\left(\frac{\ell}{f_K[\chi(z)]}, z\right),\]

where \(P(k,z)\) is the matter power spectrum and \(W^{(i)}(z)\) are the window functions.

Parameters:

l (array) – Multipoles at which to compute the shear power spectra.
windows (list of strings, default = None) – which spectra (auto and cross) must be computed. If set to None all the spectra will be computed.

Returns:

dictionary whose keys are combinations of window functions specified in windows. Each key is a 3-D array whose entries are Cl[bin i, bin j, multipole l]

load_CMB_lensing_window_functions(z, nz, z_LSS=1089.0, name='CMB lensing')

This function computes the window function for CMB lensing given a galaxy distribution. The function automatically normalizes the galaxy distribution such that the integral over redshifts is 1. The routine adds a key (specified in the name argument) to the self.window_function dictionary. Given a galxy distruibution in a redshift bin \(n^{(i)}(z)\), the equation is:

\[W^{(i)}_\mathrm{CMB}(z) = \frac{3}{2}\Omega_m \ \frac{H_0^2}{c^2} \ f_K[\chi(z)] (1+z) \ n^{(i)}(z) \frac{H(z)}{c} \ \frac{f_K[\chi(z_{LSS})]-f_K[\chi(z)]}{f_K[\chi(z_{LSS})]}\]

Parameters:

z (1-D array, default = None) – array or list of redshift at which the galaxy distribution nz is evaluated
nz (2-D array with shape (n_bins, len(z)), default = None) – 2-D array or 2-D list where each sublist is the galaxy distribution of a given redshift bin
z_LSS (float, default = 1089) – last-scattering surface redshift.
name (string, default = 'CMB lensing') – name of the key to add to the dictionary

An example call can be, for 3 bins all with a colibri.limber.limber.euclid_distribution() with default arguments for a and b but different bin edges zmin, zmax:

bin_edges = [0.00, 0.72, 1.11, 5.00]
nbins     = len(bin_edges)-1
z_w       = np.linspace(0., 6., 1001)
nz_w      = [S.euclid_distribution(z = z_w, a = 2.0, b = 1.5, zmin = bin_edges[i], zmax = bin_edges[i+1]) for i in range(nbins)]
S.load_CMB_window_functions(z = z_w, nz = nz_w, z_LSS = 1089.)

Returns:: A key of a given name is added to the self.window_function dictionary

load_HI_window_functions(z, nz, Omega_HI=0.000625, bias=1.0, name='HI')

This function computes the window function for HI brightness temperature given a galaxy distribution. The function automatically normalizes the galaxy distribution such that the integral over redshifts is 1. The routine adds a key (specified in the name argument) to the self.window_function dictionary. Given a galxy distruibution in a redshift bin \(n^{(i)}(z)\), the equation is:

\[W^{(i)}_\mathrm{HI}(z) = b(z) \ T_b(z) \ D(z) \ \ n^{(i)}(z) \frac{H(z)}{c}\]

Parameters:

z (1-D array, default = None) – array or list of redshift at which the galaxy distribution nz is evaluated
nz (2-D array with shape (n_bins, len(z)), default = None) – 2-D array or 2-D list where each sublist is the galaxy distribution of a given redshift bin
Omega_HI (float, default = 6.25e-4) – HI density parameter.
bias (float or array, same length of nz, default = 1) – Galaxy bias.
name (string, default = 'HI') – name of the key to add to the dictionary

An example call can be, for 3 bins all with a colibri.limber.limber.euclid_distribution() with default arguments for a and b but different bin edges zmin, zmax:

bin_edges = [0.00, 0.72, 1.11, 5.00]
nbins     = len(bin_edges)-1
z_w       = np.linspace(0., 6., 1001)
nz_w      = [S.euclid_distribution(z = z_w, a = 2.0, b = 1.5, zmin = bin_edges[i], zmax = bin_edges[i+1]) for i in range(nbins)]
S.load_HI_window_functions(z = z_w, nz = nz_w, Omega_HI = 0.00063, bias = 1)

Returns:: A key of a given name is added to the self.window_function dictionary

load_IA_window_functions(z, nz, A_IA=1.0, eta_IA=0.0, beta_IA=0.0, lum_IA=1.0, name='IA')

This function computes the window function for intrinsic alignment given a galaxy distribution. The function automatically normalizes the galaxy distribution such that the integral over redshifts is 1. The routine adds a key (specified in the name argument) to the self.window_function dictionary. Given a galxy distruibution in a redshift bin \(n^{(i)}(z)\), the equation is:

\[W^{(i)}_\mathrm{IA}(z) = -\frac{A_\mathrm{IA} \mathcal C_1 \Omega_\mathrm m}{D_1(k,z)}(1+z)^{\eta_\mathrm{IA}} \left[\frac{L(z)}{L_*(z)}\right]^{\beta_\mathrm{IA}} \ n^{(i)}(z) \frac{H(z)}{c}\]

Parameters:

z (1-D array, default = None) – array or list of redshift at which the galaxy distribution nz is evaluated
nz (2-D array with shape (n_bins, len(z)), default = None) – 2-D array or 2-D list where each sublist is the galaxy distribution of a given redshift bin
A_IA (float, default = 1) – Intrinsic alignment amplitude.
eta_IA (float, default = 0) – Exponent for redshift dependence of intrinsic alignment.
beta_IA (float, default = 0) – Exponent for luminosity dependence of intrinsic alignment.
lum_IA (float or callable whose only argument is \(z\), default = 1) – Relative luminosity of galaxies \(L(z)/L_*(z)\).
name (string, default = 'IA') – name of the key to add to the dictionary

An example call can be, for 3 bins all with a colibri.limber.limber.euclid_distribution() with default arguments for a and b but different bin edges zmin, zmax:

bin_edges = [0.00, 0.72, 1.11, 5.00]
nbins     = len(bin_edges)-1
z_w       = np.linspace(0., 6., 1001)
nz_w      = [S.euclid_distribution(z = z_w, a = 2.0, b = 1.5, zmin = bin_edges[i], zmax = bin_edges[i+1]) for i in range(nbins)]
S.load_IA_window_functions(z = z_w, nz = nz_w, A_IA = 1, eta_IA = 0, beta_IA = 0, lum_IA = 1)

Returns:: A key of a given name is added to the self.window_function dictionary

load_custom_window_functions(z, window, name)

This function loads a custom window function and adds the key to the dictionary The window function in input must already be normalized.

Parameters:

z (1-D array, default = None) – array or list of redshift at which the galaxy distribution nz is evaluated
window (2-D array with shape (n_bins, len(z)), default = None) – 2-D array or 2-D list where each sublist is the galaxy distribution of a given redshift bin
name (string) – name of the key to add to the dictionary

Returns:

A key of a given name is added to the self.window_function dictionary

load_galaxy_clustering_window_functions(z, nz, bias=1.0, name='galaxy')

This function computes the window function for galaxy clustering given a galaxy distribution. The function automatically normalizes the galaxy distribution such that the integral over redshifts is 1. The routine adds a key (specified in the name argument) to the self.window_function dictionary. Given a galxy distruibution in a redshift bin \(n^{(i)}(z)\), the equation is:

\[W^{(i)}_\mathrm{G}(z) = b(z) \ n^{(i)}(z) \frac{H(z)}{c}\]

Parameters:

z (1-D array, default = None) – array or list of redshift at which the galaxy distribution nz is evaluated
nz (2-D array with shape (n_bins, len(z)), default = None) – 2-D array or 2-D list where each sublist is the galaxy distribution of a given redshift bin
bias (float or array, same length of nz, default = 1) – Galaxy bias.
name (string, default = 'galaxy') – name of the key to add to the dictionary

An example call can be, for 3 bins all with a colibri.limber.limber.euclid_distribution() with default arguments for a and b but different bin edges zmin, zmax:

bin_edges = [0.00, 0.72, 1.11, 5.00]
nbins     = len(bin_edges)-1
z_w       = np.linspace(0., 6., 1001)
nz_w      = [S.euclid_distribution(z = z_w, a = 2.0, b = 1.5, zmin = bin_edges[i], zmax = bin_edges[i+1]) for i in range(nbins)]
S.load_galaxy_clustering_window_functions(z = z_w, nz = nz_w, bias = 1)

Returns:: A key of a given name is added to the self.window_function dictionary

load_lensing_window_functions(z, nz, A_IA=0.0, eta_IA=0.0, beta_IA=0.0, lum_IA=1.0, name='lensing')

This function computes the window function for lensing (comprehensive of shear and intrinsic alignment) given a galaxy distribution. The function automatically normalizes the galaxy distribution such that the integral over redshifts is 1. The routine adds a key (specified in the name argument) to the self.window_function dictionary. See colibri.limber.limber.load_shear_window_functions() and colibri.limber.limber.load_IA_window_functions() for the equations.

Parameters:

z (1-D array, default = None) – array or list of redshift at which the galaxy distribution nz is evaluated
nz (2-D array with shape (n_bins, len(z)), default = None) – 2-D array or 2-D list where each sublist is the galaxy distribution of a given redshift bin
A_IA (float, default = 1) – Intrinsic alignment amplitude.
eta_IA (float, default = 0) – Exponent for redshift dependence of intrinsic alignment.
beta_IA (float, default = 0) – Exponent for luminosity dependence of intrinsic alignment.
lum_IA (float or callable whose only argument is \(z\), default = 1) – Relative luminosity of galaxies \(L(z)/L_*(z)\).
name (string, default = 'lensing') – name of the key to add to the dictionary

An example call can be, for 3 bins all with a colibri.limber.limber.euclid_distribution() with default arguments for a and b but different bin edges zmin, zmax:

bin_edges = [0.00, 0.72, 1.11, 5.00]
nbins     = len(bin_edges)-1
z_w       = np.linspace(0., 6., 1001)
nz_w      = [S.euclid_distribution(z = z_w, a = 2.0, b = 1.5, zmin = bin_edges[i], zmax = bin_edges[i+1]) for i in range(nbins)]
S.load_lensing_window_functions(z = z_w, nz = nz_w, A_IA = 1, eta_IA = 0, beta_IA = 0, lum_IA = 1)

Returns:: A key of a given name is added to the self.window_function dictionary

load_power_spectra(k, z, power_spectra)

This routine interpolates the total matter power spectrum (using the CDM prescription) in scales (units of \(h/\mathrm{Mpc}\)) and redshifts. power_spectra and galaxy_bias must be a 2D array of shape (len(z), len(k)) which contains the power spectrum (in units of \((\mathrm{Mpc}/h)^3\)) and galaxy bias evaluated at the scales and redshifts specified above.

Parameters:

k (array) – Scales in units of \(h/\mathrm{Mpc}\).
z (array) – Redshifts at which power spectrum must be/is computed.
power_spectra (2D NumPy array) – It must be a 2D array of shape (len(z), len(k)) which contains the power spectrum (in units of \((\mathrm{Mpc}/h)^3\)) evaluated at the scales and redshifts specified above.

Returns:

Nothing, but two 2D-interpolated object self.power_spectra_interpolator and self.galaxy_bias_interpolator containing \(P(k,z)\) in units of \((\mathrm{Mpc}/h)^3\) and \(b(k,z)\) are created

load_shear_window_functions(z, nz, name='shear')

This function computes the window function for cosmic shear given the galaxy distribution in input. The function automatically normalizes the galaxy distribution such that the integral over redshifts is 1. The routine adds a key (specified in the name argument) to the self.window_function dictionary. Given a galxy distruibution in a redshift bin \(n^{(i)}(z)\), the equation is:

\[W^{(i)}_\gamma(z) = \frac{3}{2}\Omega_m \ \frac{H_0^2}{c^2} \ f_K[\chi(z)] (1+z) \int_z^\infty dx \ n^{(i)}(x) \ \frac{f_K[\chi(z-x)]}{f_K[\chi(z)]}\]

Parameters:

z (1-D array, default = None) – array or list of redshift at which the galaxy distribution nz is evaluated
nz (2-D array with shape (n_bins, len(z)), default = None) – 2-D array or 2-D list where each sublist is the galaxy distribution of a given redshift bin
name (string, default = 'shear') – name of the key to add to the dictionary

An example call can be, for 3 bins all with a colibri.limber.limber.euclid_distribution() with default arguments for a and b but different bin edges zmin, zmax:

bin_edges = [0.00, 0.72, 1.11, 5.00]
nbins     = len(bin_edges)-1
z_w       = np.linspace(0., 6., 1001)
nz_w      = [S.euclid_distribution(z = z_w, a = 2.0, b = 1.5, zmin = bin_edges[i], zmax = bin_edges[i+1]) for i in range(nbins)]
S.load_shear_window_functions(z = z_w, nz = nz_w)

Returns:: A key of a given name is added to the self.window_function dictionary