Angpow: a software for the fast computation of accurate tomographic power spectra

TL;DR

Angpow is a software tool designed to rapidly and accurately compute angular power spectra between redshift bins, facilitating cosmological analyses without relying on the Limber approximation, especially useful for large parameter space explorations.

Contribution

This paper introduces Angpow, a novel software that efficiently computes angular power spectra using Chebyshev polynomial basis and Clenshaw-Curtis quadrature, improving over traditional methods.

Findings

Angpow achieves fast computation of $C_\ell(z_1,z_2)$ spectra.

It accurately handles auto and cross-correlations without Limber approximation.

Validated results match other codes and demonstrate high performance.

Abstract

The statistical distribution of galaxies is a powerful probe to constrain cosmological models and gravity. In particular the matter power spectrum $P(k)$ brings information about the cosmological distance evolution and the galaxy clustering together. However the building of $P(k)$ from galaxy catalogues needs a cosmological model to convert angles on the sky and redshifts into distances, which leads to difficulties when comparing data with predicted $P(k)$ from other cosmological models, and for photometric surveys like LSST. The angular power spectrum $C_\ell(z_1,z_2)$ between two bins located at redshift $z_1$ and $z_2$ contains the same information than the matter power spectrum, is free from any cosmological assumption, but the prediction of $C_\ell(z_1,z_2)$ from $P(k)$ is a costly computation when performed exactly. The Angpow software aims at computing quickly and accurately the auto ($z_1=z_2$) and cross ($z_1 \neq z_2$) angular power spectra between redshift bins. We describe the developed algorithm, based on developments on the Chebyshev polynomial basis and on the Clenshaw-Curtis quadrature method. We validate the results with other codes, and benchmark the performance. Angpow is flexible and can handle any user defined power spectra, transfer functions, and redshift selection windows. The code is fast enough to be embedded inside programs exploring large cosmological parameter spaces through the $C_\ell(z_1,z_2)$ comparison with data. We emphasize that the Limber's approximation, often used to fasten the computation, gives wrong $C_\ell$ values for cross-correlations.

Tables

Table 1. Table 1: Wall time (in seconds) measured at CCIN2P3 (on Intel Xeon CPU E5-2640 v3 processors) for the test benches described in the text, according to the number of OpenMP threads used. Results are given for the intel icpc 15.0 and gcc 5.2 compilers.

# Threads	1	2	4	8	16
Linux/icpc
Test 1	0.38	0.21	0.13	0.09	0.08
Test 2	0.76	0.41	0.23	0.15	0.11
Test 3	3.72	1.96	1.05	0.64	0.44
Test 4	9.97	5.25	2.79	1.60	1.01
Linux/gcc
Test 1	0.56	0.30	0.17	0.12	0.09
Test 2	1.14	0.60	0.33	0.20	0.14
Test 3	5.01	2.59	1.38	0.81	0.50
Test 4	13.80	7.07	3.71	2.12	1.27

Equations

$$C_{\ell }(z_1,z_2)=\frac{2}{\pi }{\int }_0^{\infty }\mathrm{d}k{k}^{2}P(k){\Delta }_{\ell }(z_1,k){\Delta }_{\ell }(z_2,k),$$

$${\Delta }_{\ell }^{mat.}(z,k)=bD(z,k)j_{\ell }(kr(z)),$$

$$\begin{array}{rl}{C}_{\ell }^{\mathrm{thick}}& (z_1,z_2;{\sigma }_1,{\sigma }_2)\\ & =\frac{2}{\pi }{\iint }_0^{\infty }\mathrm{d}z\mathrm{d}{z}^{\prime }{W}_1(z;z_1,{\sigma }_1)W_2(z^{\prime };z_2,{\sigma }_2)\\ & \times {\int }_0^{\infty }\mathrm{d}k{k}^{2}P(k){\Delta }_{\ell }(z,k){\Delta }_{\ell }(z^{\prime },k)\end{array}.$$

$$\begin{array}{rl}{f}_{\ell }(z,k)& \equiv \sqrt{\frac{2}{\pi }}k\sqrt{P(k)}{\Delta }_{\ell }(z,k)\\ & =\sqrt{\frac{2}{\pi }}k\sqrt{P(k)}\left\{bD(z,k)j_{\ell }(kr(z))+\dots \right\}\\ & =\sqrt{\frac{2}{\pi }}k\sqrt{P(k)D(z,k{)}^2}\left\{b{j}_{\ell }(kr(z))+\dots \right\}\\ & \equiv \sqrt{\frac{2}{\pi }}k\sqrt{P(k,z)}{\tilde{\Delta }}_{\ell }(z,k)\end{array},$$

$$\begin{array}{rl}{C}_{\ell }^{\mathrm{thick}}& (z_1,z_2;{\sigma }_1,{\sigma }_2)\\ & ={\iint }_0^{\infty }\mathrm{d}z\mathrm{d}{z}^{\prime }{W}_1(z;z_1,{\sigma }_1)W_2(z^{\prime };z_2,{\sigma }_2)\\ & \times {\int }_0^{\infty }\mathrm{d}k{f}_{\ell }(z,k)f_{\ell }(z^{\prime },k).\end{array}$$

