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/00046361/201118463  
Published online  27 March 2012 
3DEX: a code for fast spherical FourierBessel decomposition of 3D surveys^{⋆}
^{1} Department of Physics and Astronomy, University College London, London WC1E 6BT, UK
email: boris.leistedt@gmail.com
^{2} Laboratoire d’Astrophysique, École Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland
email: anais.rassat@epfl.ch
^{3} Institute for Astronomy, ETH Zürich, WolfgangPauliStrasse 27, 8093 Zürich, Swtitzerland
^{4} Laboratoire AIM, UMR CEACNRSParis 7, Irfu, SAp/SEDI, Service d’Astrophysique, CEA Saclay, 91191 GifsurYvette Cedex, France
Received: 15 November 2011
Accepted: 6 February 2012
Context. Highprecision cosmology requires the analysis of largescale surveys in 3D spherical coordinates, i.e. spherical FourierBessel decomposition. Current methods are insufficient for future datasets from widefield cosmology surveys.
Aims. The aim of this paper is to present a public code for fast spherical FourierBessel decomposition that can be applied to cosmological data or 3D data in spherical coordinates in other scientific fields.
Methods. We present an equivalent formulation of the spherical FourierBessel decomposition that separates radial and tangential calculations. We propose to use the existing pixelisation scheme HEALPix for a rapid calculation of the tangential modes.
Results. 3DEX (3D EXpansions) is a public code for fast spherical FourierBessel decomposition of 3D allsky surveys that takes advantage of HEALPix for the calculation of tangential modes. We perform tests on very large simulations and we compare the precision and computation time of our method with an optimised implementation of the spherical FourierBessel original formulation. For surveys with millions of galaxies, computation time is reduced by a factor 4–12 depending on the desired scales and accuracy. The formulation is also suitable for precalculations and external storage of the spherical harmonics, which allows for additional speed improvements. The 3DEX code can accommodate data with masked regions of missing data. 3DEX can also be used in other disciplines, where 3D data are to be analysed in spherical coordinates.
Key words: cosmology: observations / largescale structure of Universe / methods: statistical / methods: data analysis / galaxies: statistics
The code and documentation can be downloaded at http://ixkael.com/blog/3dex
© ESO, 2012
1. Introduction
In the last few decades, cosmology has become a datadriven field, where highprecision 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 largescale structure with galaxy surveys. In order to be beneficial, cosmological studies of these surveys need to use highprecision statistical methods, such as a full 3D analysis on the sky where allsky 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 FourierBessel 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 timeconsuming.
In Sect. 2.1, we present the theory behind the 3D FourierBessel 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 FourierBessel 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 FourierBessel 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 FourierBessel normalisation.
2. Theory
2.1. The spherical FourierBessel 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 − k^{2}. For an homogeneous threedimensional field f(r) = f(r,θ,φ) in a flat geometry, the spherical FourierBessel decomposition (Fisher et al. 1995; Heavens 2003; Castro et al. 2005) is (1)with the inverse relation (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 (3)a naïve estimator of which is (4)This can be seen as an extension of the usual 2D power spectrum . The latter arises from the spherical harmonic transform of a 2D field given on the sphere f(θ,φ) = ∑ _{ℓm}f_{ℓm}Y_{ℓ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 FourierBessel decomposition becomes (Erdoğdu et al. 2006b; Fisher et al. 1995) (5)which is exact if the ranges of ℓ,m and n are infinite. The FourierBessel 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 nonlinearities, but these are not the purpose of this work. In this work we only consider linear or quasilinear scales (ℓ < 50, k < 0.2 h Mpc^{1}).
If the only information about the field f is a list of coordinates r_{p} = (r_{p},θ_{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) 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 Nn_{max}(ℓ_{max} + 1)^{2}/2, Eq. (7) will become highly timeconsuming when applied to larger surveys or when precise decomposition is required (large n_{max} 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 FourierBessel 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 (8)for which the inversion formula gives harmonics coefficients f_{ℓm}(r) depending on the radius r(9)It is then possible to perform a spherical Bessel transform (10)leading to the final FourierBessel coefficients f_{ℓm}(k) (11)This formulation hence extends the notion of 2D spherical harmonics to threedimensional 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 (12)for which the inversion formula gives (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: (14)Again, using the inversion formula, we obtain the FourierBessel coefficients (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 sphericalharmonicBessel, Eqs. (8) to (11)) and the reverse (denoted by SBH for sphericalBesselharmonic, Eqs. (12) to (15)) methods are equivalent formulations of the spherical FourierBessel decomposition of any threedimensional field f(r,θ,φ). This is summarised in the following schematic description of each method: (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 k_{ln}, bounding each integral and adapting normalisation. The formulas arising from this adaptation are used in the next sections.
2.3. Estimating FourierBessel 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 r_{1},...,r_{Nlayers}. In each layer, the coefficients f_{ℓm}(r_{i}) 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 N_{pix}(i) pixels for the ith shell) and to take advantage of the related highperformance algorithms. Angular space is hence discretised into nodes (r_{i},θ_{p},φ_{p}) = (r_{i},γ_{q}) and the field is approximated on each node, giving . The spherical harmonic decomposition in the ith shell becomes (17)and the final coefficients are obtained by performing the following spherical Bessel decomposition: (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 N_{pix} zones (pixels) related to angular nodes γ_{q}. If G_{q} 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 (19)and each of these intermediate maps is decomposed into spherical harmonic (spherical harmonics transform), which gives the FourierBessel coefficients (20)With the reverse method, one can avoid to discretise radial space. Moreover, this oneshell 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 FourierBessel approach, whereas the SHB formulation naturally handles the dependence on k of the weighting function.
The three methods to estimate the spherical FourierBessel decomposition can therefore also be expressed for a discrete 3D survey, summarised schematically below:
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,n_{max} , 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 FourierBessel coefficients. However, since each point contributes to the calculation of every coefficient (∀ l,m,n), computation time is proportional to N·n_{max}(ℓ_{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 highperformance 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 N_{pix}, 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 Ω(r_{i}). 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 FourierBessel coefficients. Furthermore, designing a meaningful radial discretisation is a difficult task. For equalarea pixelisations, the area of each pixel on the ith shell is . With HEALPix, the n_{side} angular parameter may only be increased by a factor 2, which changes the number of pixels by a factor of 4 (since N_{pix}(i) = 12n_{side}(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 equalvolume voxels. Extending 2D schemes to 3D is difficult and may even require a novel formalism for an equalvoxel 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 highperformance 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 precomputed and stored in external files. This is a particularly useful feature (incompatible with the original formulation of spherical FourierBessel), 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 highresolution fullsky horizon simulation (Teyssier et al. 2009). Horizon is a Nbody 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 × 10^{5},3.1 × 10^{6} and 1 × 10^{7} halos, which correspond to three different depths (z_{max} = 0.1,0.2 and 0.3, respectively) in the horizon simulation. The HEALPix angular parameter is given by n_{side}.
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 n_{side} 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 n_{side}. We also compared computation times of the two methods by observing the ratio T = t_{reverse}/t_{raw}. Given a survey and a method, computation time denotes the CPU time required to compute the k_{ln}’s (from the Bessel functions) and the final coefficients f_{ℓmn} without using precomputed 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 precomputed and stored in external files, which decreases computation time and complexity when working with 3DEX.
The reverse method is about an order of magnitude faster than the raw method, but this depends on the choice of n_{side}. For n_{side} = 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(ℓ,k_{n}) 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 n_{side} 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 n_{max} (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 n_{side} = 512,1024, respectively. The dashed line corresponds to the general rule that the raw method scales as , whereas the points are all estimated from calculations. With ℓ_{max} = n_{max} = 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, k_{max} is determined by the choice of R for the boundary condition and by the bandlimit n_{max} for spherical Bessel coefficients. For each multipole ℓ we have k_{max} = k_{ℓnmax} = q_{ℓnmax}/R where q_{ℓnmax} is the n_{max}the root of the ℓth spherical Bessel function. Because R is usually imposed by the problem or the data, one must increase n_{max} 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 k_{max}, given that we provide the complexity for ℓ_{max} and n_{max}.
Fig. 1 Speed results (raw and reverse methods) for increasing summation limits ℓ_{max} = n_{max}, for a survey of N = 9.7 × 10^{6} halos. Dashed lines are the fitted complexity curves. The reverse formulation is suitable for precalculations and external storage of the spherical harmonics, which was not performed here but enables for additional speed improvements. 
Fig. 2 Speed results (raw and reverse methods) for increasing survey size, for ℓ_{max} = n_{max} = 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} = n_{max} = 128 with n_{side} = 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 modemode space (ℓ − n) and in modescale space (ℓ − k_{ℓn}). For this comparison we decomposed a survey of N = 4.2 × 10^{5} halos with z_{max} = 0.1. Figures in modemode 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 kscales. When comparing the results from R = 1000 and R = 2000 in modescale 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 n_{side} = 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 n_{side} (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.
Fig. 3 Relative error on the power spectrum in modemode space C(l,n). We compare the original formulation of spherical FourierBessel decomposition with the reverse formulation, testing n_{side} = 512,1024 (rows) and R = 1000,2000 (columns). Only a few zones (white spots) are outside the range [−0.3%, + 0.3%]. 
Fig. 4 Relative error on the power spectrum in modescale space C(l,k) (k is in h Mpc^{1}). We compare the original formulation of spherical FourierBessel decomposition with the reverse formulation, testing n_{side} = 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 HEALPixlike 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 FourierBessel 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 (f_{lmn} and C_{ln}), or reconstructing radial maps from FourierBessel coefficients.
Three executable programmes are generated:

survey2almn performs the spherical FourierBesseldecomposition (reverse method) of a discrete survey with inputparameters l_{max}, n_{max}, r and n_{side}. Outputs are the f_{lmn} 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 FourierBessel 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 f_{lmn} coefficients from a FITS file and reconstructs the field (HEALPix map) at a given radius. Inputs are the resolution, the radius and summation limits l_{max} and n_{max}, 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,n_{side},R. The output is both the coefficient values (almn.fis) and the FourierBessel 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 r_{max}, n_{side}, ℓ, n, R.
5. Conclusion
Highprecision cosmology from galaxy and weak lensing surveys will require the analysis of 3D data in spherical coordinates, a situation for which spherical FourierBessel 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 FourierBessel 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 FourierBessel 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 FourierBessel decomposition.
We tested the 3DEX code on linear and quasilinear scales (ℓ < 50 and k_{ℓn} < 0.2 h Mpc^{1}) using the horizon halo simulation for redshifts z < 0.3. For n_{side} = 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 FourierBessel 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 FourierBessel 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 (ERC228261). This research is in part supported by the Swiss National Science Foundation (SNSF).
References
 Abramo, L. R., Reimberg, P. H., & Xavier, H. S. 2010, PRD, 82, 043510 [NASA ADS] [CrossRef] [Google Scholar]
 Annis, J., Bridle, S., Castander, F. J., et al. 2005, unpublished [arXiv:astroph/0510195] [Google Scholar]
 Baddour, N. 2010, J. Opt. Soc. Am., 27, 10 [Google Scholar]
 Binney, J., & Quinn, T. 1991, MNRAS, 249, 678 [NASA ADS] [Google Scholar]
 Castro, P. G., Heavens, A. F., & Kitching, T. D. 2005, Phys. Rev. D, 72, 023516 [NASA ADS] [CrossRef] [Google Scholar]
 Crittenden, R. G. 2000, Astrophys. Lett. Comm., 37, 377 [NASA ADS] [Google Scholar]
 Crittenden, R. G., & Turok, N. G. 1998 [arXiv:astroph/9806374] [Google Scholar]
 Doroshkevich, A., Naselsky, P., Verkhodanov, O., et al. 2005, IJMPD, 14, 275 [NASA ADS] [CrossRef] [Google Scholar]
 Erdoğdu (a), P., Huchra, J. P., Lahav, O., et al. 2006a, MNRAS, 368, 1515 [Google Scholar]
 Erdoğdu (b), P., Lahav, O., Huchra, J. P., et al. 2006b, MNRAS, 373, 45 [NASA ADS] [CrossRef] [Google Scholar]
 Fisher, K. B., Lahav, O., Hoffman, Y., LyndenBell, D., & Zaroubi, S. 1995, MNRAS, 272, 885 [NASA ADS] [Google Scholar]
 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]
 Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Heavens, A. 2003, MNRAS, 343, 1327 [NASA ADS] [CrossRef] [Google Scholar]
 Heavens, A. F., & Taylor, A. N. 1995, MNRAS, 275, 483 [NASA ADS] [Google Scholar]
 Huchra, J. P., Macri, L. M., Masters, K. L., et al. 2012, ApJS, 199, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Lanusse, F., Rassat, A., & Starck, J.L. 2012, A&A, in press, DOI: 10.1051/00046361/201118568 [Google Scholar]
 Larson, D., Dunkley, J., Hinshaw, G., et al. 2011, ApJS, 192, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Laureijs, R., Amiaux, J., Arduini, S., et al. 2011 [arXiv:astroph/1110.3193] [Google Scholar]
 Percival, W. J., Cole, S., Eisenstein, D. J., et al. 2007b, MNRAS, 381, 1053 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Planck Collaboration 2006 [arXiv:astroph/0604069] [Google Scholar]
 Rassat, A., & Refregier, A. 2012, A&A, in press, DOI: 10.1051/0004.6361/201118638 [Google Scholar]
 Refregier, A., Amara, A., Kitching, T. D., et al. 2011, A&A, 528, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schlegel, D. J., Blanton, M., Eisenstein, D., & et al. 2007, in AAS Meeting Abstracts, 211, 132.29 [Google Scholar]
 Schrabback, T., Hartlap, J., Joachimi, B., et al. 2010, A&A, 516, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Spergel, D. N., Bean, R., Doré, O., Nolta, M. R., et al. 2007, ApJS, 170, 377 [NASA ADS] [CrossRef] [Google Scholar]
 Strauss, M. A., Huchra, J. P., Davis, M., et al. 1992, ApJS, 83, 29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tadros, H., Ballinger, W. E., Taylor, A. N., et al. 1999, MNRAS, 305, 527 [NASA ADS] [CrossRef] [Google Scholar]
 Teyssier, R., Pires, S., Prunet, S., et al. 2009, A&A, 497, 335 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tyson, J. A., & LSST. 2004, in BAAS, 36, 108.01 [Google Scholar]
Appendix A: Normalisation and discrete radial spectrum
The basis functions kj_{l}(kr)Y_{lm}(θ,φ) 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 (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 {k_{ln}}. The FourierBessel coefficients become a set f_{lmn} = f_{lm}(k_{ln}), and the complete description of the field in the socalled FourierBessel 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 {k_{ln}} such that, for all l and n, (A.2)If q_{ln} denotes the nth root of j_{l}(z), the closure relation of the Bessel basis is (A.3)which gives with k_{ln} = q_{ln}/R, (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) (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 is such that (A.6)and normalisation constraint becomes (A.7)
Appendix B: Angular masks
3DEX takes into account optional angular masks under the form of either an equatorial cut or an input allsky 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 N_{pix} pixels at resolution n_{side} (which must be identical to FourierBessel resolution parameter) (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 (B.2)The FourierBessel 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 n_{max} HEALPix maps f(k_{ln},γ_{q}). The mask is then applied to each of these maps (B.3)and the modified spherical harmonics transform gives the final f_{lmn} coefficients.
All Tables
All Figures
Fig. 1 Speed results (raw and reverse methods) for increasing summation limits ℓ_{max} = n_{max}, for a survey of N = 9.7 × 10^{6} halos. Dashed lines are the fitted complexity curves. The reverse formulation is suitable for precalculations and external storage of the spherical harmonics, which was not performed here but enables for additional speed improvements. 

In the text 
Fig. 2 Speed results (raw and reverse methods) for increasing survey size, for ℓ_{max} = n_{max} = 30. Dashed lines are the fitted complexity curves. 

In the text 
Fig. 3 Relative error on the power spectrum in modemode space C(l,n). We compare the original formulation of spherical FourierBessel decomposition with the reverse formulation, testing n_{side} = 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 
Fig. 4 Relative error on the power spectrum in modescale space C(l,k) (k is in h Mpc^{1}). We compare the original formulation of spherical FourierBessel decomposition with the reverse formulation, testing n_{side} = 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 (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.