© 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  − k². 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 (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(θ,φ) =  ∑ _ℓ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) (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 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, Eq. (7) will become highly time-consuming 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 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 (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 Fourier-Bessel coefficients f_ℓm(k) (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 (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 Fourier-Bessel 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 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: (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 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 r₁,...,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 high-performance 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 Fourier-Bessel coefficients (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:

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 Fourier-Bessel 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, 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 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 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 . 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)²). 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 × 10⁵,3.1 × 10⁶ and 1 × 10⁷ 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 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.

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 band-limit 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 = 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.

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 = 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 mode-mode space (ℓ − n) and in mode-scale space (ℓ − k_ℓn). For this comparison we decomposed a survey of N = 4.2 × 10⁵ halos with z_max = 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 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 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%].

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 (f_lmn and C_ln), 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 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 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 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 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 r_max, n_side, ℓ, 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 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 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

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 Fourier-Bessel coefficients become a set f_lmn = f_lm(k_ln), 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  {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 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 N_pix pixels at resolution n_side (which must be identical to Fourier-Bessel 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 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 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 = 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

	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

	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

	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