$$\begin{array}{rl}{C}_{\ell }^{\mathrm{thick}}(z_1& ;{\sigma }_1)={\iint }_0^{\infty }\mathrm{d}z\mathrm{d}{z}^{\prime }W(z;z_1,{\sigma }_1)W(z^{\prime };z_1,{\sigma }_1)\\ & \times {\int }_0^{\infty }\mathrm{d}k{f}_{\ell }(z,k)f_{\ell }(z^{\prime },k).\end{array}$$

$$C_{\ell }^{\delta }(z_1)=\frac{2}{\pi }{\int }_0^{\infty }\mathrm{d}k{k}^{2}P(k,z_1){\tilde{\Delta }}_{\ell }^2(z_1,k).$$

$$C_{\ell }^{\mathrm{thick}}(z_1,z_2)\approx \sum _{i=0}^{N_{{\mathrm{z}}_1}-1}\sum _{j=0}^{N_{{\mathrm{z}}_2}-1}{w}_i{w}_j{W}_1(z_i,z_1)W_2(z_j,z_2){\hat{P}}_{\ell }(r_i,r_j)$$

$${\hat{P}}_{\ell }(z_i,z_j)={\int }_0^{\infty }\mathrm{d}k{f}_{\ell }(z_i,k)f_{\ell }(z_j,k),$$

$${\hat{P}}_{\ell }(z_i,z_j)\approx \sum _{p=0}^{N_{\mathrm{k}}-1}{I}_{\ell }(k_p^{\ell },k_{p+1}^{\ell };z_i,z_j),$$

$$\begin{array}{rl}{C}_{\ell }^{\mathrm{thick}}(z_1,z_2)\approx \sum _{i=0}^{N_{{\mathrm{z}}_1}}\sum _{j=0}^{N_{{\mathrm{z}}_2}}& w_i{w}_j{W}_1(z_i,z_1)W_2(z_j,z_2)\\ & \times \sum _{p=0}^{N_{\mathrm{k}}-1}{I}_{\ell }(k_p^{\ell },k_{p+1}^{\ell };z_i,z_j).\end{array}$$

$$I_{\ell }(k_p^{\ell },k_{p+1}^{\ell };r_i,r_j)={\int }_{k_p^{\ell }}^{k_{p+1}^{\ell }}\mathrm{d}k{f}_{\ell }(k,z_i)f_{\ell }(k,z_j),$$

$$C_{\ell }^{\mathrm{obs}}(z_1,z_2)\approx \sum _{p=0}^{N_{\mathrm{k}}-1}\sum _{i=0}^{N_{{\mathrm{z}}_1}}\sum _{j=0}^{N_{{\mathrm{z}}_2}}{w}_i{w}_{j}W(z_i,z_1)W(z_j,z_2)I_{\ell }(k_p,k_{p+1};z_i,z_j).$$

$$j_{\ell }(x)\approx \sqrt{\frac{\pi }{2\ell +1}}{\delta }^D\left(x-\left(\ell +\frac{1}{2}\right)\right).$$

$$r(z_{\ell }(k))=\frac{l+1/2}{k}=d_{\mathrm{H}}{\int }_0^{z_{\ell }(k)}\frac{\mathrm{d}z}{E(z)}.$$

$$\begin{array}{rl}{C}_{\ell }^{\mathrm{thick}}(z_1,z_2;{\sigma }_1,{\sigma }_2)& \approx \frac{2}{d_{\mathrm{H}}^2(2\ell +1)}{\int }_0^{\infty }\mathrm{d}k{W}_1(z_{\ell }(k);z_1,{\sigma }_1)\\ & \times W_2(z_{\ell }(k);z_2,{\sigma }_2)E^2(z_{\ell }(k))P(k,z_{\ell }(k)).\end{array}$$

$$C(\theta ;z_1,z_2)=\frac{1}{4\pi }\sum _{\ell =0}^{{\ell }_{\max }}(2\ell +1){\tilde{C}}_{\ell }(z_1,z_2)P_{\ell }(\cos \theta ),$$

$${\tilde{C}}_{\ell }(z_1,z_2)=C_{\ell }(z_1,z_2)e^{-\ell (\ell +1)/{\ell }_a^2}.$$

$$I={\int }_{-1}^1\mathrm{d}xf(x)g(x).$$

$$I\approx \sum _{k=0}^{N_{\mathrm{cc}}}{w}_{k}f(x_k)g_{(}{x}_k)=\sum _{k=0}^{N_{\mathrm{cc}}}{w}_{k}h(x_k),$$

$$f_N(x)=\frac{a_0}{2}+\sum _{k=1}^{N-1}{a}_k{T}_k(x).$$

$${\mathbf{f}}^{(N)}=(f(t_{\mu }^{(N)}){)}^T\mathrm{with}{t}_{\mu }^{(N)}\equiv \cos \frac{\mu \pi }{N};\mu =0,\dots ,N,$$

$${\mathbf{a}}^{(N)}=\frac{2}{N}{\mathbf{C}}_N^I{\mathbf{f}}^{(N)},$$

$${\mathbf{b}}^{(M)}=\frac{2}{M}{\mathbf{C}}_M^I{\mathbf{g}}^{(M)}.$$

$${\mathbf{h}}^{(P)}=(h(t_{\mu }^{(P)}){)}^T;\mu =0,\dots ,P.$$

$${\mathbf{h}}^{(P)}={\mathbf{f}}^{(P)}\odot {\mathbf{g}}^{(P)}.$$

$${\tilde{\mathbf{a}}}^{(P)}$$

