Free Access
Issue
A&A
Volume 540, April 2012
Article Number A60
Number of page(s) 8
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/201118463
Published online 27 March 2012

© ESO, 2012

1. Introduction

In the last few decades, cosmology has become a data-driven field, where high-precision measurements of the cosmic microwave background (CMB, e.g., Larson et al. 2011), weak lensing (e.g., Schrabback et al. 2010) and galaxy surveys (e.g., Percival et al. 2007b) have permitted the establishment of a standard cosmological model in which the Universe is composed of 4% baryons, 22% dark matter and 74% dark energy. Some major questions remain, the nature of dark matter and dark energy in particular is still not understood. Similarly, the initial conditions of the Universe are yet to be established and alternative models of gravity are still to be tested in comparison with Einstein’s general relativity.

New surveys are underway with these science objectives, e.g. Planck for the CMB (Planck Collaboration 2006), DES (Dark Energy Survey, Annis et al. 2005), BOSS (Baryon Oscillation Spectroscopic Survey, Schlegel et al. 2007), LSST (Large Synoptic Survey Telescope, Tyson & LSST 2004) and Euclid (Laureijs et al. 2011; Refregier et al. 2011) for weak lensing and the study of large-scale structure with galaxy surveys. In order to be beneficial, cosmological studies of these surveys need to use high-precision statistical methods, such as a full 3D analysis on the sky where all-sky 3D surveys are available.

Several tools have been developed to analyse data on the sphere, which is required for a 2D spherical harmonic CMB analysis (Crittenden & Turok 1998; Crittenden 2000; Górski et al. 2002, 2005; Doroshkevich et al. 2005). Weak lensing and galaxy survey data can also be analysed tomographically (i.e. in 2D slices), but unlike for the CMB, a full 3D spherical Fourier-Bessel analysis can also be sought (Fisher et al. 1995; Heavens & Taylor 1995; Heavens 2003; Castro et al. 2005; Rassat & Refregier 2012). Previous 3D data analyses were on relatively small data sets (Fisher et al. 1995; Heavens & Taylor 1995; Erdoğdu et al. 2006b,a), but future surveys like Euclid and LSST will provide surveys with billions of galaxies, making previous methods for calculating the 3D spectra unfeasibly time-consuming.

In Sect. 2.1, we present the theory behind the 3D Fourier-Bessel decomposition for infinite and finite continuous fields as well as the usual method for a discrete survey (e.g. galaxy survey). In Sect. 2.2, we present two additional equivalent formulations of the spherical Fourier-Bessel decomposition, one of which is central to the 3DEX code. In Sect. 3, we compare the accuracy and calculation time for the usual method used for calculating Fourier-Bessel coefficients and methods with the 3DEX code presented in this paper. In Sect. 4, we describe the 3DEX library and give examples of how to use it. In Sect. 5 we present our conclusions. We also include an appendix, where we discuss the subtleties of the Fourier-Bessel normalisation.

2. Theory

2.1. The spherical Fourier-Bessel decomposition

In observational cosmology, spherical coordinates (where the observer is at the origin) are a natural choice for the analysis of cosmological fields. In this system of coordinates, eigenfunctions of the Laplacian operator are products of spherical Bessel functions and spherical harmonics, i.e. functions j(kr)Yℓm with eigenvalues  − k2. For an homogeneous three-dimensional field f(r) = f(r,θ,φ) in a flat geometry, the spherical Fourier-Bessel decomposition (Fisher et al. 1995; Heavens 2003; Castro et al. 2005) is f(r,θ,φ)=2πdkℓmfℓm(k)kj(kr)Yℓm(θ,φ),\begin{eqnarray} \hspace*{-0mm} f(r,\theta,\phi) = \sqrt{\frac{2}{\pi}}\int {\rm d}k \sum_{\ell m} f_{\ell m}(k) k j_\ell (kr) Y_{\ell m}(\theta,\phi) ,\label{thraw1} \end{eqnarray}(1)with the inverse relation fℓm(k)=2πd3rf(r,θ,φ)kj(kr)Yℓm(θ,φ).\begin{eqnarray} f_{\ell m}(k) = \sqrt{\frac{2}{\pi}} \int {\rm d}^3\vec{r}\ f(r,\theta,\phi)k j_\ell (kr)Y^{*}_{\ell m}(\theta,\phi) .\label{thraw2} \end{eqnarray}(2)Note that this decomposition uses the same notation as Rassat & Refregier (2012) and Castro et al. (2005), which is slightly different from the one used in Lanusse et al. (2012). The coefficients may be used to calculate the 3D power spectrum C(ℓ,k), defined by fℓm(k)fm(k)=C(l,k)δD(kk)δδmm,\begin{equation} \left< f_{\ell m}(k)f^*_{\ell' m'}(k')\right>=C(l,k) \delta_D(k-k')\delta_{\ell \ell'}\delta_{mm'}, \end{equation}(3)a naïve estimator of which is C(k)=12l+1m|fℓm(k)|2.\begin{eqnarray} C_\ell (k) = \frac{1}{2l+1} \sum_m |f_{\ell m}(k)|^2. \end{eqnarray}(4)This can be seen as an extension of the usual 2D power spectrum fℓmfm=Clδδmm\hbox{$\left< f_{\ell m}f^*_{\ell' m'}\right>=C_l \delta_{\ell \ell'}\delta_{mm'}$}. The latter arises from the spherical harmonic transform of a 2D field given on the sphere f(θ,φ) =  ∑ ℓmfℓmYℓm(θ,φ).