$${\tilde{\mathbf{b}}}^{(P)}$$

$${\mathbf{h}}^{(N_{\mathrm{cc}})}=({\mathbf{C}}_{N_{\mathrm{cc}}}^I{\tilde{\mathbf{a}}}^{(N_{\mathrm{cc}})})\odot ({\mathbf{C}}_{N_{\mathrm{cc}}}^I{\tilde{\mathbf{b}}}^{(N_{\mathrm{cc}})}),$$

Full text

*Methods. *The Angpow software aims at quickly and accurately computing the auto ($z_{1}=z_{2}$) and cross ($z_{1}\neq z_{2}$) angular power spectra between redshift bins. We describe the developed algorithm based on developments on the Chebyshev polynomial basis and on the Clenshaw-Curtis quadrature method. We validate the results with other codes, and benchmark the performance.

*Results. *Angpow is flexible and can handle any user-defined power spectra, transfer functions, and redshift selection windows. The code is fast enough to be embedded inside programs exploring large cosmological parameter spaces through the $C_{\ell}(z_{1},z_{2})$ comparison with data. We emphasize that the Limber’s approximation, often used to speed up the computation, gives incorrect $C_{\ell}$ values for cross-correlations. The C++ code is available from https://gitlab.in2p3.fr/campagne/AngPow.

1 Introduction

Cosmology is entering the era of wide and deep surveys of galaxies, such as, for example, with the Dark Energy Spectroscopic Instrument (DESI) (Levi et al., 2013), the Large Synoptic Survey Telescope (LSST) (Ivezic et al., 2008), and the Euclid satellite (Laureijs et al., 2011), in order to investigate the mechanisms of cosmic acceleration (for a review see Weinberg et al., 2013). Cosmological models can be tested, that is, compared against actual measurements, by studying the statistical properties of galaxy clustering. Several methods exist for this, the most classical ones computing correlations in real (Landy & Szalay, 1993) or Fourier space (Feldman et al., 1994). However, for wider and deeper surveys, one may also try to condense the clustering information into redshift bins (”shells”) and compute the auto- and cross-correlations between redshift shells (Asorey et al., 2012). This is known as tomography, and allows for a more precise understanding of potential systematic errors in different redshift regions. Several studies compare the merits of this kind of approach with the more classical ones (Asorey et al., 2012, 2014; Di Dio et al., 2014; Nicola et al., 2014; Lanusse et al., 2015) and try to optimize the binning to keep most of the cosmological information. All previous studies have been based on the Fisher formalism, which considers observables as Gaussian; unrepresentative, however, of real life data.

To go on further and prepare the future tomographic analyses, one needs to implement a full pipeline and test the method with, for instance, a Monte-Carlo Markov Chain (MCMC) exploration of the cosmological parameter space. But there is a technical bottleneck; running a typical MCMC algorithm in cosmology is already very lengthy and requires computing typically a few $10^{5}$ models. Each model is the result of a numerical code that solves the cosmological equations (known as ”Boltzmann solvers”), such as CLASS111http://class-code.net (Blas et al., 2011), which takes typically 5-10 s on eight cores. Today this amounts to several days of computation.

For a tomographic method, one further needs to transform the output of the Boltzmann solver, the matter power spectrum, into the observable space, represented as $C_{\ell}(z_{i},z_{j})$ angular power spectra between two redshift shells located at $z_{i}$ and $z_{j}$ . This transformation is numerically challenging because of overlapping integrals between very oscillating spherical Bessel functions.

One then often makes use of the Limber’s approximation, which essentially replaces the Bessel functions by a single Dirac value. However, as we show here, this leads to a poor approximation for auto-correlations and may even be incorrect for cross-correlations, since it cannot capture any anti-correlation.

We therefore address here the issue of accurately and quickly computing the integrals required to derive the correlations between tomographic bins. Our goal, in computational terms, is that this computation be faster than one typical Boltzmann code computation time, that is, essentially at or below the 1s level (on eight cores). Another aspect of this work is to provide a stand-alone library that offers a generic interface where the user can plug any matter power spectrum. This is a different approach from CLASSgal (Di Dio et al., 2013) which also provides some theoretical computations related to tomography that are deeply rooted within the CLASS software.

The integrals defining the $C_{\ell}(z_{i},z_{j})$ angular power spectra are introduced in section 2. In section 3 we detail the algorithm implemented in Angpow, while we address some numerical tests in section 4. Section 5 provides insight into the code design and we conclude in section 6.

2 The position of the problem

Our aim is to compute the angular over density power spectrum $C_{\ell}(z_{1},z_{2})$ as a cross-correlation between two $z$-shells with mean values ($z_{1},z_{2}$) and also the auto-correlation $C_{\ell}(z_{1})$ with $z_{1}=z_{2}$, taking into account, in both cases, possible redshift selection functions. Following notations of reference (Di Dio et al., 2013), for a couple of redshift $(z_{1},z_{2})$ one computes $C_{\ell}(z_{1},z_{2})$ according to

[TABLE]

with on one hand $P(k)$ the non-normalized primordial power spectrum, and on the other hand, $\Delta_{\ell}(z,k),$ a general function used to describe physical processes down to redshift $z$ (Di Dio et al., 2013; Bonvin & Durrer, 2011). At the lowest order, $\Delta_{\ell}(z,k)$ can be expressed as the product of the bias $b$ and a growth factor $D(z,k)$ to account for matter density contribution:

[TABLE]

where $j_{\ell}(x)$ is a first kind spherical Bessel function of parameter $\ell$, and $r(z)$ is the radial comoving distance of the shell located at redshift $z$.

For thick redshift shells, one introduces two normalized redshift selection functions $W_{1}(z;z_{1},\sigma_{1})$ and $W_{2}(z^{\prime};z_{2},\sigma_{2})$ around $z_{1}$ and $z_{2}$ with typical width $\sigma_{1}$ and $\sigma_{2}$, respectively. Then, one extends equation 1 by

[TABLE]

It is convenient to introduce the auxiliary function justified in the following section

[TABLE]

where we have used the factorization of the growth factor $D(k,z)$ from the matter density contribution to introduce the notation $P(k,z)$ and $\widetilde{\Delta}_{\ell}(z,k)$. The dots signify that other contributions may be introduced as the redshift distortions and lensing effects that we ignore here for simplicity. Then,

[TABLE]

The auto-correlation is a particular case where the redshift selection function $W_{2}(z^{\prime};z_{2},\sigma_{2})$ is reduced to $W_{1}(z^{\prime};z_{1},\sigma_{1})$ and we can use a single $W$ function, which leads to

[TABLE]

To account for infinite redshift precision at $z=z_{1}$, the use of a Dirac selection function for $W$ yields

[TABLE]

Angpow is designed to efficiently compute these $C_{\ell}$ expressions once one provides access to the power spectra $P(k,z)$, the $\tilde{\Delta}_{\ell}(z,k)$ extra function, and the cosmological distance $r(z)$. To simplify the notation, we do not write the width $\sigma$ of the selection functions if not explicitly needed.

3 A brief description of the computational algorithm

The redshift integral computations of Eq. 5 can be conducted in practice inside the rectangle $[z_{1\mathrm{min}},z_{1\mathrm{max}}]\times[z_{2\mathrm{min}},z_{2\mathrm{max}}]$ given by the $W$ selection functions using a Cartesian product of a one-dimensional (1D) quadrature defined by the set of sample nodes $z_{i}$ and weights $w_{i}$. In practice, we use the Clenshaw-Curtis quadrature. The corresponding sampling points $(z_{1i},z_{2j})$ are weighted by the product $w_{i}w_{j}$ using the 1D quadrature sample points and weights on both redshift regions with $i=0,\dots,N_{\mathrm{z}_{1}}-1$ and $j=0,\dots,N_{\mathrm{z}_{2}}-1$. Then, one gets the following approximation:

[TABLE]

with the notations $z_{i}=z_{1i}$, $z_{j}=z_{2j}$ and $r_{i}=r(z_{1i})$, $r_{j}=r(z_{2j})$ and

[TABLE]

defined with the $f_{\ell}(z,k)$ function of equation 4.

We use a piecewise integration over a user-defined range $[k_{\mathrm{min}},k_{\mathrm{max}}]$ such that

[TABLE]

where the $k^{\ell}_{p}$ bounds are related to the roots of $j_{\ell}(x)$ noted $q_{\ell p}$ and the user-defined number of roots $q_{\ell p}$ per sub-interval $[k^{\ell}_{p},k^{\ell}_{p+1}]$ (see appendix A). Then, equation 8 may be rewritten as

[TABLE]

The integral $I_{\ell}(k^{\ell}_{p},k^{\ell}_{p+1};r_{i},r_{j})$ defined as

[TABLE]

is computed using the 3C-algorithm of appendix A. Investigating the different steps of the algorithm, one notices that the sampling of the function $f_{\ell}(k,z_{i})$ along the $k$ axis for a given $[k^{\ell}_{p},k^{\ell}_{p+1}]$ interval depends on $z_{i}$ but is independent from $z_{j}$ and vice versa for the $f_{\ell}(k,z_{j})$ function (those samplings are independent from $z_{i}$). So, one may proceed to $k$-sampling before performing the double sum over $(i,j),$ which is particularly efficient as the CPU bottleneck is the computation of the spherical Bessel function $j_{\ell}$. Angpow uses a spherical Bessel function implementation based on Numerical Recipes (Press et al., 1992). The brute force Cartesian quadrature evolves as $O(N_{\mathrm{z}_{1}}\times N_{\mathrm{z}_{2}}),$ while the optimized version reduces the CPU times scaling to $O(N_{\mathrm{z}_{1}}+N_{\mathrm{z}_{2}})$. Therefore, as a matter of efficiency, it is more appropriate to exchange the order of the $p$ and $(i,j)$ summations leading to

[TABLE]

As a first hint for the 3C-algorithm, we use typically $2^{8}-2^{9}$ polynomial approximations for $\ell_{\mathrm{max}}\approx 500-1000$ of the $f_{\ell}(k,r_{i})$ and $f_{\ell}(k,r_{j})$ functions, $100$ spherical Bessel roots per sub-interval, $99$-point Clenshaw-Curtis quadrature nodes, and weights for $z$ integration in case of $\sigma=0.02$.

4 Numerical tests

4.1 Comparison with other codes