In practice, surveys will only cover a finite amount of volume, limiting the analysis to a sphere of radius R. These boundary conditions lead to a discrete spectrum  { kℓn } , which is detailed in the appendices. In this paper, we assumed as a boundary condition that f vanishes at r = R. The spherical Fourier-Bessel decomposition becomes (Erdoğdu et al. 2006b; Fisher et al. 1995) f(r,θ,φ)=ℓmnκℓnfℓm(kℓn)kℓnj(kℓnr)Yℓm(θ,φ),\begin{eqnarray} \hspace*{-0mm}f(r,\theta,\phi) = \sum_{\ell mn} \kappa_{\ell n} f_{\ell m}(k_{\ell n}) k_{\ell n} j_\ell (k_{\ell n}r)Y_{\ell m}(\theta,\phi), \label{fbreconstr} \end{eqnarray}(5)which is exact if the ranges of ,m and n are infinite. The Fourier-Bessel coefficients are denoted by fℓmn = fℓm(kℓn), and κℓn is the normalisation constant (see Appendices for more details).

In various applications, though, the continuous field f cannot be directly observed. This is notably the case in cosmology where galaxy surveys give indirect information about the underlying matter density field through their spacial positions. Note that these tracers are subject to various distortions and non-linearities, but these are not the purpose of this work. In this work we only consider linear or quasi-linear scales ( < 50, k < 0.2   h   Mpc-1).

If the only information about the field f is a list of coordinates rp = (rp,θp,φp) with p = 1,...,N (where N is the number of galaxies in the latter example), the survey may be considered as a superposition of 3D Dirac deltas and each coefficient fℓmn can simply be estimated with a sum (Heavens & Taylor 1995; Fisher et al. 1995; Erdoğdu et al. 2006b; Abramo et al. 2010) ˜f(r)=p=1Nδ(3)(rrp),˜fℓmn=\begin{eqnarray} \tilde{f}(\vec{r}) & = & \sum^N_{p=1} \delta^{(3)}(\vec{r}-\vec{r}_p), \\ \tilde{f}_{\ell mn} & = & \sum^N_{p=1} k_{\ell n}j_\ell (k_{\ell n}r_p)Y^{*}_{\ell m}(\theta_p,\phi_p). \label{rawestimate} \end{eqnarray}This formulation has been used for the analysis of shallow galaxy surveys such as the IRAS 1.2 mJ survey (~6   k galaxies, Strauss et al. 1992; Fisher et al. 1995; Heavens & Taylor 1995), and the 2MRS survey (2MASS Redshift Survey,  ~45 k galaxies, Huchra et al. 2012; Erdoğdu et al. 2006b,a). Since the time to calculate Eq. (7) is proportional to Nnmax(max + 1)2/2, Eq. (7) will become highly time-consuming when applied to larger surveys or when precise decomposition is required (large nmax and max).

2.2. Three equivalent formulations

In spherical coordinates, since 3D space can be viewed as an infinite series of closed shells Ω(r), the spherical Fourier-Bessel decomposition may also arise from repeated 2D spherical harmonic transforms to which spherical Bessel transforms are applied (Abramo et al. 2010). Formally, the field f given on each shell Ω(r) is first expanded into spherical harmonics f(r,θ,φ)=ℓmfℓm(r)Yℓm(θ,φ),\begin{eqnarray} f(r,\theta,\phi) = \sum_{\ell m} f_{\ell m}(r) Y_{\ell m}(\theta,\phi) ,\label{thdirect1} \end{eqnarray}(8)for which the inversion formula gives harmonics coefficients fℓm(r) depending on the radius rfℓm(r)=Ω(r)f(r,θ,φ)Yℓm(θ,φ).\begin{eqnarray} f_{\ell m}(r) = \int_{\Omega(r)} {\rm d}\Omega \ f(r,\theta,\phi) Y^{*}_{\ell m}(\theta,\phi). \label{thdirect2} \end{eqnarray}(9)It is then possible to perform a spherical Bessel transform fℓm(r)=2πdkfℓm(k)kj(kr),\begin{eqnarray} f_{\ell m}(r) = \sqrt{\frac{2}{\pi}} \int {\rm d}k \ f_{\ell m}(k) kj_\ell (kr), \label{thdirect3} \end{eqnarray}(10)leading to the final Fourier-Bessel coefficients fℓm(k) fℓm(k)=2πdrr2fℓm(r)kj(kr).\begin{eqnarray} f_{\ell m}(k) = \sqrt{\frac{2}{\pi}} \int {\rm d}r \ r^2 f_{\ell m}(r) kj_\ell (kr) .\label{thdirect4} \end{eqnarray}(11)This formulation hence extends the notion of 2D spherical harmonics to three-dimensional fields.