We proceed now to a numerical comparison of the estimation of $C_{\ell}(z_{1})$ and $C_{\ell}(z_{1},z_{2})$ computed by CLASSgal (Di Dio et al., 2013) and Mathematica (Wolfram Research Inc., 2016) with Dirac redshift selection functions. Concerning CLASSgal (v1.1.3), we have started with the provided explanatory.ini file where we have modified the cosmological parameters to: $h=0.679$, $\Omega_{\mathrm{b}}=0.0483$, $\Omega_{\mathrm{cdm}}=0.2582$ and $\Omega_{\mathrm{k}}=\Omega_{\mathrm{fld}}=0$. We have also set k_scalar_max_tau0_over_l_max to fix the upper bound of the $k$-integration taking into account the maximal $\ell$ value and the redshift mean value. Concerning Angpow we have taken advantage of the possibility to read an external file produced by CLASSgal as an input $P(k)$ computed at $z=0$ in association to the growth function computed internally given by (Lahav et al., 1991; Carroll et al., 1992). Finally, to avoid the Limber’s approximation for CLASSgal, we have set the ”Limber” threshold much higher than the $\ell$ upper limit.

Results of the auto-correlation $C_{\ell}^{\delta}(z)$ computations at $z=1$ using Dirac selection functions are shown in Fig. 1. As the three softwares use the same primordial power spectrum, all the results agree with one other within a maximal relative error of $0.06\%$ on the whole $\ell$ range.

The $C_{\ell}^{\mathrm{thick}}(z)$ auto-correlation computation within a thick redshift band, selected by a Gaussian of mean $z=1$ and $\sigma=0.03$ cut at $\pm 5\sigma$, has been used as a test bench. But, for this test we have neglected the Mathematica software, which is too slow, and have restricted testing to comparison of Angpow to CLASSgal. Figure 2 shows computation results. The orange, violet, and red curves are produced by Angpow using Dirac selection function in the range $\pm 5\sigma$ around the mean redshift $z=1$, while the green, forest green, purple, and blue curves are results of the above mentioned Gaussian selection function but sampled using different grid sizes: $9\times 9$, $19\times 19$, $39\times 39$ and $159\times 159$ Clenshaw-Curtis sample points. As the number of points, or equivalently the quadrature order, increases, the $C_{\ell}^{\mathrm{thick}}(z)$ computation converges towards the CLASSgal result (cyan points). We also address the cross-correlation computations performed by Angpow and CLASSgal ($C_{\ell}^{\mathrm{thick}}(z_{1},z_{2})$) using Gaussian selection functions ($\sigma=0.01$ and a $\pm 5\sigma$ cut). The results are shown in Figure 3. One should not be surprised by negative values since we are cross-correlating two different quantities. In both tests we have used the power spectrum computed at $z=0$ by CLASSgal as input to Angpow. The agreement between the two software codes is very good, keeping the relative residuals at values less than $0.02\%$.

4.2 Note on Limber’s approximation

The Angpow library can also be used to compute, if desired, the first order Limber’s approximation (Loverde & Afshordi, 2008). In such an approximation, the spherical Bessel function is formally reduced to

[TABLE]

In such conditions, the product $kr(z)$ is constrained, and if one uses the following notation for the comoving distance computation with $d_{\mathrm{H}}=c/H_{0}$ , the Hubble distance ($H_{0}=100h\,\mathrm{km}.\mathrm{s}^{-1}\mathrm{Mpc}^{-1}$ and $c$ the speed of light) and $E(z),$ the dimensionless Hubble parameter,

[TABLE]

Then, Eq. 5 is transformed to the following expression

[TABLE]

This integral can be computed using a divide-and-conquer recursive method with the Gauss-Kronrod quadrature (Laurie, 1997). The Gauss sample points are a subset of the Gauss-Kronrod sample points and can easily be used to set up an error estimate to drive the recursive algorithm.

Looking at Eq. 16 one realizes that all the terms are positive, indicating that such approximation is not suitable for cross-correlation where $C_{\ell}(z_{1},z_{2})$ is not guaranteed to be positive as can be explicitly seen in Figure 3. We have proceeded to the computation of $C_{\ell}(z)$ in the case of a Gaussian selection function of mean $z=1$ and $\sigma=0.03$ for both Angpow and CLASSgal, with/without the Limber’s approximation. The results are shown in Figure 4. The two software codes agree very well and show that the Limber’s approximation can give sizeable errors compared to the exact computation; of the order of the cosmic variance in the given example. So, this Limber’s approximation, even if it runs 100 times faster than the exact computation, should then be used with extreme care not only for cross-correlation but also for auto-correlation.

4.3 Correlations in real space

Angpow can also quickly compute the correlation function in real space from the power spectrum

[TABLE]

where $P_{\ell}$ denotes the $\ell-\mathrm{th}$ order Legendre polynomial. Because the $C_{\ell}(z_{1},z_{2})$ values are generally cut at a given $\ell_{\mathrm{max}}$, one needs to introduce an apodization to avoid ringing due to a sharp cut-off. We introduce a Gaussian one (which is the smoothest in both real and harmonic spaces) as in (Di Dio et al., 2013) so that the $\tilde{C}_{\ell}(z_{1},z_{2})$ term in Eq. 17 is

[TABLE]

The apodization length $\ell_{a}$ may depend on the signal but for the standard cosmology shown here (around $z=1$) we noticed that using $l_{a}\simeq 0.4\ell_{max}$ gives good results. Correlations in real space are generally easier to comprehend as is shownin Figure 5, which represents the analogue of Figure 3 but in real space. Here one may identify a peak, named the ”Baryonic Acoustic Oscillation” (e.g., Weinberg et al., 2013) that decreases in the cross-correlation when the distance between shells increases and is finally washed out.

4.4 Speed tests

Angpow is designed and written in C++ and parallelization is achieved through OpenMP. To qualify the code we provide four input parameter files and their corresponding outputs obtained in one of our runs. In all tests, we used a $\Lambda$CDM reference cosmology, a $P(k)$ at $z=0$ provided by CLASSgal, $\ell_{\mathrm{max}}=1000$, a 3C-algorithm with $2^{9}$ Chebyshev polynomial order, and 100 roots per sub-$k$-interval:

Test 1

: Auto-correlation with a Dirac selection function at $z=1$ and $k_{\mathrm{max}}=10\ \mathrm{Mpc}^{-1}$;

Test 2

: Cross-correlation with two Dirac selection functions at $z=1$ and $z=1.05$ and $k_{\mathrm{max}}=10\ \mathrm{Mpc}^{-1}$;

Test 3

: Auto-correlation with a Gaussian selection function with ($z_{\mathrm{mean}}=1$, $\sigma_{z}=0.02$, $5\sigma_{z}$-cut) and $k_{\mathrm{max}}=1\ \mathrm{Mpc}^{-1}$ and a radial quadrature based on $N_{\mathrm{pts}}=139$ sample points;

Test 4

: Cross-correlation with two Gaussian selection functions with ($z_{\mathrm{mean,1}}=0.90$, $z_{\mathrm{mean,2}}=1.00$) both with ($\sigma_{z}=0.02$, $5\sigma_{z}$-cut) and $k_{\mathrm{max}}=1\ \mathrm{Mpc}^{-1}$ and a radial quadrature based on $N_{\mathrm{pts}}=139$ sample points;

We have tested the code both on laptop (Linux, MacOSX) as well as on Computer Center (CCIN2P3 in France and NERSC in the USA). We use OpenMP to distribute the computation of a given $C_{\ell}$ on a single thread. Table 1 gives Central Processing Unit (CPU) execution times averaged over ten processes. The code wall time decreases reasonably well with the number of threads and a wall time at the $\mathcal{O}$(1s) level can be reached to reconstruct these accurate spectra. Such performances are much higher than those obtained with CLASSgal when not using the Limber approximation. For example, on our Test-3 setup, running the latter takes about 15s (on 16 threads), which is to be compared to about 0.5s in Table 1.

Figure 6 shows the dependence of the wall time with respect to the width of the selection function ($\sigma$) in the conditions of Test 3 using 16 threads. For a given $\sigma$ , we have chosen the minimal radial_order value such that the relative accuracy on the $C_{\ell}$ is of the order of 0.01% compared to a computation with a much larger radial_order value (see section 5). If one uses a looser criteria on the accuracy of the $C_{\ell}$ or if the number of sigma is lower than 5, then one may use a lower radial_order and gain on the wall time.

5 Code design and input parameters

Angpow is written in C++ which allows both good CPU performances and keeps the code flexible. A front end interface to Python is also foreseen and the code is distributed publicly at https://gitlab.in2p3.fr/campagne/AngPow. The angpow.cc file is an example of the library usage to perform $C_{\ell}$ and $C(\theta)$ computations. We also provide the limber.cc file if one wants to test the Limber’s approximation (Sect. 4.2). The different files are located in self-explained directories: src, inc/Angpow, lib, data. Finally, a README file provides details for the installation and build procedures.

The two main classes Pk2Cl and KIntegrator (located in angpow_pk2cl.h and angpow_kinteg.h files) are generic codes using abstract base classes. They define interface to the power spectrum function $P(\ell,k,z)$ (class PowerSpecBase), the generalisation of $P(k,z)$ used in Eq.  5; to the comoving distance computation $r(z)$ (class CosmoCoorBase); and to the radial/redshift selection functions $W(z)$ (class RadSelectBase). The user can derive their own concrete classes to access a third party library or use the ones implemented by default. For instance, we have implemented a file access to a $(k,P(k))$-tuple saved by the CLASSgal output. In this implementation, we have coded the growth factor defined in (Lahav et al., 1991; Carroll et al., 1992) as minimal $\tilde{\Delta}_{\ell}(z,k)$ function (Eq. 2) .

To run the executable, one provides an ascii file defining the input parameters that drive the computational conditions of the algorithm and define the I/O locations. Some ready-made input parameter files are also provided (angpow_bench<n>.ini) as well as the $C_{\ell}$ output files (angpow_bench<n>_cl.txt.REF) corresponding to the icpc outputs of Table 1.

Among the different input parameters some are more sensitive than others, as those that deal with the radial/redshift 1D quadrature and the 3C algorithm. Here is a closer look at these parameters:

• •

cl_kmax: This is the maximal value of $k$ in the $k$-integral. We have not set up an internal algorithm to determine this upper bound. As a hint, one may consider a relation with the factor $\ell_{\mathrm{max}}/r(z_{\mathrm{min}})$. The lower bound on $k$ is internally fixed using the cut-off $x_{\mathrm{min}}(\ell)$ defined as $x<x_{\mathrm{min}}(\ell)\Leftrightarrow j_{\ell}(x)<\mathrm{cut}$ (see the input parameters jl_xmin_cut and Lmax_for_xmin set as default to $5\ 10^{-10}$ and 2000, respectively).

• •