It is also possible to conceive the reverse approach, i.e. to perform the spherical Bessel transform first and subsequently expand the resulting coefficients into spherical harmonics. Formally, the th order spherical Bessel transform of f (similar to its Hankel transform) is f(r,θ,φ)=2πdkf(k,θ,φ)kj(kr),\begin{eqnarray} f(r,\theta,\phi) = \sqrt{\frac{2}{\pi}} \int {\rm d}k \ f_{\ell }(k,\theta,\phi) kj_\ell (kr) ,\label{threverse1} \end{eqnarray}(12)for which the inversion formula gives f(k,θ,φ)=2πdrr2f(r,θ,φ)kj(kr).\begin{eqnarray} f_\ell (k,\theta,\phi) = \sqrt{\frac{2}{\pi}} \int {\rm d}r \ r^2 f(r,\theta,\phi) kj_\ell (kr) .\label{threverse2} \end{eqnarray}(13)The result is then expanded into spherical harmonics but with an unusual formulation since f(k,θ,φ) and Yℓm(θ,φ) (as well as the basis functions j(kr) and Yℓm(θ,φ)) have now the parameter in common: f(k,θ,φ)=mfℓm(k)Yℓm(θ,φ).\begin{eqnarray} f_\ell (k,\theta,\phi) = \sum_{m} f_{\ell m}(k) Y_{\ell m}(\theta,\phi) .\label{threverse3} \end{eqnarray}(14)Again, using the inversion formula, we obtain the Fourier-Bessel coefficients fℓm(k)=Ωf(k,θ,φ)Yℓm(θ,φ).\begin{eqnarray} f_{\ell m}(k) = \int_{\Omega} {\rm d}\Omega \ f_\ell (k,\theta,\phi) Y^{*}_{\ell m}(\theta,\phi). \label{threverse4} \end{eqnarray}(15)Due to the closed domains of shells Ω(r) and thus the relative independence of angular and radial dimensions, the raw (Eqs. (1) and (2)), the forward (denoted by SHB for spherical-harmonic-Bessel, Eqs. (8) to (11)) and the reverse (denoted by SBH for spherical-Bessel-harmonic, Eqs. (12) to (15)) methods are equivalent formulations of the spherical Fourier-Bessel decomposition of any three-dimensional field f(r,θ,φ). This is summarised in the following schematic description of each method: RAW:f(r)0three-dimensionalintegrationfℓm(k)SHB:f(r)0SHTfℓm(r)0SBTfℓm(k)SBH:f(r)0SBTf(k,θ,φ)0SHTfℓm(k).\begin{eqnarray} \textnormal{RAW: } \ & & f(\vec{r}) \ \xrightarrow[]{\textrm{three-dimensional integration}} \ f_{\ell m}(k) \nonumber \\ \textnormal{SHB: } \ & & f(\vec{r}) \ \xrightarrow[]{\textrm{SHT}} \ \ \ \ f_{\ell m}(r) \ \ \ \ \xrightarrow[]{\textrm{SBT}} \ \ f_{\ell m}(k) \\ \textnormal{SBH \!: } \ & & f(\vec{r}) \ \xrightarrow[]{\textrm{SBT}} \ \ f_{\ell }(k,\theta,\phi) \ \ \xrightarrow[]{\textrm{SHT}} \ \ f_{\ell m}(k) \nonumber. \end{eqnarray}(16)Note that this section is related to the ideal case R = ∞, but all equations can be straightforwardly rewritten for a finite R by replacing k by kln, bounding each integral and adapting normalisation. The formulas arising from this adaptation are used in the next sections.

2.3. Estimating Fourier-Bessel coefficients from a real survey

Although the three approaches described in Sect. 2.2 are theoretically equivalent, their estimates and numerical implementations take different forms.

2.3.1. Forward approach (SHB)

Estimating the fℓmn coefficients using the forward method naturally requires the radial dimension to be discretised. Indeed, the first step is to compute the spherical harmonic transform on a set of shells located at radial values r1,...,rNlayers. In each layer, the coefficients fℓm(ri) are estimated. Although it is possible to perform a raw estimate for the later harmonics transform, it is often advisable to use a robust 2D discretisation scheme (of Npix(i) pixels for the ith shell) and to take advantage of the related high-performance algorithms. Angular space is hence discretised into nodes (ri,θp,φp) = (ri,γq) and the field is approximated on each node, giving ˜f(ri,γp)\hbox{$\tilde{f}(r_i,{\boldsymbol{\gamma}}_p)$}. The spherical harmonic decomposition in the ith shell becomes ˜fℓm(ri)=p=1Npix(i)˜f(ri,γp)Yℓm(γp),\begin{eqnarray} \tilde{f}_{\ell m}(r_i) = \sum^{Npix(i)}_{p=1} \tilde{f}(r_i,{\boldsymbol{\gamma}}_p)Y^{*}_{\ell m}({\boldsymbol{\gamma}}_p), \end{eqnarray}(17)and the final coefficients are obtained by performing the following spherical Bessel decomposition: ˜fℓmn=i=1Nlayers˜fℓm(ri)kℓnj(kℓnri).\begin{eqnarray} \tilde{f}_{\ell mn} = \sum_{i=1}^{N_{\rm layers}}\tilde{f}_{\ell m}(r_i) k_{\ell n}j_\ell (k_{\ell n}r_i). \end{eqnarray}(18)With this method, radial and angular spaces are discretised and both transforms are approximated.

2.3.2. Reverse approach (SBH)