radial_order: If noted $n$, this fixes the number of sample points along one $z$ direction, that is, $N_{\mathrm{pts}}=2n-1$ . The accuracy on the selection function as well as the CPU increase with $n$ but we keep a $O(n)$ complexity of the $k$-sampling ( see Figure 6 );

• •

chebyshev_order: If noted $N$, this fixes the degree $d$ of the Chebyshev polynomial approximation of the $f_{\ell}(k,z_{i})$ and $f_{\ell}(k,z_{j})$ functions (Eq. 4); that is, $d=2^{N}$. Keeping the same degree of approximation for both functions guarantees a power of 2 for the product approximation. Even if not mandatory, this helps in getting CPU performance for the DCT-I transform (using the FFTW library). When $\ell_{\mathrm{max}}$ increases it may be worth updating this parameter by 1 unit step. For $\ell_{\mathrm{max}}=500$ chebyshev_order is set to 8 by default. Increasing the angular spectrum computation up to $\ell_{\mathrm{max}}=1000$ keeping this default value leads to absolute error of the order of $10^{-6}$ for Tests 1 and 2 and $5\ 10^{-10}$ for Test 3 and 4, then to get better accuracy in this case we use $\mathtt{chebyshev\_order}=9$.

• •

n_bessel_roots_per_interval: This is the number of Bessel roots $q_{\ell p}$ used to define the bounds of the integral $I_{\ell}(k^{\ell}_{p},k^{\ell}_{p+1};r_{i},r_{j})$ (Eq. 12). By default it is set to 100. There is an interplay with the chebyshev_order parameter as a lower n_bessel_roots_per_interval value is coherent with a lower chebyshev_order.

• •

total_weight_cut and deltaR_cut: These two threshold parameters are used to avoid unnecessary 3C algorithm processing (especially for the $k$-sampling of the $f_{\ell}(k,z_{i})$ or $f_{\ell}(k,z_{j})$ functions). So, we do not consider a couple $(z_{i},z_{j})$ for which either the product $w_{i}w_{j}W(z_{i},z_{1})W(z_{j},z_{2})$ is too low (total_weight_cut cut) or the radial distance $|r(z_{i})-r(z_{j})|$ is too large to produce a sizeable contribution to the final $C_{\ell}$. The deltaR_cut cut is in Mpc units and is used in conjunction with has_deltaR_cut set to 1. These two threshold parameters depend on the use case under consideration and for preliminary tests we recommend to set both total_weight_cut and has_deltaR_cut to 0.

6 Summary and outlooks

We have set up a fast and generic software to compute the tomographic $C_{\ell}(z_{1},z_{2})$ with redshift selection functions. Angpow is versatile enough to accept user-defined matter power spectrum $P(k)$, transfer functions $\tilde{\Delta}_{\ell}(k,z),$ and cosmology. The code provides an accurate computation of the auto and cross-correlation power spectra, checked by comparison with other codes, which is fast enough to be included inside MCMC cosmology softwares. The rapidity of the software relies on the use of the 3C-algorithm, adapted to the computation of integrals over spherical Bessel functions, while other codes rely on the Limber’s approximation. We emphasize that the Limber’s approximation can lead to incorrect $C_{\ell}(z_{1},z_{2}),$ especially in the case of cross-correlations, as Limber’s Dirac functions cannot model the interferences between two spherical Bessel functions at different redshifts.

This code is thus fast and accurate enough to be used to test cosmological parameters, and perform tomographic analysis of the galaxy distribution. The definition of the $\tilde{\Delta}_{\ell}(k,z)$ function is general and can accept zero order function as $\tilde{\Delta}^{\rm mat}_{\ell}(k,z)$, but also relativistic corrections such as the redshift space distortions or the gravitational lensing; despite these corrections, it can also contain spherical Bessel functions (see e.g., Di Dio et al., 2013). Because Angpow provides the correct angular power spectra for cross-correlations, it can be a key tool to perform an integrated approach to cosmology, as advertised in (Nicola et al., 2016). In particular, we propose this tool for deriving the auto- and cross-correlation angular power spectra for galaxy clustering, but also with angular power spectra from cosmic shear and the cosmological microwave background.

Angpow is publicly available222from https://gitlab.in2p3.fr/campagne/AngPow. and can be interfaced to other codes; a Python interface is foreseen. At the moment the code only accepts two redshift bins but soon it will be generalized to any number of bins. Feedback from Angpow users would be greatly accepted.

Acknowledgements.

The authors want to thank M. Reinecke who kindly provided pieces of the code, and J. D. McEwen for fruitful discussion on the Chebyshev transform.

Appendix A Clenshaw-Curtis-Chebyshev algorithm (3C-algorithm)

Each integral type of Eq. 12 involves the product of ”highly” oscillatory functions. The purpose of this section is not to provide a review of all the integration methods used in the different fields of physics to tackle such a difficult task.. To focus on our case, where we have to deal with (at least) the product of spherical Bessel functions, the authors point out that the reader may find either specific integral solving rules as in (Glasser & Montaldi 1993) or general methods as in (Lucas & Stone 1995). However we need a precise and also very fast method. We cannot rely on methods that solve the problem of a product of spherical Bessel functions multiplied by a regular function. In fact, both the primordial power spectrum and the extension beyond the matter density $\tilde{\Delta}^{\mathrm{mat.}}(z,k)$ may show oscillation features in the form of derivative of spherical Bessel functions. So, we have searched and found a general method that meets our requirements of precision and speed.

Eq. 12 is a special case of the following generic integral after a proper change of variable

[TABLE]

To get an approximate value of this integral, we use in this section the Clenshaw-Curtis quadrature at order $N_{\mathrm{cc}}$ (noting $h=f\times g$):

[TABLE]

where the sampling points are defined as $x_{k}=\cos k\pi/N_{\mathrm{cc}}$ ($k=0,\dots,N_{\mathrm{cc}}$) and the corresponding weights $w_{k}$ are addressed later in the section. But, if the functions $f$ and $g$ have a highly oscillatory behavior, one needs, in principle, to use large values of $N_{\mathrm{cc}}$ to reach a sufficient accuracy level. In that case, dealing with the above sum may not be computationally efficient. The idea is then to use Chebyshev series to approximate both functions $f$ and $g$. Then, one performs the product of both Chebyshev series, which is also a Chebyshev series but with a higher order, and finally one uses a fast Clenshaw-Curtis weights computation to perform the last weighted sum. We briefly describe those steps leaving the details of the demonstration that the interested reader can find in (Baszenski & Tasche 1997).

Let $f_{N}$ be a polynomial approximation of $f$ of degree $N-1$. We expend $f_{N}$ onto the following basis of the first kind of Chebyshev polynomials $\{T_{n};n=0,\dots,N-1\}$ which have the property $T_{n}(\cos\theta)=\cos n\theta$:

[TABLE]

To determine the $a_{k}$ coefficients one uses the following sampling vector

[TABLE]

of length $N+1$ and related to the vector $\mathbf{a}^{(N)}=(a_{0},\dots,a_{N-1},0)^{T}$ by the linear algebra relation

[TABLE]

with $\mathbf{C}^{I}_{N}$ being a discrete cosine transform of type-I (DCT-I) matrix of dimension $(N+1)^{2}$. Similarly, we note $g_{M}$ a polynomial approximation of $g$ of degree $M-1$ from which we determine the sampling vector $\mathbf{g}^{(M)}$ using the sample points $t_{\mu}^{(M)}$. The coefficient vector $\mathbf{b}^{(M)}=(b_{0},\dots,b_{M-1},0)^{T}$ is related to $\mathbf{g}^{(M)}$ using a relation similar to Eq. 23, namely

[TABLE]

By combining the polynomial approximations $f_{N}$ and $g_{M}$, the function $h$ is then approximated by a Chebyshev series of degree $N+M-2$ with coefficient vector $\mathbf{c}^{(P)}$ of length $P+1$ with $P\geq N+M-1$. Using a relation of type Eq. 23, the vector $\mathbf{c}^{(P)}$ is related to the sampling vector

[TABLE]

To get $\mathbf{h}^{(P)}$ it is not necessary to compute $\mathbf{c}^{(P)}$ and proceed to an inversion of a DCT-I matrix. The key point is that if we note $\odot,$ the component-wise multiplication, one has

[TABLE]

Moreover, $\mathbf{f}^{(P)}$ ($\mathbf{g}^{(P)}$) is obtained from $\mathbf{a}^{(N)}$ ($\mathbf{b}^{(M)}$) of length $N+1$ ($M+1$) by an extension to a larger vector at least of length $N+M$ noted $\tilde{\mathbf{a}}^{(P)}$ ($\tilde{\mathbf{b}}^{(P)}$) by appending with zeros:

[TABLE]

Then, the sampling vector used in Eq. 20 where one identifies $N_{\mathrm{cc}}=P$ is determined by

[TABLE]

using $\mathbf{C}^{I}_{N_{\mathrm{cc}}}$ the DCT-I matrix of dimension $(N_{\mathrm{cc}}+1)^{2}$ and the inversion property $(\mathbf{C}^{I}_{P})^{-1}=(2/P)\mathbf{C}^{I}_{P}$. In some sense, for both $f$ and $g$ approximation sampling vectors, we have performed a Chebyshev basis change to a larger parameter space compatible with the polynomial degree involved in the product $f^{(N)}\times g^{(M)}$.

The second key point is that the Clenshaw-Curtis weights associated to $\mathbf{h}^{(N_{\mathrm{cc}})}$ in Eq. 20 can also be computed with a DCT-I transform from the vector $(2/N_{\mathrm{cc}})(1,0,-1/3,0,-1/15,\dots,((-1)^{k}+1)/2(1-k^{2}),\dots)$ of length $N_{cc}+1$ (Waldvogel 2006) (the normalization depends on the exact definition of the DCT-I used).

So, to perform the integral given by Eq. 20, one needs $4+1$ DCT-I transforms, 1 for the Clenshaw-Curtis weights and 4 to transform the Chebyshev coefficients vectors. The DCT-I transform may be implemented using $O(n\log n)$ efficient algorithm, for example, the FFTW library (Frigo & Johnson 2005) used by Angpow. Angpow uses a power of $2$ for $N$, $M,$ and also $P$ (keeping $P\geq N+M-1$) for a fast implementation. The 3C-algorithm is a special case of a more general class of algorithms dealing with the product of polynomials. We note that according to reference (Giorgi 2012) an even faster algorithm (although moderate) might be implemented in a future version of Angpow if necessary. We note finally that this general method can be applied to use cases beyond the power spectrum computation in other fields of interest.