For the reverse approach, a 2D scheme on the sphere was required as well. As previously, this scheme defines a set of Npix zones (pixels) related to angular nodes γq. If Gq denotes the points of the survey located in the solid angle corresponding to the qth zone of the scheme, we perform the spherical Bessel transform (raw estimate) in each zone ˜fℓn(γq)=˜f(kℓn,γq)=pGqkℓnj(kℓnrp),\begin{eqnarray} \tilde{f}_{\ell n}({\boldsymbol{\gamma}}_q) = \tilde{f}_\ell (k_{\ell n},{\boldsymbol{\gamma}}_q) = \sum_{p \in G_q} k_{\ell n}j_\ell (k_{\ell n}r_p), \label{revdiscr1} \end{eqnarray}(19)and each of these intermediate maps is decomposed into spherical harmonic (spherical harmonics transform), which gives the Fourier-Bessel coefficients ˜fℓmn=q=1Npix˜fℓn(γq)Yℓm(γq).\begin{eqnarray} \tilde{f}_{\ell mn} = \sum^{Npix}_{q=1} \tilde{f}_{\ell n}({\boldsymbol{\gamma}}_q)Y^{*}_{\ell m}({\boldsymbol{\gamma}}_q) \label{revdiscr2} . \end{eqnarray}(20)With the reverse method, one can avoid to discretise radial space. Moreover, this one-shell pixelisation of the sky (thus based on physical solid angles) allows for a natural treatment of radial distortions (redshift, relativistic) and masking effects. Using multiple resolutions at different radial values, as would be possible with the forward method, is much more questionable. The SHB method also proves to be a powerful tool for weighting the data prior to estimating the power spectrum. For instance, Tadros et al. (1999) used a fiducial power spectrum to derive an optimal weighting operation. This operation is quite complex when using the raw Fourier-Bessel approach, whereas the SHB formulation naturally handles the dependence on k of the weighting function.

The three methods to estimate the spherical Fourier-Bessel decomposition can therefore also be expressed for a discrete 3D survey, summarised schematically below:

RAW:{rp}0Rawsum,bestestimateofFBcoefficients˜fℓmnSHB:{rp}0ApproxSHT˜fℓm(ri)0ApproxSBT˜fℓmnSBH:{rp}0ExactSBT˜fℓn(γp)0ApproxSHT˜fℓmn.\begin{eqnarray} \textnormal{RAW: } &\{\vec{r}_p\} \xrightarrow[]{\textrm{Raw sum, best estimate of FB coefficients}} \ \tilde{f}_{\ell mn} \nonumber \\ \textnormal{SHB: } & \{\vec{r}_p\} \xrightarrow[]{\textrm{Approx SHT}} \ \tilde{f}_{\ell m}(r_i) \ \ \xrightarrow[]{\textrm{Approx SBT}} \ \ \tilde{f}_{\ell mn}\nonumber \\ \textnormal{SBH: } & \{\vec{r}_p\} \xrightarrow[]{\textrm{Exact SBT}} \ \ \ \tilde{f}_{\ell n}(\boldsymbol{\gamma}_p) \ \ \xrightarrow[]{\textrm{Approx SHT}} \ \tilde{f}_{\ell mn} \nonumber. \end{eqnarray}Note that in practice, the range of (l,m,n) is finite, which introduces an additional approximation. Here, and n are restricted to  [0,max  and  [1,nmax , respectively. Given , m goes from  − to .

3. Method comparison

3.1. Complexity, accuracy and discretisation grids

For a survey that probes a field by N discrete points, the raw method is the natural estimate of the Fourier-Bessel coefficients. However, since each point contributes to the calculation of every coefficient ˜fℓmn\hbox{$\tilde{f}_{\ell mn}$} (∀  l,m,n), computation time is proportional to N·nmax(max + 1)2/2, which can be highly problematic for large surveys.

In the forward method, the repeated spherical harmonic transforms take advantage of tesselation schemes and high-performance algorithms such as those provided by HEALPix Górski et al. (2005), IGLOO Crittenden & Turok (1998) or GLESP Doroshkevich et al. (2005). Roughly speaking, the number of nodes to be considered is reduced from N to Npix, and the use of fast spherical harmonic transforms on these schemes significantly decreases computation time.

However, this approach requires the 3D space to be divided into shells Ω(ri). Both radial and angular dimensions are discretised, and the survey is approximated on an actual 3D grid. In practice, this approximation deteriorates the accuracy of the estimated Fourier-Bessel coefficients. Furthermore, designing a meaningful radial discretisation is a difficult task. For equal-area pixelisations, the area of each pixel on the ith shell is 4πri2/Npix(i)\hbox{$4\pi r_i^2/N_{\rm pix}(i)$}. With HEALPix, the nside angular parameter may only be increased by a factor 2, which changes the number of pixels by a factor of 4 (since Npix(i) = 12nside(i)2). This means that pixel areas cannot be stabilised for subsequent shells as r increases. Consequently, it is not possible to adapt a resolution to obtain a 3D scheme with equal-volume voxels. Extending 2D schemes to 3D is difficult and may even require a novel formalism for an equal-voxel 3D grid.

In the reverse approach, though, the use of angular 2D schemes is possible, but radial space does not need to be discretised. The points of the survey are grouped according to angular zones instead of being approximated on a 3D grid. An estimate of the spherical Bessel transform is computed in every solid angle, and the result is then expanded in spherical harmonics on the 2D spherical grid. In the final account, this method naturally leads to more accurate coefficients than the forward method and also takes advantage of high-performance 2D schemes. For these reasons, we focus on the reverse approach and its implementation, using HEALPix for the angular transform.

Finally, for both forward and reverse methods the spherical harmonics discretised basis (coefficients Yℓm(γp)) may be fully pre-computed and stored in external files. This is a particularly useful feature (incompatible with the original formulation of spherical Fourier-Bessel), which significantly eases and speeds up the use of these methods.

3.2. Speed and accuracy of the reverse method

To test the accuracy and speed of the reverse method compared to the raw method, we considered the high-resolution full-sky horizon simulation (Teyssier et al. 2009). Horizon is a N-body simulation covering a 2 h-1   Gpc periodic box using 70 billion dark matter particles using a WMAP3 cosmology (Spergel et al. 2007). A derived halo catalogue is available, which we used to calculate fℓmn and C(kℓn) values using both methods (raw and reverse). Since we are interested only in comparing the speed of each method, we simply considered each halo to have equal weight.

We performed the raw and the reverse estimates on three “surveys” of N = 4.2 × 105,3.1 × 106 and 1 × 107 halos, which correspond to three different depths (zmax = 0.1,0.2 and 0.3, respectively) in the horizon simulation. The HEALPix angular parameter is given by nside.

The results of the accuracy and speed tests are given in Table 1. The third (fourth) column gives the percentage f coefficients for which the relative accuracy ϵ(fℓmn) (ϵ(Cℓn)) is lower than 0.3% for given values of nside and N. We considered the intervals (l,n) ∈ ([0,20], [1,20]) and (l,n) ∈ ([21,50], [21,50]) separately, since the estimation of higher coefficients depends more on the value of nside. We also compared computation times of the two methods by observing the ratio T = treverse/traw. Given a survey and a method, computation time denotes the CPU time required to compute the kln’s (from the Bessel functions) and the final coefficients fℓmn without using pre-computed quantities. Note that we performed this analysis by distributing the calculations on five machines and simply adding the individual contributions to computation time since our method is linear with survey size. With the reverse method, though, the roots of the Bessel functions as well as the spherical harmonics may be pre-computed and stored in external files, which decreases computation time and complexity when working with 3DEX.

Table 1

Estimation of Fourier-Bessel coefficients: comparison of the new method, the reverse formulation (Eqs. (19) and (20) using HEALPix discretisation) with the original, raw formulation (Eq. (7)).

The reverse method is about an order of magnitude faster than the raw method, but this depends on the choice of nside. For nside = 1024 almost all fℓmn coefficients in the range  [0,20  (for and n) have relative error below 1%, and 90% have it below 0.3%, whereas over 99% of C(ℓ,kn) coefficients are accurate to  <0.3%. In the range  [20,50], the accuracy is somewhat degraded due to the extension of the HEALPix formalism to 3D surveys. Indeed, for data distributed on the sphere, 3D space is very sparse even for large surveys. Increasing nside to 1024 or 2048 strongly improves the accuracy for higher orders . Note that comparisons for  > 50 are limited by the amount of time the raw method takes.

Figures 1 and 2 show the time taken for calculations as a function of max and nmax (Fig. 1), and as a function of number of halos (Fig. 2). The boxes correspond to the raw method, the circles and diamonds to the reverse method with nside = 512,1024, respectively. The dashed line corresponds to the general rule that the raw method scales as Nnmaxmax(+1)2/2\hbox{$N n_{\rm max} \left(\ell_{\rm max}+1\right)^2/2$}, whereas the points are all estimated from calculations. With max = nmax = 100 and N = 9.7e6, the raw decomposition took a few days of calculations, whereas the reverse method only took 12 hours. In our formalism, kmax is determined by the choice of R for the boundary condition and by the band-limit nmax for spherical Bessel coefficients. For each multipole we have kmax = knmax = qnmax/R where qnmax is the nmax-the root of the th spherical Bessel function. Because R is usually imposed by the problem or the data, one must increase nmax to probe smaller radial scales. In fact, a reasonable approximation (or even a simple plot) shows that qℓn ≈ ( + 3n). This observation enabled one to predict which radial scales are probed and how computation time scales with kmax, given that we provide the complexity for max and nmax.

thumbnail Fig. 1

Speed results (raw and reverse methods) for increasing summation limits max = nmax, for a survey of N = 9.7 × 106 halos. Dashed lines are the fitted complexity curves. The reverse formulation is suitable for pre-calculations and external storage of the spherical harmonics, which was not performed here but enables for additional speed improvements.

thumbnail Fig. 2

Speed results (raw and reverse methods) for increasing survey size, for max = nmax = 30. Dashed lines are the fitted complexity curves.

One of the main advantages of the reverse method is that it is naturally suited to parallel computing because it uses HEALPix fast spherical harmonics transform routines. All previous tests were performed on a recent computer (i7 processor, 8 Go RAM) and take advantage of OpenMP (with four threads). More advanced computing means (larger RAM and more processors) significantly decrease calculation time. For example, max = nmax = 128 with nside = 2048 took about an hour with 128 cores and 512 Go RAM, whereas computation time for the raw method was estimated to several days on the same machine. Note that the raw method is also suited to parallelisation: galaxies may be treated separately by different threads. In all experiments, we took advantage of this property and performed both raw and reverse decompositions with four threads to perform relevant comparisons between the two.

In terms of the power spectrum, Figs. 3 and 4 show the relative error between the raw and the reverse methods both in mode-mode space ( − n) and in mode-scale space ( − kℓn). For this comparison we decomposed a survey of N = 4.2 × 105 halos with zmax = 0.1. Figures in mode-mode space naturally differ according to the choice of the boundary R because the latter determines the discrete radial scale spectrum  {kℓn}, and hence mode n computed with two different R’s corresponds to different k-scales. When comparing the results from R = 1000 and R = 2000 in mode-scale space, we observe that the boundary condition fixes the explored scales. The left column is thus complementary to the right column to explore higher values of k. Although Fig. 3 gives information about the final coefficients, Fig. 4 is hence more appropriate to see which scales are probed and with what accuracy.

In view of the  − kℓn space, we see that no fluctuations are observed along the k axis up to k = 0.03   h   Mpc-1. In this range, fluctuations occur in space, which are accurately probed with nside = 512 until  = 25 but naturally require a more precise scheme for  > 25, k > 0.03   h   Mpc-1 (smaller scales in physical space). In conclusion, parameter nside (as well as R) must be chosen depending on the scales one wishes to probe. Figures 3 and 4 provide accuracy results that are complementary to Table 1.

thumbnail Fig. 3

Relative error on the power spectrum in mode-mode space C(l,n). We compare the original formulation of spherical Fourier-Bessel decomposition with the reverse formulation, testing nside = 512,1024 (rows) and R = 1000,2000 (columns). Only a few zones (white spots) are outside the range  [−0.3%, + 0.3%].

thumbnail Fig. 4

Relative error on the power spectrum in mode-scale space C(l,k) (k is in h Mpc-1). We compare the original formulation of spherical Fourier-Bessel decomposition with the reverse formulation, testing nside = 512,1024 (rows) and R = 1000,2000 (columns). Only a few zones (white spots) are outside the range  [−0.3%, + 0.3%].

4. The 3DEX library

The 3DEX library requires the HEALPix package (v2.12 or later) and the CFITSIO library. 3DEX may either by installed with an HEALPix-like procedure (configure and make commands) or using CMake. The Fortran modules, the 3DEX dynamic library and the related executables will be created in the relevant directories (see README file for more information).

In addition to the numerical procedures required to compute Fourier-Bessel coefficients, various other routines are provided in the library, such as those converting redshift to comoving distance, computing spherical Bessel functions and their zeros, reading and writing 3D structures (flmn and Cln), or reconstructing radial maps from Fourier-Bessel coefficients.

Three executable programmes are generated:

  • survey2almn performs the spherical Fourier-Besseldecomposition (reverse method) of a discrete survey with inputparameters lmax, nmax, r and nside. Outputs are the flmn coefficients and the power spectrum;

  • survey2almn_interactive is very similar to the previous programme, but converts redshift values into comoving distance before performing the spherical Fourier-Bessel decomposition. In particular, the routine takes a .txt file as input, taking into account parameters on the cosmology and on the decomposition;

  • almn2rmap extracts the flmn coefficients from a FITS file and reconstructs the field (HEALPix map) at a given radius. Inputs are the resolution, the radius and summation limits lmax and nmax, which allows one to reconstruct several maps at different scales and resolutions.

The corresponding calls are given by the examples below.

 > survey2almn survey_thetaphir.dat almn.fits cln.fits 20 20 256 2000.0,

where survey_thetaphir.dat is a survey with columns representing θ,φ,r, and the keywords correspond to values of ℓ,n,nside,R. The output is both the coefficient values (almn.fis) and the Fourier-Bessel spectrum (cln.fits).

 > survey2almn_interactive parameters.txt,

where parameters.txt is an external file containing input parameters for the survey and the cosmology (which allows for more flexible use). Finally, for the map reconstruction, we can use:

 > almn2rmap almn.fits map.fits 400.0 256 10 10 2000.0,

where the keywords correspond to rmax, nside, , n, R.

5. Conclusion

High-precision cosmology from galaxy and weak lensing surveys will require the analysis of 3D data in spherical coordinates, a situation for which spherical Fourier-Bessel decomposition is most suited. Current methods will be inadequate for future planned cosmological surveys, which will provide for example galaxy surveys with billions of galaxies, compared to millions today.

We have reviewed the forward or SHB formalism of the spherical Fourier-Bessel decomposition which first calculates the tangential, then the radial decomposition. We also introduced the reverse or SBH formalism that inverses this order. Only the latter approach can take advantage of existing fast codes for the calculation of tangential modes. (To do the same, the former would require a a new voxelisation scheme.)

We presented a public code 3DEX (3D EXpansions) for the fast calculation of Fourier-Bessel coefficients and spectra, which uses the HEALPix pixelisation scheme for calculating the tangential modes. The 3DEX code is based on the reverse/SBH formulation of the Fourier-Bessel decomposition.

We tested the 3DEX code on linear and quasi-linear scales ( < 50 and kℓn < 0.2   h   Mpc-1) using the horizon halo simulation for redshifts z < 0.3. For nside = 1024 the 3DEX method for calculating the power spectrum C(ℓ,k) is accurate to 0.3% on these scales.

For surveys with  < 10 million galaxies, computation time is reduced by a factor 4–12 depending on the desired scales and accuracy. For larger surveys the gain in time will be even greater. Finally, the use of the 3DEX code is not restricted to cosmological calculations, and can be used in any other discipline that requires a spherical Fourier-Bessel analysis of 3D data.

Acknowledgments

The authors are grateful to Ofer Lahav, Pirin Erdoğdu and François Lanusse for useful discussions regarding the spherical Fourier-Bessel theory, as well as to Romain Teyssier and Nicolas Clerc for access and help using the horizon simulation. The authors thank Yves Revaz for computational help at EPFL. The 3DEX library uses Healpix software (Górski et al. 2002, 2005). This work is supported by the European Research Council grant SparseAstro (ERC-228261). This research is in part supported by the Swiss National Science Foundation (SNSF).

References

  1. Abramo, L. R., Reimberg, P. H., & Xavier, H. S. 2010, PRD, 82, 043510 [NASA ADS] [CrossRef] [Google Scholar]
  2. Annis, J., Bridle, S., Castander, F. J., et al. 2005, unpublished [arXiv:astro-ph/0510195] [Google Scholar]
  3. Baddour, N. 2010, J. Opt. Soc. Am., 27, 10 [Google Scholar]
  4. Binney, J., & Quinn, T. 1991, MNRAS, 249, 678 [NASA ADS] [Google Scholar]
  5. Castro, P. G., Heavens, A. F., & Kitching, T. D. 2005, Phys. Rev. D, 72, 023516 [NASA ADS] [CrossRef] [Google Scholar]
  6. Crittenden, R. G. 2000, Astrophys. Lett. Comm., 37, 377 [NASA ADS] [Google Scholar]
  7. Crittenden, R. G., & Turok, N. G. 1998 [arXiv:astro-ph/9806374] [Google Scholar]
  8. Doroshkevich, A., Naselsky, P., Verkhodanov, O., et al. 2005, IJMPD, 14, 275 [NASA ADS] [CrossRef] [Google Scholar]
  9. Erdoğdu (a), P., Huchra, J. P., Lahav, O., et al. 2006a, MNRAS, 368, 1515 [Google Scholar]
  10. Erdoğdu (b), P., Lahav, O., Huchra, J. P., et al. 2006b, MNRAS, 373, 45 [NASA ADS] [CrossRef] [Google Scholar]
  11. Fisher, K. B., Lahav, O., Hoffman, Y., Lynden-Bell, D., & Zaroubi, S. 1995, MNRAS, 272, 885 [NASA ADS] [Google Scholar]
  12. Górski, K. M., Banday, A. J., Hivon, E., & Wandelt, B. D. 2002, in Astronomical Data Analysis Software and Systems XI, ed. D. A. Bohlender, D. Durand, & T. H. Handley, ASP Conf. Ser., 281, 107 [Google Scholar]
  13. Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]
  14. Heavens, A. 2003, MNRAS, 343, 1327 [NASA ADS] [CrossRef] [Google Scholar]
  15. Heavens, A. F., & Taylor, A. N. 1995, MNRAS, 275, 483 [NASA ADS] [Google Scholar]
  16. Huchra, J. P., Macri, L. M., Masters, K. L., et al. 2012, ApJS, 199, 26 [NASA ADS] [CrossRef] [Google Scholar]
  17. Lanusse, F., Rassat, A., & Starck, J.-L. 2012, A&A, in press, DOI: 10.1051/0004-6361/201118568 [Google Scholar]
  18. Larson, D., Dunkley, J., Hinshaw, G., et al. 2011, ApJS, 192, 16 [NASA ADS] [CrossRef] [Google Scholar]
  19. Laureijs, R., Amiaux, J., Arduini, S., et al. 2011 [arXiv:astro-ph/1110.3193] [Google Scholar]
  20. Percival, W. J., Cole, S., Eisenstein, D. J., et al. 2007b, MNRAS, 381, 1053 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  21. Planck Collaboration 2006 [arXiv:astro-ph/0604069] [Google Scholar]
  22. Rassat, A., & Refregier, A. 2012, A&A, in press, DOI: 10.1051/0004.6361/201118638 [Google Scholar]
  23. Refregier, A., Amara, A., Kitching, T. D., et al. 2011, A&A, 528, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Schlegel, D. J., Blanton, M., Eisenstein, D., & et al. 2007, in AAS Meeting Abstracts, 211, 132.29 [Google Scholar]
  25. Schrabback, T., Hartlap, J., Joachimi, B., et al. 2010, A&A, 516, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Spergel, D. N., Bean, R., Doré, O., Nolta, M. R., et al. 2007, ApJS, 170, 377 [NASA ADS] [CrossRef] [Google Scholar]
  27. Strauss, M. A., Huchra, J. P., Davis, M., et al. 1992, ApJS, 83, 29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Tadros, H., Ballinger, W. E., Taylor, A. N., et al. 1999, MNRAS, 305, 527 [NASA ADS] [CrossRef] [Google Scholar]
  29. Teyssier, R., Pires, S., Prunet, S., et al. 2009, A&A, 497, 335 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Tyson, J. A., & LSST. 2004, in BAAS, 36, 108.01 [Google Scholar]

Appendix A: Normalisation and discrete radial spectrum

The basis functions kjl(kr)Ylm(θ,φ) form a set of eigenfunctions of the Laplacian operator in spherical coordinates. In particular, these functions are orthonormalised in the continuous case, thanks to the orthogonality relation dΩdrr2jl(kr)jl(kr)Ylm(θ,φ)Ylm(θ,φ)=π2kkδD(kk)δllKδmmK,\appendix \setcounter{section}{1} \begin{eqnarray} \nonumber \int {\rm d}\Omega {\rm d}r \ r^2 j_l(kr) j_{l'}(k' r) Y_{lm}(\theta,\phi) Y^{*}_{l'm'}(\theta,\phi) =\\\qquad \frac{\pi}{2kk'}\delta^D(k-k')\delta^K_{ll'}\delta^K_{mm'} \label{closure1}, \end{eqnarray}(A.1)(Baddour 2010) where δK is Kronecker’s delta notation and δD Dirac’s function.

A common approach to simplify the problem is to assume some boundary conditions for the field f. Different conditions have been explored in the literature (Fisher et al. 1995; Heavens & Taylor 1995), including potential or gradient continuity. In this paper, we used a condition that derives from the classical formulation of the discrete spherical Bessel transform: space is assumed to be finite and limited to a sphere of radius R. In this case, the spherical Bessel functions are not normalised and the boundary effect leads to a discrete spectrum  {kln}. The Fourier-Bessel coefficients become a set flmn = flm(kln), and the complete description of the field in the so-called Fourier-Bessel basis (Binney & Quinn 1991) is summarised in Eq. (5).

As a consequence, a natural choice for the boundary condition is to impose the field to vanish at r = R (Abramo et al. 2010), which constrains the Bessel functions and generates the radial spectrum  {kln}  such that, for all l and n, jl(klnR)=0.\appendix \setcounter{section}{1} \begin{eqnarray} j_{l}(k_{ln}R) = 0. \end{eqnarray}(A.2)If qln denotes the nth root of jl(z), the closure relation of the Bessel basis is 01dzz2jl(qlnz)jl(qlnz)=12[jl+1(qln)]2δllδnn,\appendix \setcounter{section}{1} \begin{equation} \int_0^1{\rm d}z \ z^2j_l(q_{ln}z)j_{l'}(q_{l'n'}z) = \frac{1}{2}[j_{l+1}(q_{ln})]^2 \delta_{ll'}\delta_{nn'}, \end{equation}(A.3)which gives with kln = qln/R, 0Rdrr2klnklnjl(klnr)jl(klnr)=kln2[jl+1(qln)]22R-3δllδnn.\appendix \setcounter{section}{1} \begin{equation} \int_0^R{\rm d}r \ r^2 k_{ln}k_{l'n'}j_l(k_{ln}r)j_{l'}(k_{l'n'}r) = \frac{k^2_{ln}[j_{l+1}(q_{ln})]^2}{2R^{-3}} \delta_{ll'}\delta_{nn'}. \end{equation}(A.4)The discrete spectrum is thus fixed by the zeros of the spherical Bessel functions. We obtain the normalisation coefficients κln (Fisher et al. 1995) κln-1=R32[klnjl+1(klnR)]2,\appendix \setcounter{section}{1} \begin{eqnarray} \kappa_{ln}^{-1} = \frac{R^{3}}{2}[k_{ln}j_{l+1}(k_{ln}R)]^2, \end{eqnarray}(A.5)which are used for field reconstruction (Eq. (5)).

Other approaches are possible to tackle boundary conditions in radial space, notably those imposing potential continuity at r = R (Fisher et al. 1995). Then, the discrete spectrum kln\hbox{$k'_{ln}$} is such that jl1(klnR)=0,\appendix \setcounter{section}{1} \begin{eqnarray} j_{l-1}(k'_{ln}R) = 0, \end{eqnarray}(A.6)and normalisation constraint becomes κln-1=R32[klnjl(klnR)]2.\appendix \setcounter{section}{1} \begin{eqnarray} {\kappa'}^{-1}_{ln} = \frac{R^{3}}{2}[k_{ln}j_{l}(k_{ln}R)]^2. \end{eqnarray}(A.7)

Appendix B: Angular masks

3DEX takes into account optional angular masks under the form of either an equatorial cut or an input all-sky FITS map.

In the first case, supplying θcut defines the latitude (in degrees) of a straight symmetric cut around the equator. Pixels located within that cut (l = cos(θcut)) are ignored.

In the second case, the supplied mask must be an HEALPix map (ring ordering) of Npix pixels at resolution nside (which must be identical to Fourier-Bessel resolution parameter) {w(γq)}q=1,...,Npix.\appendix \setcounter{section}{2} \begin{eqnarray} \left\{w(\boldsymbol{\gamma}_q)\right\}_{q=1,\dots,N_{\rm pix}}\!. \end{eqnarray}(B.1)In the forward method, the first step is to apply the mask to each discrete shell before performing the spherical harmonics transform. Hence for the ith shell, field f is weighted by the mask at each pixel f(ri,γq)=w(γq)f(ri,γq).\appendix \setcounter{section}{2} \begin{eqnarray} f'(r_i,\boldsymbol{\gamma}_q) = w(\boldsymbol{\gamma}_q) f(r_i,\boldsymbol{\gamma}_q). \end{eqnarray}(B.2)The Fourier-Bessel coefficients are obtained after performing SH and SB transforms.

In the reverse method, the first step is still the spherical Bessel transform, which gives a set of nmax HEALPix maps f(kln,γq). The mask is then applied to each of these maps f(kln,γq)=w(γq)f(kln,γq).\appendix \setcounter{section}{2} \begin{eqnarray} f'(k_{ln},\boldsymbol{\gamma}_q) = w(\boldsymbol{\gamma}_q) f(k_{ln},\boldsymbol{\gamma}_q). \end{eqnarray}(B.3)and the modified spherical harmonics transform gives the final flmn coefficients.

All Tables

Table 1

Estimation of Fourier-Bessel coefficients: comparison of the new method, the reverse formulation (Eqs. (19) and (20) using HEALPix discretisation) with the original, raw formulation (Eq. (7)).

All Figures

thumbnail Fig. 1

Speed results (raw and reverse methods) for increasing summation limits max = nmax, for a survey of N = 9.7 × 106 halos. Dashed lines are the fitted complexity curves. The reverse formulation is suitable for pre-calculations and external storage of the spherical harmonics, which was not performed here but enables for additional speed improvements.

In the text
thumbnail Fig. 2

Speed results (raw and reverse methods) for increasing survey size, for max = nmax = 30. Dashed lines are the fitted complexity curves.

In the text
thumbnail Fig. 3

Relative error on the power spectrum in mode-mode space C(l,n). We compare the original formulation of spherical Fourier-Bessel decomposition with the reverse formulation, testing nside = 512,1024 (rows) and R = 1000,2000 (columns). Only a few zones (white spots) are outside the range  [−0.3%, + 0.3%].

In the text
thumbnail Fig. 4

Relative error on the power spectrum in mode-scale space C(l,k) (k is in h Mpc-1). We compare the original formulation of spherical Fourier-Bessel decomposition with the reverse formulation, testing nside = 512,1024 (rows) and R = 1000,2000 (columns). Only a few zones (white spots) are outside the range  [−0.3%, + 0.3%].

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.