SuperFaB: a fabulous code for Spherical Fourier-Bessel decomposition
SSuperFaB: a fabulous code for Spherical Fourier-Bessel decomposition
Henry S. Grasshorn Gebhardt
1, 2, ˚ and Olivier Doré
1, 2 Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA California Institute of Technology, Pasadena, CA 91125, USA
The spherical Fourier-Bessel (SFB) decomposition is a natural choice for the radial/angular sepa-ration that allows optimal extraction of cosmological information from large volume galaxy surveys.In this paper we develop a SFB power spectrum estimator that allows the measurement of thelargest angular and radial modes with the next generation of galaxy surveys. The code measuresthe pseudo-SFB power spectrum, and takes into account mask, selection function, pixel window,and shot noise. We show that the local average effect is significant only in the largest-scale mode,and we provide an analytical covariance matrix. By imposing boundary conditions at the minimumand maximum radius encompassing the survey volume, the estimator does not suffer from the nu-merical instabilities that have proven challenging in the past. The estimator is demonstrated onsimplified
Roman -like,
SPHEREx -like, and
Euclid -like mask and selection functions. For intuitionand validation, we also explore the SFB power spectrum in the Limber approximation. We releasethe associated public code written in
Julia . Keywords: cosmology; large-scale structure
I. INTRODUCTION
One of the aims of future galaxy surveys such as the
Nancy Grace Roman Space Telescope , SPHEREx , Eu-clid , DESI, PFS (Spergel et al. et al. et al. et al. et al. et al. k ´ in the simplestmodels (Salopek and Bond 1990, Komatsu and Spergel2001, Dalal et al. et al. angular scales are difficult to exploit fully with astandard 3D power spectrum analysis due to line-of-sighteffects such as redshift-space distortions (RSD). Whenthe angular separation between galaxies is large, the as-sumption that a single line of sight can be used for bothgalaxies breaks down, which results in a loss of informa-tion from the measurement. On medium large scales theproblem can be mitigated by choosing a common line ofsight for each pair of galaxies (Yamamoto et al. et al. ˚ [email protected]; NASA Postdoctoral ProgramFellow Large radial scales pose a different kind of challenge,and in the past have mostly been treated by splittingsurveys into redshift bins (e.g., Beutler et al. et al. (2016).The spherical Fourier-Bessel transform for the analysisof galaxy surveys has been considered multiple times inthe past. Binney and Quinn (1991) used a SFB decompo-sition to characterize overdensities deep in the nonlinearregime. Lahav (1993) applied the SFB analysis to lo-cal galaxies on larger scales. Heavens and Taylor (1995)applied the SFB analysis to the IRAS 1.2-Jy galaxy cata-logue, Tadros et al. (1999) applied it to the PSCz galaxy a r X i v : . [ a s t r o - ph . C O ] F e b catalogue, and Percival et al. (2004) use it in the con-text of the 2dF Galaxy Redshift Survey. Leistedt et al. (2012) have provided a public SFB code, 3DEX, that per-forms the SFB decomposition first in the radial directionfor each galaxy individually, then performs the angulartransform using HEALPix. More recently, Wang et al. (2020) have built a combined SFB/ P p k q estimation codethat uses SFB on very large scales, and a Cartesian mul-tipole power spectrum estimator on smaller scales.SFB power spectrum measurements tended to beplagued by numerical instabilities and computationalcomplexity. The source of the numerical instability isthe incomplete coverage of the analysis volume by thesurvey. For example, typically a boundary condition isapplied at some distance r max , and the analysis is per-formed for the entire volume inside a sphere of radius r max . However, most surveys will leave most of that vol-ume unexplored, and the de-convolution of the windowfunction or the inversion of the covariance matrix be-come numerically unstable. The numerical complexitystems from the large number of modes that need to becalculated even when a large fraction of the analysis vol-ume remains empty. Another source of computationalcomplexity is the combined estimation of the real-spacepower spectrum and redshift-space distortion parametersthat requires repeated estimations of the power spectrum(Wang et al. SuperFaB , combines several ap-proaches to address these problems. For the first time,we limit the redshift range by introducing a boundarycondition at r min as suggested by Samushia (2019). Wealso use the 3DEX approach by Leistedt et al. (2012) thatdoes not suffer from pixel window effects in the radial di-rection. The 3DEX approach also allows separation ofthe angular and radial transforms, and for the angulartransform we use HEALPix (Górski et al. et al. et al. et al. et al. et al. et al. (2015) conclude that SFB yields similarconstraints as SHT, but when it comes to marginalizingover systematic biases such as evolving scale-dependentgalaxy bias, SFB performs better. Additionally, Cas-torina and White (2018) developed various approachesto incorporating wide-angle effects in Fourier based es-timators. Beutler et al. (2019) implement a small-angleexpansion for the standard multipole power spectrum.In section Section II we review the SFB power spec-trum and develop intuition in the Limber approxima-tion. Section III details the approach taken for theSFB decomposition, window de-convolution, shot noise,bandpower binning, local average effect, and covari- ance matrix. We show comparisons with log-normalsimulations in Section IV for Roman , SPHEREx , and
Euclid , and we conclude in Section V. We leave tothe appendices a collection of useful formulae in Ap-pendix A, review the Laplacian in an expanding uni-verse in Appendix B, derive the radial potential bound-ary conditions in Appendix C, and simplify the covari-ance matrix in Appendices D and E. Our
SuperFaB codeis publically available at https://github.com/hsgg/SphericalFourierBesselDecompositions.jl . II. SFB POWER SPECTRUM
In this section we develop the SFB formalism. Westart with the basic transformation between configura-tion space and SFB space as well as between Fourierspace and SFB space. We then briefly show the powerspectrum in a completely homogeneous and isotropic uni-verse before adding in selection function, linear growthfactor, galaxy bias, and RSD. We develop intuition byapplying Limber’s approximation.The spherical Fourier-Bessel decomposition expressesa field δ p r q in terms of eigenfunctions of the Laplacianin spherical coordinates. For more details, we refer thereader to Section III A. We define the spherical Fourier-Bessel modes δ (cid:96)m p k q by δ p r q “ ż d k ÿ (cid:96)m «c π k j (cid:96) p kr q Y (cid:96)m p θ, φ q ff δ (cid:96)m p k q , (1) δ (cid:96)m p k q “ ż d r «c π k j (cid:96) p kr q Y ˚ (cid:96)m p ˆ r q ff δ p r q , (2)where r “ r ˆ r is the position vector, r is the comovingangular diameter distance from the origin, and ˆ r is the di-rection on the sky. The orthogonality relations Eqs. (A6)and (A7) for the spherical Bessel functions and spheri-cal harmonics are used to prove that Eqs. (1) and (2)are inverses of each other. The factor k { π can be splitbetween Eqs. (1) and (2) as pleased. Here we use theconvention in Nicola et al. (2014), because for a non-evolving, homogeneous, and isotropic universe the SFBpower spectrum then equals P p k q , see Eq. (7) below.The relation between the SFB coefficients δ (cid:96)m p k q andthe Fourier modes δ p k q is obtained by expressing δ p r q interms of its Fourier transform in Eq. (2), δ (cid:96)m p k q “ c π k ż d r j (cid:96) p kr q Y ˚ (cid:96)m p ˆ r q ż d q p π q e i q ¨ r δ p q q . (3)With Rayleigh’s formula Eq. (A10) this turns into δ (cid:96)m p k q “ ż d q p π q c π q π ÿ (cid:96) m i (cid:96) Y ˚ (cid:96) m p ˆ q q δ p q qˆ kqπ ż d r r j (cid:96) p kr q j (cid:96) p qr qˆ ż d ˆ r Y ˚ (cid:96)m p ˆ r q Y (cid:96) m p ˆ r q“ k p π q i (cid:96) ż d ˆ k Y ˚ (cid:96)m p ˆ k q δ p k q , (4)where the orthogonality relations Eqs. (A6) and (A7)were used. Eq. (4) shows that SFB is a spherical har-monic transform of Cartesian Fourier modes with an ad-ditional phase factor i (cid:96) . And δ p k q “ p π q k ÿ (cid:96)m i ´ (cid:96) Y (cid:96)m p ˆ k q δ (cid:96)m p k q (5)is the inverse of Eq. (4). A. The Homogeneous and Isotropic Universe
In a homogeneous and isotropic universe in real-space(with no line-of-sight effects), we have @ δ p k q δ ˚ p k q D “ p π q δ D p k ´ k q P p k q . (6)Therefore, applying Eq. (4) gives the SFB power spec-trum as @ δ (cid:96)m p k q δ ˚ (cid:96) m p k q D “ δ D p k ´ k q δ K(cid:96)(cid:96) δ Kmm P p k q , (7)where we used Eq. (A2) for the three-dimensional Dirac-delta function in spherical coordinates. That is, in a ho-mogeneous and isotropic universe with no observationaleffects the SFB power spectrum equals the 3D powerspectrum P p k q . B. The Linear Universe
We now generalize to include line-of-sight effects, alinearly evolving power spectrum, and a radial windowfunction. The galaxy density contrast we consider is δ obs g p r q “ W p r q D p r q ż d q p π q e i q ¨ r r A RSD p µ, qµ, r qˆ b p r, q q δ p q q , (8)where δ p q q is the matter density contrast in Fourierspace, W p r q is the survey window function, D p r q isthe linear growth factor, b p r, q q is the possibly scale-dependent linear galaxy bias, µ “ ˆ q ¨ ˆ r , and the redshift-space distortions are encoded in r A RSD p µ, qµ, r q “ ` ` βµ ˘ r A FoG p qµ q , (9)with β “ f { b , where f “ d ln D { d ln a is the linear growthrate, and we assume a Gaussian fingers-of-God term r A FoG p qµ q “ e ´ σ u q µ , (10)with σ u “ σ v { aH the pair-wise velocity dispersion inunits of length. The tilde on A RSD signifies that it is aFourier-space function.The RSD term r A RSD in Eq. (8) can be expressed as afunction of derivatives on the complex exponential. Thatis, by performing a Taylor series expansion we can replace µ Ñ ´ i B qr , or r A RSD p µ, qµ, r q e i q ¨ r “ ÿ n a n p q, r q n ! µ n e iqrµ “ ÿ n a n p q, r q n ! ˆ ´ i BBp qr q ˙ n e iqrµ “ r A RSD p´ i B qr , ´ iq B qr , r q e i q ¨ r . (11)Furthermore, the complex exponential is expanded us-ing Rayleigh’s formula Eq. (A10) so that the derivativesin r A RSD only act on the spherical Bessel function fromRayleigh’s formula. Further expressing the Fourier-spacedensity contrast in terms of its SFB modes Eq. (5), theobserved density contrast Eq. (8) now becomes δ obs g p r q “ W p r q D p r q ż d q p π q b p r, q q « r A RSD p´ i B qr , ´ iq B qr , r q π ÿ L M i L j L p qr q Y ˚ L M p ˆ q q Y L M p ˆ r q ff ˆ p π q q ÿ LM i ´ L Y LM p ˆ q q δ LM p q q (12) “ W p r q D p r q c π ż d q q b p r, q q « r A RSD p´ i B qr , ´ iq B qr , r q ÿ LM j L p qr q Y LM p ˆ r q ff q δ LM p q q . (13)Using Eq. (2) to transform into SFB space, δ g, obs (cid:96)m p k q “ ż d q ÿ LM W LM(cid:96)m p k, q q δ LM p q q , (14)where W LM(cid:96)m p k, q q “ ż d ˆ r Y LM p ˆ r q Y ˚ (cid:96)m p ˆ r q W L(cid:96) p k, q, ˆ r q , (15)and W L(cid:96) p k, q, ˆ r q “ qkπ ż d r r W p r q D p r q b p r, q q j (cid:96) p kr qˆ r A RSD p´ i B qr , ´ iq B qr , r q j L p qr q . (16)The SFB correlation function is, therefore, A δ g, obs (cid:96)m p k q δ g, obs , ˚ (cid:96) m p k q E “ ż d q ÿ LM W LM(cid:96)m p k, q q W LM, ˚ (cid:96) m p k , q q P p q q , (17)where we used Eqs. (7) and (14).Here we will only consider a radial selection function,as the angular mask will be handled in the estimator.Then, W p r q “ φ p r q , (18)and we define the simplification of Eq. (16) W (cid:96) p k, q q “ W (cid:96)(cid:96) p k, q, ˆ r q , (19)which is then independent of the direction ˆ r . Eq. (15)and Eq. (17) then simplify to A δ g, obs (cid:96)m p k q δ g, obs , ˚ (cid:96) m p k q E “ δ K(cid:96)(cid:96) δ Kmm C (cid:96) p k, k q , (20)with the SFB power spectrum defined as C (cid:96) p k, k q “ ż d q W (cid:96) p k, q q W ˚ (cid:96) p k , q q P p q q . (21)Eqs. (16), (19) and (21) show that RSD and linear growthcan be taken into account by a change in the radial win-dow function.Note that for a homogeneous and isotropic universewithout selection function, W p r q “ D p r q “ r A RSD “ and b p r, q q “ const , the window becomes W (cid:96) p k, q q “ δ D p k ´ q q , and Eq. (7) is reproduced. Also, note that W (cid:96) p k, q q is real, because the imaginary parts only appearwith even powers in r A RSD . C. Limber’s Approximation
In this section we aim to gain intuition for the SFBpower spectrum in Eq. (21) by applying Limber’s approx-imation. To do so, we will also approximate the effect of the FoG. We write Eqs. (9) and (10) acting on a sphericalBessel function as r A RSD p´ iq B qr , ´ i B qr , r q j (cid:96) p qr q“ e σ u q B qr ` ´ β B qr ˘ j (cid:96) p qr q . (22)The second derivative is obtained exactly via a recursionrelation for the derivative of Spherical Bessel function, ` ´ β B qr ˘ j (cid:96) p qr q “ ` ´ βf (cid:96) ˘ j (cid:96) p qr q ´ βf (cid:96) ´ j (cid:96) ´ p qr q´ βf (cid:96) j (cid:96) ` p qr q (23) “ ÿ ∆ (cid:96) ` δ K ∆ (cid:96), ´ βf (cid:96) ∆ (cid:96) ˘ j (cid:96) ` ∆ (cid:96) p qr q , (24)where the only non-zero f (cid:96) ∆ (cid:96) are f (cid:96) ´ “ (cid:96) p (cid:96) ´ qp (cid:96) ´ qp (cid:96) ` q , (25) f (cid:96) “ ´ (cid:96) ` (cid:96) ´ p (cid:96) ´ qp (cid:96) ` q , (26) f (cid:96) “ p (cid:96) ` qp (cid:96) ` qp (cid:96) ` qp (cid:96) ` q . (27)The FoG term acting on the spherical Bessel function isa convolution r A FoG p´ iq B qr q j (cid:96) p qr q “ ż d k π e ikqr r A FoG p qk q r j (cid:96) p k q“ ż d y A FoG p r ´ y q j (cid:96) p qy q (28) “ ż d y ? π σ u e ´ p r ´ y q σ u j (cid:96) p qy q , (29)where the tilde signify Fourier transforms, and we tookthe inverse transform of Eq. (10) (see e.g. Grasshorn Geb-hardt and Jeong 2020). As a first approximation, if thefrequency q is low, then the convolution will have littleeffect. If the frequency is high, the convolution will erasethe oscillations to vanish. That is, we approximate r A FoG p´ iq B qr q j (cid:96) p qr q « e ´ σ u q j (cid:96) p qr q . (30)We are now in a position to apply Limber’s approxi-mation (LoVerde and Afshordi 2008) j (cid:96) p kr q Ñ c π kr δ D ˆ kr ´ (cid:96) ´ ˙ . (31) r in h − Mpc01234 D ( r ) f ( r ) σ u ( r ) / [ h − Mpc] b ( r ) φ ( r ) − − − k in h Mpc − C ‘ ( k ) i n [ h − M p c ] ‘ = 5 ‘ = 10 ‘ = 20 ‘ = 40 ‘ = 80 ‘ = 160 ‘ = 320 ≤ r Mpc /h ≤ P ( k )10 − − − − k in h Mpc − C ‘ ( k ) i n [ h − M p c ] r = ‘ +0 . k P ( k ) r = 1000 h − Mpc ( z = 0 . r = 2000 h − Mpc ( z = 0 . r = 3000 h − Mpc ( z = 1 . r = 4000 h − Mpc ( z = 2 . − − − k in h Mpc − ‘ r = M p c / h r = M p c / h C ‘ ( k ) . . × . × . × . × . × . × . × . × FIG. 1. Top left: Linear growth factor D p r q , linear growth rate f p r q , velocity dispersion σ u “ σ v {p aH q , galaxy bias b p r q “ b { D p r q where b “ , and selection function φ p r q defined in Eq. (18). Top right: The SFB power spectrum in theLimber approximation closely traces the 3D power spectrum. However, for a given perpendicular mode (cid:96) and redshift range,only a part of the power spectrum is measured. The horizontal lines show the range of k modes that a given (cid:96) mode is ableto measure for a survey within ď rh ´ Mpc ď , and the shaded bands show an estimate for the σ measurementuncertainty for that particular (cid:96) -mode. Bottom left: Here, each line fixes the redshift, and all (cid:96) modes are used. The Kaisereffect is not visible due to the Limber approximation becoming invalid on large scales. Bottom right: Here we show the SFBpower spectrum on a grid of (cid:96) - k modes. The SFB power spectrum can be measured within a band such that r » ` (cid:96) ` ˘ { k iswithin the survey. Then, Eqs. (16) and (19) become W (cid:96) p k, q q “ qkπ ż d r r φ p r q D p r q b p r, q q j (cid:96) p kr qˆ e ´ σ u q ÿ ∆ (cid:96) ` δ K ∆ (cid:96), ´ βf (cid:96) ∆ (cid:96) ˘ j (cid:96) ` ∆ (cid:96) p qr q (32) “ c qk φ ˆ (cid:96) ` k ˙ D ˆ (cid:96) ` k ˙ b ˆ (cid:96) ` k , q ˙ ˆ e ´ σ u q ÿ ∆ (cid:96) ` δ K ∆ (cid:96), ´ βf (cid:96) ∆ (cid:96) ˘ ˆ δ D ˆ q ´ (cid:96) ` ∆ (cid:96) ` (cid:96) ` k ˙ . (33) Since the Limber approximation is only applicable forlarge (cid:96) , we further assume ∆ (cid:96) ! (cid:96) . Then, W (cid:96) p k, q q “ δ D p q ´ k q φ ˆ (cid:96) ` k ˙ D ˆ (cid:96) ` k ˙ b ˆ (cid:96) ` k , k ˙ ˆ e ´ σ u k ÿ ∆ (cid:96) ` δ K ∆ (cid:96), ´ βf (cid:96) ∆ (cid:96) ˘ . (34)Therefore, the SFB power spectrum Eq. (21) in the Lim- r in h − Mpc01234 D ( r ) f ( r ) σ u ( r ) / [ h − Mpc] b ( r ) φ ( r ) − − − k in h Mpc − C ‘ ( k ) i n [ h − M p c ] ‘ = 5 ‘ = 10 ‘ = 20 ‘ = 40 ‘ = 80 ‘ = 160 ‘ = 320 ≤ r Mpc /h ≤ P ( k )10 − − − − k in h Mpc − C ‘ ( k ) i n [ h − M p c ] r = ‘ +0 . k P ( k ) r = 1000 h − Mpc ( z = 0 . r = 2000 h − Mpc ( z = 0 . r = 3000 h − Mpc ( z = 1 . r = 4000 h − Mpc ( z = 2 . − − − k in h Mpc − ‘ r = M p c / h r = M p c / h C ‘ ( k ) . . × . × . × . × . × . × . × . × FIG. 2. Same as Fig. 1, except that the linear galaxy bias is now constant b p r q “ . . As a result, the redshift evolution of thelinear growth factor D p r q is no longer cancelled by the bias. Top right: At a fixed (cid:96) , larger-scale modes are probed at higherredshift where the growth factor is smaller. Thus, compared to Fig. 1, each (cid:96) -segment appears tilted. Bottom left: The redshiftevolution of the linear growth factor D p r q causes a shift in the power spectrum amplitude with the Limber ratio r “ p (cid:96) ` q { k .Bottom right: High- (cid:96) modes are suppressed because they primarily probe the high redshifts. ber approximation is C (cid:96) p k, k q “ P p k q e ´ σ u k δ D ` k ´ k ˘ ˆ φ ˆ (cid:96) ` k ˙ D ˆ (cid:96) ` k ˙ b ˆ (cid:96) ` k , k ˙ ˆ “ ´ β ` f (cid:96) ´ ` f (cid:96) ` f (cid:96) ˘‰ . (35)The exponential is the suppression due to the FoG.The Dirac-delta function shows that even with redshift-evolution most of the power is on the diagonal k “ k , asfor a non-evolving universe. Redshift evolution manifestsitself mainly through the interplay between (cid:96) and k suchthat in the Limber approximation the ratio r “ (cid:96) ` k (36)is the comoving angular diameter distance. For example,if the scale k is fixed, then changing the angular scale (cid:96) corresponds to changing the redshift. We note thatprimordial non-Gaussianity will lead to a scale-dependentbias that will be absorbed directly in the bias term onthe second line. Finally, the last line in Eq. (35) accountsfor the linear Kaiser effect.In Fig. 1 we show the SFB power spectrum in theLimber approximation. We define C (cid:96) p k, k q “ δ D p k ´ k q C (cid:96) p k q . For the galaxy bias we choose b p r, k q “ b { D p r q with b “ so that the bias and linear growth factor can-cel, and our selection function is a constant φ p r q “ forillustration. The top left panel shows the inputs to ourcalculation.The top right panel of Fig. 1 shows the SFB powerspectrum for several fixed (cid:96) modes. The shaded areascorrespond to the σ measurement uncertainty estimated − − − − k in h Mpc − C ‘ ( k ) i n [ h − M p c ] ‘ = ‘ = ‘ = ‘ = ‘ = ‘ = ‘ = ‘ = ‘ = ‘ = ‘ = ‘ = ≤ r Mpc /h ≤ r = 1000 h − Mpc r = 4000 h − Mpc
FIG. 3. Here we combine the top right and bottom left plotsin Fig. 2 to show the region that the SFB power spectrumprobes in the Limber approximation for a given redshift rangeand linear power spectrum evolution. via ∆ C (cid:96) p k q “ c (cid:96) ` ˆ C (cid:96) p k q ` n ˙ , (37)(Pratten and Munshi 2013). Since (cid:96) corresponds to per-pendicular wave modes, at small r the corresponding k modes are large, and at large r the corresponding k modes are small. Therefore, at a constant (cid:96) , the SFBpower spectrum as a function of k sweeps through bothredshift and k modes, measuring a redshift correspond-ing to r » h ´ Mpc at lower k , and a redshift corre-sponding to r » h ´ Mpc at higher k .Consequently, in the (cid:96) - k plane, only a band of modescan be measured, as illustrated in the bottom right plot ofFig. 1. Fixing the redshift, which is possible in the Lim-ber approximation, results in a bona-fide power spectrummeasurement, as illustrated in the bottom left panel.In the Limber approximation, the SFB power spectrumdoes not exhibit a strong Kaiser effect on large scales. Weattribute this to the fact that Limber’s approximation isvalid only in the small-scale limit where the Kaiser effectis negligible.Because Eq. (36) relates the (cid:96) and k modes to a defi-nite redshift, we can only measure a band of modes. Weillustrate this in the bottom right panel of Fig. 1.The choice b p r, q q 9 D ´ p r q results in linear bias andlinear growth cancelling each other. In general, this maynot be the case, and we illustrate redshift evolution bysetting the bias constant, b p r, q q “ . , in Fig. 2. Becausefor fixed (cid:96) larger-scale modes k correspond to higher red-shift, the linear growth evolution tilts each (cid:96) -segment ofthe power spectrum in the top right panel. The bottomleft panel shows the power spectrum at fixed redshift ac-cording to the Limber ratio Eq. (36), sweeping through (cid:96) . Each segment in the top right panel of Fig. 2 crosses thelines in the bottom left panel. We illustrate this furtherin Fig. 3. III. SFB DECOMPOSITION
This section describes our SFB decomposition for agalaxy survey with mask and selection function. Welargely follow Samushia (2019) for the radial basis func-tions and Leistedt et al. (2012) for the angular/radialsplit in the estimator.We start by giving a review of the basis functions, thenwe add window and selection functions, we model thediscrete galaxy distribution, and estimate the covariancematrix.
A. Spherical Fourier-Bessel basis with potentialboundary conditions
As noted in the introduction, we choose the eigenbasisof the Laplacian because that makes scalar modes inde-pendent to first order in cosmological perturbation theory(Bernardeau et al. f becomes ∇ f “ r BB r ˆ r B f B r ˙ ` r sin θ BB θ ˆ sin θ B f B θ ˙ ` r sin θ B f B φ , (38)where r , θ , and φ are the comoving angular diameterdistance, zenith angle, and azimuthal angle, respectively(see Appendix B for a derivation). The eigenbasis toEq. (38) that satisfies ∇ f “ ´ k f (39)for some mode k is of the form (see e.g. Samushia 2019) f (cid:96)µ p k ; r, θ, φ q “ “ c j j (cid:96) p kr q ` c y y (cid:96) p kr q ‰ ˆ “ c p P µ(cid:96) p cos θ q ` c q Q µ(cid:96) p cos θ q ‰ ˆ “ c ` e iµφ ` c ´ e ´ iµφ ‰ , (40)where the c i are constants, and j (cid:96) and y (cid:96) are sphericalBessel functions of the first and second kind, and P µ(cid:96) and Q µ(cid:96) are Legendre functions of the first and second kind.The constants c i are set by boundary conditions. First,the spherical Bessel of the second kind, y (cid:96) diverges withvanishing argument; hence, typically c y “ . Typi-cally, the functions also need to be periodic about theazimuthal angle φ ; therefore, µ “ , , , . . . is an integer.
500 600 700 800 900 1000 r in Mpc /h − . − . . . . g n ‘ ( r ) i n [ h / M p c ] × − n = 1 n = 2 n = 3 n = 4 ‘ = 0 ‘ = 1 ‘ = 2 ‘ = 3 r min in Mpc /h . . . . . . k n ‘ i n h / M p c ‘ = 10 FIG. 4. Left: The radial basis functions for potential boundary conditions as a function of r . Color indicates the modes n ,line style indicates (cid:96) as shown in the legend. Here, r min “ h ´ Mpc and r max “ h ´ Mpc . Right: k n(cid:96) for potentialboundary conditions as a function of r min when (cid:96) “ . The grey lines are for r min “ , and we fix r max “ h ´ Mpc . Then, the functions also need to be finite for cos θ “ ˘ ,typically; therefore, c q “ , and (cid:96) “ , , , . . . is an inte-ger, and ´ (cid:96) ď µ ď (cid:96) .Effectively, the preceding paragraph imposed bound-ary conditions at r min “ and assumed coverage of thewhole sky. Typically (e.g., Heavens and Taylor 1995,Tadros et al. et al. et al. r max such that the survey volume is con-tained within a sphere of radius r max . This restricts theSFB volume, i.e., the volume on which the SFB transformis performed, as the region from ď r ď r max . Demand-ing the basis functions to be orthogonal then leads to adiscrete spectrum of modes k “ k n(cid:96) .Realistic galaxy surveys do not occupy the entire SFBvolume, but are restricted in both redshift and angulararea, and, therefore, they leave large fractions of the SFBvolume unobserved. This leads to the deconvolution ofthe window function to be numerically unstable. It alsoresults in wasted computational resources if the surveycovers only a (potentially thick) shell at high redshift.The analogous picture for a standard Fourier transformwould be to have a transform box that is much larger thanthe survey volume. Therefore, the selection function willvanish for part of the SFB volume, and, because in thatcase some modes are not well constrained, the inversionof the window function becomes numerically unstable.In this paper, we employ two strategies to deal withthis problem. First, we follow Hivon et al. (2002), Alonso et al. (2019) and bin the pseudo SFB power spectrum intobandpowers. This combines several poorly-constrainedmodes into one well-constrained mode. We rely on thisstrategy especially for the angular mask so that we canleverage the full-sky spherical harmonic algorithms fromthe HealPy software (Górski et al. et al. r min so that theSFB volume extends from r min ď r ď r max . For galaxysurveys that start at some minimum redshift this elimi-nates from the SFB tranform volume a hole around theorigin. As a result, the inversion of the window func-tion is numerically well behaved even without resortingto bandpower binning. Furthermore, the number of SFBmodes is reduced not just by the boundary condition at r max , but the boundary condition at r min also reducesthe number of modes further by the fraction r { r .In all cases considered in this paper, this eliminates theneed for bandpower binning in the radial direction.We differ from Samushia (2019) in that we use poten-tial boundary conditions (Fisher et al. k n(cid:96) , as shownin Appendix C. In the appendix we also derive that theradial basis functions with such boundary conditions be-come a linear combination of spherical Bessels of the firstand second kind, g n(cid:96) p r q “ c n(cid:96) j (cid:96) p k n(cid:96) r q ` d n(cid:96) y (cid:96) p k n(cid:96) r q , (41)which satisfy an orthonormality relation ż r max r min d r r g n(cid:96) p r q g n (cid:96) p r q “ δ Knn , (42)where δ Knn is a Kronecker delta, and the coefficients c n(cid:96) and d n(cid:96) are derived in Appendix C. With Eq. (42), theFourier pair Eqs. (1) and (2) remains a Fourier pair withthe discrete k n(cid:96) -spectrum, and the pair becomes δ p r q “ ÿ n(cid:96)m ” g n(cid:96) p r q Y (cid:96)m p ˆ r q ı δ n(cid:96)m , (43) δ n(cid:96)m “ ż d r ” g n(cid:96) p r q Y ˚ (cid:96)m p ˆ r q ı δ p r q , (44)where δ n(cid:96)m “ δ (cid:96)m p k n(cid:96) q . Note that our choice to normal-ize g n(cid:96) p r q as in Eq. (42) changes the units of δ n(cid:96)m com-pared to Eqs. (1) and (2). In effect, this choice of unitstakes into account the survey volume at this stage ratherthan at the stage of forming the correlation function.Examples of the resulting basis functions and modes k n(cid:96) are show in Fig. 4. We point out that the (cid:96) “ modes are closely related to taking the average of thetransformed field δ p r q . Also, a larger r min results in asmaller volume, and, therefore, fewer modes that can beconstrained. B. Window and selection function
The observed number density n p r q of galaxies is sub-ject to the window and selection function W p r q of thesurvey, which we define as the fraction of galaxies ob-served at position r . For a random catalogue subject tothe same window function, with density n r p r q , and with { α as many galaxies as the survey, we then have α x n r p r qy “ W p r q ¯ n , (45)where x n r p r qy “ α ´ ¯ n p r q is the average number densityof the ensemble of random catalogs, and ¯ n is the aver-age number density in the survey. Note that Eq. (45)can equivalently be expressed in terms of the limit lim α Ñ α n r p r q “ ¯ n p r q “ W p r q ¯ n . With this definitionof the window function, we define the effective volume as V eff “ ż d r W p r q , (46)so that the average number density ¯ n becomes ¯ n “ N obsgal V eff , (47)and N obsgal is the observed number of galaxies in the sur-vey. Any variation across the survey in the actual averagenumber density, e.g., due to an evolving luminosity func-tion, is absorbed into W p r q . Our treatment is in linewith Taruya et al. (2021), and our W p r q takes the role ofthe function G p r q in Feldman et al. (1994), except thatwe do not at present include a weighting scheme.In a sense, there are two window functions here: first,the one defined by the SFB procedure and limited by r min ď r ď r max , and second, W p r q , which defines thegeometry and selection of the survey. However, the firstone should be irrelevant as long as the survey volume isentirely inside r min ď r ď r max and as along as sufficientnumber of modes are included in the SFB analysis. The observed density fluctuation field is, then, δ obs p r q “ n p r q ´ α n r p r q ¯ n “ n p r q ¯ n ´ W p r q , (48)where Eq. (45) was used in the limit that the randomcatalogue has an infinite number of galaxies, or α Ñ .Because the observed density n p r q is also subject to thewindow function W p r q , the observed and true densitycontrasts are related by δ obs p r q “ W p r q δ A p r q , (49)where we attach the superscript ‘A’ to refer to the localaverage effect (see Section III H below). Transformingto SFB-space and expressing δ A p r q in terms of its SFBdecomposition Eqs. (43) and (44), we get δ obs n(cid:96)m “ ÿ n (cid:96) m W n (cid:96) m n(cid:96)m δ An (cid:96) m , (50)where W n (cid:96) m n(cid:96)m “ ż d r r g n(cid:96) p r q g n (cid:96) p r qˆ ż d ˆ r Y ˚ (cid:96)m p ˆ r q Y (cid:96) m p ˆ r q W p r, ˆ r q . (51)
1. Properties and implementation
From Eq. (51) follows the symmetry W n ,(cid:96) , ´ m n,(cid:96),m “ p´ q m ` m W n ,(cid:96) ,m , ˚ n,(cid:96), ´ m , (52)and the Hermitian property W n (cid:96) m n(cid:96)m “ W n(cid:96)m, ˚ n (cid:96) m . (53)In the special case that W p r q “ everywhere, W n (cid:96) m n(cid:96)m “ δ Knn δ K(cid:96)(cid:96) δ Kmm , (54)which follows from Eq. (42) and Eq. (A7).In all generality, Eq. (51) can be simplified for compu-tational convenience by expressing the window functionin terms of an angular transform. That is, introduce W LM p r q “ ż d ˆ r Y ˚ LM p ˆ r q W p r, ˆ r q . (55)Then, W n (cid:96) m n(cid:96)m “ p´ q m ÿ L G (cid:96)(cid:96) L ´ m,m ,m ´ m ˆ ż d r r g n(cid:96) p r q g n (cid:96) p r q W L,m ´ m p r q , (56)where we used Eq. (A12) and introduced the Gaunt fac-tor Eq. (A14). In writing Eq. (56) we performed the an-gular transform of the window function only as that leadsto a computationally suitable form. Had we performed afull SFB transform, we would have been left with an infi-nite sum over n that converges slowly, in addition to theneed of computing integrals over three spherical Besselfunctions.0 C. Discrete points
Now we specialize the SFB decomposition to the casethat we have galaxies represented by discrete points.That is, we assume the number density is given by n p r q “ ÿ p δ D p r ´ r p q , (57)where the sum is over all points (galaxies) in the survey.In the 3DEX approach (Leistedt et al. δ obs n(cid:96)m “ ż d Ω ˆ r Y ˚ (cid:96)m p ˆ r q δ obs n(cid:96) p ˆ r q , (58)where δ obs n(cid:96) p ˆ r q “ ż r max r min d r r g n(cid:96) p r q δ obs p r, ˆ r q (59)represents an angular field for each n and (cid:96) , and g n(cid:96) isdefined in Eq. (41). For the density contrast Eq. (48)with number density Eq. (57), δ obs n(cid:96) p ˆ r q “ n ÿ p δ D p ˆ r ´ ˆ r p q g n(cid:96) p r p q ´ W n(cid:96) p ˆ r q , (60)where W n(cid:96) p ˆ r q “ ż r max r min d r r g n(cid:96) p r q W p r q . (61)Eq. (60) is an exact expression for the observed densitycontrast δ obs n(cid:96) p ˆ r q . However, for the angular transform wewish to make use of the fast HEALPix scheme , , andwe need the density contrast in pixel i averaged over thepixel area ∆Ω i , ¯ δ obs n(cid:96) p ˆ r i q “ i ż ∆Ω i dΩ ˆ r δ obs n(cid:96) p ˆ r q (62) “ n ∆Ω i ÿ p P ∆Ω i g n(cid:96) p r p q ´ W n(cid:96) p ˆ r i q , (63)where we assumed that W n(cid:96) p ˆ r q varies slowly over thesize of an angular pixel. Then, the angular transform isperformed: ¯ δ obs n(cid:96)m “ ÿ i ∆Ω i Y ˚ (cid:96)m p ˆ r i q ¯ δ obs n(cid:96) p ˆ r i q . (64)In what follows we will generally drop the bar indicatingthe angular-pixel-averaging. https://healpix.jpl.nasa.gov/index.shtml https://healpy.readthedocs.io/en/latest/ D. Power spectrum estimation
Hivon et al. (2002), Wandelt et al. (2001) use a pseudo- C (cid:96) method to estimate the power spectrum. Translat-ing to the SFB decomposition, the pseudo- C (cid:96) methodassumes that much of the information about the powerspectrum is contained in the pseudo-power spectrum ˆ C obs (cid:96)nn “ (cid:96) ` ÿ m δ obs n(cid:96)m δ obs , ˚ n (cid:96)m . (65)That is, we ignore off-diagonal terms L ‰ (cid:96) and M ‰ m , and average over m . The effect of the window isthen described by a mixing matrix between the ˆ C obs (cid:96)nn and ˆ C A(cid:96)nn , ˆ C obs (cid:96)nn “ ÿ LNN M LNN (cid:96)nn ˆ C ALNN , (66)where we used Eq. (50) and defined M LNN (cid:96)nn “ (cid:96) ` ÿ mM W NLMn(cid:96)m W N LM, ˚ n (cid:96)m , (67)and the index ‘A’ on C A(cid:96)nn indicates the local averageeffect, see Section III H. Next, with Eqs. (55) and (56)we get M LNN (cid:96)nn “ L ` π ÿ L ˆ (cid:96) L L ˙ ÿ M ˆ ż d r r g n(cid:96) p r q g NL p r q W L M p r qˆ ż d r r g n (cid:96) p r q g N L p r q W ˚ L M p r q , (68)and we used the orthogonality of the Gaunt factorEqs. (A15) and (A16). (The sum over M could be per-formed first. However, that approach is much more mem-ory intensive, so that computing the integrals first endsup being faster. We have also avoided expressing theresult in terms of a full SFB transform, as that wouldrequire a slowly-converging sum over n .) Note that thematrix p L ` q ´ M LNN (cid:96)nn is symmetric, but M by itselfis not.
1. Separable mask and radial selection
It is quite common that the window function is sepa-rable into a radial and an angular term, W p r q “ φ p r q M p ˆ r q . (69)If the flux limit in a blind survey is near L ˚ , then the se-lection could change dramatically as a function of angulardepth variations that are due to, e.g., atmospheric varia-tions, and the separation of angular and radial selectionwould be a poor approximation. However, eBOSS, for1example, had more targets selected in regions where twoor more plates overlapped (e.g. de Mattia et al. et al. W LM p r q “ φ p r q W LM , (70)where W LM “ ż d ˆ r Y ˚ LM p ˆ r q M p ˆ r q , (71)and Eq. (68) becomes M LNN (cid:96)nn “ L ` π ÿ L ˆ (cid:96) L L ˙ ÿ M | W L M | ˆ ż d r r g n(cid:96) p r q g NL p r q φ p r qˆ ż d r r g n (cid:96) p r q g N L p r q φ p r q , (72)which is also separable, and therefore significantly re-duces computation cost. Eq. (56) simplifies in a similarmanner.In the special case that W p r q “ everywhere, we re-cover the unit matrix M LNN (cid:96)nn “ δ K(cid:96)L δ KnN δ Kn N , (73)as expected. E. Shot noise
The sampling of the density field by a limited num-ber of points leads to a shot noise component in thepower spectrum. To estimate the shot noise, we startwith (Feldman et al. @ n p r q n p r q D “ ¯ n p r q ¯ n p r q “ ` ξ p r , r q ‰ ` ¯ n p r q δ D p r ´ r q , (74) @ n p r q n r p r q D “ α ´ ¯ n p r q ¯ n p r q , (75) @ n r p r q n r p r q D “ α ´ ¯ n p r q ¯ n p r q ` α ´ ¯ n p r q δ D p r ´ r q . (76)The density contrast is given by Eq. (48), and the en-semble average becomes @ δ obs p r q δ obs p r q D “ W p r q W p r q ξ p r , r q` p ` α q W p r q δ D p r ´ r q ¯ n , (77)where we used Eq. (45). Therefore, the SFB transformof the shot noise term becomes (see Eq. (44)) N obs “ n SFB r W p r q δ D p r ´ r qs (78) “ n W n (cid:96) m n(cid:96)m , (79) in the limit α Ñ , and the W matrix is defined inEq. (51). The window-corrected shot noise, therefore,is, in matrix form, W ´ { ¯ n .For the pseudo-SFB-power-spectrum estimator theshot noise simplifies significantly. Averaging over themodes m “ m and assuming (cid:96) “ (cid:96) , Eq. (79) becomes N obs (cid:96)nn “ n ? π ż d r r g n(cid:96) p r q g n (cid:96) p r q W p r q , (80)where we used Eq. (56). Eq. (80) can be implementedvery efficiently. F. Pixel window
The pixel window refers to a distortion of the powerspectrum due to binning galaxies into pixels. In the ra-dial direction, we do not bin the galaxies, see Eq. (60),and, therefore, we do not have a radial pixel window(Leistedt et al. pixwin functionof
HealPy to correct for the pixel window. We confirmthe accuracy of this procedure with simulations in Sec-tion IV.
G. Bandpowers
We use a similar approach as Hivon et al. (2002),Alonso et al. (2019) to bin the SFB power spectruminto bandpowers. This is necessary as the mixing ma-trices in Eqs. (51) and (66) are, in general, not invertiblewith finite-precision arithmetic. Compared to those au-thors our situation is complicated, but not significantlychanged, by the fact that we may need to bin not onlyin (cid:96) , but also in the k -modes n and n .We define the bandpower-binned pseudo- C (cid:96) SFBpower spectrum as a weighted sum over modes, ˆ B obs LNN “ ÿ (cid:96)nn r w (cid:96)nn LNN ˆ C obs (cid:96)nn , (81)where r w (cid:96)nn LNN is typically a rectangular sparse matrixthat takes the average of neighboring modes p (cid:96)nn q „p LN N q . The operation Eq. (81) is a type of compres-sion, where the compression matrix r w must satisfy thenormalization ÿ (cid:96)nn r w (cid:96)nn LNN “ . (82)In matrix notation, we write the compression operationEq. (81) and the corresponding decompression operation B W “ r wC W , C W » r vB W , (83) B “ wC , C » vB , (84)2where r w and r v are rectangular matrices operating onwindow-convolved power spectra, and w and v are rect-angular matrices operating on cleaned power spectra. Weuse the index ‘W’ to indicate that we are only consideringthe window convolution. That is, C W “ M C , B W “ N B , (85)where M is given by Eq. (67), and we can use the firstof Eq. (83), the first of Eq. (85), and the last of Eq. (84)to get N “ r w M v . (86)The compression matrix w is obtained by inverting thesecond of Eq. (85), using the first of Eq. (83), and thefirst of Eq. (85) to get B “ N ´ B W “ N ´ r wC W “ N ´ r w M C , (87)or (Alonso et al. w “ N ´ r w M . (88)Similarly, we find r v “ M v N ´ . (89)Eq. (86) then implies that wv “ r w r v “ I . (90)Eq. (90) is equivalent to assuming that a decompression-then-recompression cycle is lossless. That is, the com-pressed representation is unaffected by decompression.The opposite of compression-then-decompression, how-ever, will in general incur losses in the compression step,so that vw ‰ I except in special cases.Furthermore, once the information is lost, repeatedcompression-then-decompression cycles do not lose moreinformation. That is, p vw q n “ vw for integer n ě .Note that w and r v must be calculated via Eqs. (88)and (89), because in the general case we have p vw q : ‰ vw , and, therefore, they are not unique in satisfyingEq. (90). That is, they are not the unique Moore-Penroseinverses of the matrices v and r w (Dresden 1920, Penrose1955).Since w and r v can be expressed in terms of r w , v ,and the window mixing matrix, our procedure consistsof choosing a compression matrix r w and a decompres-sion matrix v .How to choose r w and v ? For r w we have already sug-gested that its operation on a power spectrum shall weighneighboring modes equally and satisfy the normalizationEq. (82). Since the modes n and n refer to modes k n(cid:96) ,their spacing is not independent, and we need to bin fourmodes, nine modes, or similar at a time. In this paper,our binning strategy is completely specified by the twonumbers ∆ (cid:96) and ∆ n . A natural choice for the decompression v is, then, asthe Moore-Penrose inverse of r w . Indeed, for the afore-mentioned choice for r w that takes the average of neigh-boring modes, this means that v is a step function thatassigns the same value or a value proportional to (cid:96) p (cid:96) ` q to all modes within a bin (Hivon et al. et al. H. Local Average Effect
In this section we recognize that the average numberdensity ¯ n in Eq. (48) must in practice be measured fromthe survey itself. This is often called the integral con-straint (Beutler et al. local average effect (de Putter et al. et al. ¯ n “ ` ` ¯ δ ˘ ¯ n true , (91)where ¯ n true is the underlying density contrast in thewhole universe, and the average density contrast in thesurvey volume is ¯ δ “ V eff ż d r W p r q δ p r q , (92)with the effective volume defined in Eq. (46). Therefore,with our model in Eq. (48), the measured density contrastis (Taruya et al. δ obs p r q ” δ W,A p r q “ W p r q δ p r q ´ ¯ δ ` ¯ δ , (93)where δ p r q is the true density contrast, the superscript‘A’ refers to the local average effect, and the superscript‘W’ refers to the effect of the window convolution.The SFB transform of Eq. (93) is δ W,An(cid:96)m “ ÿ n (cid:96) m W n (cid:96) m n(cid:96)m δ n (cid:96) m ´ d n (cid:96) m ¯ δ ` ¯ δ , (94)where we used Eqs. (44) and (50), and we defined d n (cid:96) m “ ? π δ K(cid:96) δ Km ż d r r g n p r q . (95)Using Eqs. (43), (44) and (50), Eq. (92) can be written ¯ δ “ V eff ÿ n (cid:96) m d W, ˚ n (cid:96) m δ n (cid:96) m , (96)3where we used Eq. (53) to define d Wn (cid:96) m as d Wn (cid:96) m “ ÿ n(cid:96)m W n(cid:96)mn (cid:96) m d n(cid:96)m . (97)Then, expanding Eq. (94) for small ¯ δ we get δ W,An(cid:96)m “ ÿ n (cid:96) m W n (cid:96) m n(cid:96)m ” δ n (cid:96) m ´ d n (cid:96) m ¯ δ ´ δ n (cid:96) m ¯ δ ` d n (cid:96) m ¯ δ ` δ n (cid:96) m ¯ δ ` O p ¯ δ q ı . (98)Since d n(cid:96)m „ ? V and ¯ δ „ {? V we expand in the vol-ume V . We also assume C (cid:96)nn ! V eff . Then, the correla-tion function becomes A δ W,ANLM δ W,A, ˚ N L M E “ ÿ n(cid:96)m W n(cid:96)mNLM ÿ n (cid:96) m W n (cid:96) m , ˚ N L M ˆ ” x δ n(cid:96)m δ ˚ n (cid:96) m y ´ @ δ n(cid:96)m ¯ δ D d n (cid:96) m ´ d n(cid:96)m @ δ ˚ n (cid:96) m ¯ δ D ` d n(cid:96)m d n (cid:96) m @ ¯ δ D ` O ´ V ´ ¯ ı , (99)where we assume the field δ to be Gaussian. The 2-point terms entering this expression are x δ n(cid:96)m δ ˚ n (cid:96) m y “ δ K(cid:96)(cid:96) δ Kmm C (cid:96)nn ` n p W ´ q n (cid:96) m n(cid:96)m , (100) @ δ n(cid:96)m ¯ δ D “ V eff ÿ n d Wn (cid:96)m C (cid:96)nn ` nV eff d n(cid:96)m , (101) @ ¯ δ D “ V ÿ (cid:96) n n D W(cid:96) n n C (cid:96) n n ` nV ÿ n (cid:96) m d n (cid:96) m d Wn (cid:96) m , (102)where we included the shot noise term Eq. (79), and wedefined the unnormalized power spectrum of a constantfield D W(cid:96) n n “ ÿ m d Wn (cid:96) m d W, ˚ n (cid:96) m (103) “ p (cid:96) ` q ÿ n ÿ n M n n (cid:96) n n d n d n , (104)and Eq. (67) was used in the last line. Now the correla-tion function becomes A δ An(cid:96)m δ A, ˚ n (cid:96) m E “ δ K(cid:96)(cid:96) δ Kmm C (cid:96)nn ` n p W ´ q n (cid:96) m n(cid:96)m ´ d n (cid:96) m V eff ÿ n d Wn (cid:96)m C (cid:96)nn ´ d n(cid:96)m V eff ÿ n d W, ˚ n (cid:96) m C (cid:96) n n ` „ nV eff ` @ ¯ δ D d n(cid:96)m d n (cid:96) m ` O ´ V ´ ¯ , (105)where we corrected for the window function. Only the first and fifth terms are proportional to δ K(cid:96)(cid:96) δ Kmm , and so wecannot take the pseudo- C (cid:96)nn power spectrum at this stage. To do so, we now add back the two window functionsin Eq. (99) to get a prediction for the observed pseudo power spectrum C W,ALNN “ ÿ (cid:96)nn M (cid:96)nn LNN C (cid:96)nn ` N obs LNN ´ V eff L ` ÿ M d W, ˚ N LM ÿ n(cid:96)m W n(cid:96)mNLM ÿ n d Wn (cid:96)m C (cid:96)nn ´ V eff L ` ÿ M d WNLM ÿ n (cid:96) m W n (cid:96) m , ˚ N LM ÿ n d W, ˚ n (cid:96) m C (cid:96) n n ` ÿ (cid:96)nn M (cid:96)nn LNN δ K(cid:96) d n d n „ nV eff ` @ ¯ δ D ` O ´ V ´ ¯ , (106)where we used Eq. (80). The third and fourth terms are the same except for N Ø N . To simplify these two terms,we express them in terms of chains of window functions W k that we define in Eq. (E1) and study in Appendix E. We If d Wn(cid:96)m δ K(cid:96) δ Km would be a good approximation, which is thecase in the absence of a window function, then the pseudo- C (cid:96) ap-proach works well. This is an argument for using eigenfunctions tailored to the survey geometry. That is, if the window effectsare already captured by the eigenfunctions, then the calculationhere would simplify dramatically. ÿ M d W, ˚ N LM ÿ n(cid:96)m W n(cid:96)mNLM ÿ n d Wn (cid:96)m C (cid:96)nn “ ÿ M ÿ n (cid:96) m ÿ n(cid:96)m ÿ n ÿ n (cid:96) m W n (cid:96) m , ˚ N LM W n(cid:96)mNLM W n (cid:96) m n (cid:96)m d n (cid:96) m d n (cid:96) m C (cid:96)nn “ ÿ (cid:96)nn ÿ (cid:96) n n ÿ Mmm W N LMn (cid:96) m W n(cid:96)mNLM W n (cid:96) m n (cid:96)m δ K(cid:96) d n d n C (cid:96)nn “ ÿ (cid:96) n n δ K(cid:96) d n d n ÿ (cid:96)nn C (cid:96)nn W ¨˝ L (cid:96) (cid:96) N n n N n n ˛‚ , (107)where W is defined in Eq. (E1). Eq. (106) now becomes C W,ALNN “ ÿ (cid:96)nn M (cid:96)nn LNN „ C (cid:96)nn ` δ K(cid:96) d n d n ˆ nV eff ` @ ¯ δ D˙ ` N obs LNN ´ V eff L ` ÿ (cid:96) n n δ K(cid:96) d n d n ÿ (cid:96)nn C (cid:96)nn »– W ¨˝ L (cid:96) (cid:96) N n n N n n ˛‚ ` @ N Ø N Dfifl . (108)The window de-convolved power spectrum is then C A “ M ´ C W,A .In the absence of a window function and assuming C (cid:96)nn δ Knn , as well as using Eq. (95), we get C A(cid:96)nn » ˆ ` A C V eff ˙ C (cid:96)nn ´ δ K(cid:96) d n d n V eff B nn , (109)where we included further terms from the expansionEq. (99), and we defined A “ ÿ n d n V eff C n n C , (110) B nn “ C nn ` C n n ´ A C ´ C nn C n n V eff . (111)To a good approximation A » , and the last term in B nn can be neglected if the effective volume is sufficientlylarge and the shot noise sufficiently low. Furthermore,a good approximation is d n δ Kn . Therefore, only the p (cid:96), n, n q “ p , , q mode is significantly affected. That is,the main effect of super-sample variance on the measuredpower spectrum is to reduce the power in the largestmode. I. Covariance matrix of power spectrum
In this section we provide a covariance matrix for theSFB power spectrum. Several approaches have been usedpreviously. Percival et al. (2004), Wang et al. (2020) tracethe likelihood function either on a grid or using MarkovChain Monte Carlo techniques. Wang et al. (2020), e.g.,use simulations to measure the covariance matrix fromsuites of mock catalogues. An analytical approach for the3D power spectrum multipoles is presented in Wadekar and Scoccimarro (2020). In this paper, we get an ana-lytical estimate for the SFB power spectrum assumingthat the density contrast is Gaussian, and we compare to100 log-normal simulations. Non-Gaussian terms in theform of the disconnected trispectrum could be includedsimilarly to Taruya et al. (2021), Sugiyama et al. (2020).Super-sample variance (e.g., de Putter et al. et al.
Beat coupling is modemixing due to the window function with correlation be-tween pairs of non-linear modes and one large mode. The local average effect is due to the large-scale mode modu-lating the average number density inside the survey vol-ume. Both these effects can be treated in the manner ofSection III H and Eq. (93).The covariance matrix on the observed SFB powerspectrum is V LNN , obs (cid:96)nn ” A ˆ C obs (cid:96)nn ˆ C obs LNN E ´ C obs (cid:96)nn C obs LNN “ p (cid:96) ` qp L ` q ÿ mM ” A δ W,An(cid:96)m δ W,ANLM
E A δ W,A, ˚ n (cid:96)m δ W,A, ˚ N LM E ` A δ W,An(cid:96)m δ W,A, ˚ N LM E A δ W,ANLM δ W,A, ˚ n (cid:96)m E ı , (112)where we used Wick’s theorem for a Gaussian densitycontrast. We simplify Eq. (112) in Appendix D. How-ever, an analytical calculation remains computationallyexpensive.To get the covariance matrix for the window-correctedpower spectrum, we write the matrix equation V “ N ´ V obs N ´ ,T . (113)where N is the bandpower-binned window coupling ma-trix given in Eq. (86), and the binning of the covariancematrix is implied.5 FIG. 5. The panels show the covariance matrix V defined in Eq. (113). Top row: V A as measured from 100 simulations,including the local average effect. Bottom row: The analytical prediction, without local average effect. Left column: Shot noiseonly. Right column: Linear power spectrum with shot noise. Here we use a simulation with
50 % sky coverage and bandpowerbinning with ∆ (cid:96) “ . The order of the p (cid:96)nn q indices is such that each block of ten indices is for one (cid:96) -bin starting with (cid:96) “ , for the block starting with index 0 and ending with (cid:96) “ , for the block starting at index 40. Within each block n “ n increases from one to 10. Note that the nonlinear color scheme amplifies small elements. A reasonably precise estimate can be obtained bycounting modes and assuming the covariance matrix isdiagonal. That is, V LNN (cid:96)nn » δ K(cid:96)L N modes “ C binned (cid:96)nN C binned Ln N ` C binned (cid:96)nN C binned Ln N ‰ , (114)where the power spectrum includes the shot noise, C binned (cid:96)nn “ C signal (cid:96)nn ` N shot (cid:96)nn , and N modes “ f vol p (cid:96) ` q ∆ (cid:96) ∆ n , (115)where ∆ (cid:96) and ∆ n are the bin widths for modes k n(cid:96) , and f vol is the fraction of the SFB transform-volume that isoccupied by the survey, defined by f vol ” V SFB ż d r τ r W p r q ´ W threshold s , (116) where τ p x q is a step function and W threshold is a thresholdof the window function.The shot noise takes into account the variation ofthe number density across the survey, and it enters inEq. (114) as part of the power spectrum. The incompletevolume coverage enters as a reduction in the number ofmodes, and it is needed for the stability of the window-deconvolution when there are large unobserved regions inthe SFB volume.In Fig. 5 we show the covariance matrices for a set ofsimulations that contain only shot noise (top left) as wellas for a set of simulations with a physical galaxy powerspectrum with bias b “ . at effective redshift z eff “ (top right). In the figure we also show the analyticalresult from Eq. (112) (bottom panels).The colorbar in the figure is nonlinear. As a result,small elements appear amplified. To provide a more use-6 FIG. 6. Relative comparison between the covariance matrices as in Eq. (117). This avoids amplifying small deviations farfrom the center diagonal, and it shows that the analytic result largely agrees with the simulations. The very largest mode isset to zero, because it is affected by the local average effect that is not included in the analytical result (see Fig. 7). ‘nn ) index − . − . − . . . . R e l a t i v e d i ff e r e n ce ∆ ρ = V s i m u l a t i o n − V t h e o r y d i ag ( V t h e o r y ) Window-corrected covariance matrix, shot noise only ‘ = L‘ = L + 1 ‘ = L , mode counting ‘ = L + 1, mode counting 0 10 20 30 40 50( ‘nn ) index − . − . − . − . . . . R e l a t i v e d i ff e r e n ce ∆ ρ = V s i m u l a t i o n − V t h e o r y d i ag ( V t h e o r y ) Window-corrected covariance matrix ‘ = L‘ = L + 1 ‘ = L , mode counting ‘ = L + 1, mode counting FIG. 7. Comparison of the diagonal and the (cid:96) “ L ` off-diagonal window-corrected covariance between simulations andanalytical result. Left shows the shot noise only, and right including shot noise and a non-zero power spectrum. In theanalytical result we do not include local average effect. Thus, the first mode in the simulations is suppressed compared to theanalytical result. The black dashed and dash-dotted lines show the approximate result from mode-counting Eq. (114). ful comparison, we introduce the difference between twocovariance matrices, scaled to the center diagonal. Thatis, we introduce the relative difference ∆ ρ ij “ C Aij ´ C Bij b C Bii C Bjj , (117)and we choose C B to refer to the analytic result. ∆ ρ does not suffer from amplification of small differences farfrom the diagonal.Therefore, in Fig. 6 we show the relative difference be-tween the covariance matrix as obtained from simulationsand the analytical result. However, in the figure we re-move the largest mode, since we have not included the lo- cal average effect in the analytical calculation. All othermodes are statistically essentially equal between simula-tion result and analytics.To show this more clearly, we present Fig. 7, wherewe compare the main diagonal and the (cid:96) “ L ` diag-onal of the covariance matrices using the same statisticEq. (117). Within the noise, we find good agreementbetween simulations and analytical result. J. Performance scaling
In this section we give a brief overview of the perfor-mance behavior of the code. We consider the scaling7 ‘ n r min = 500 h − Mpc, r max = 1000 h − Mpc k max = 0 . h Mpc − k max = 0 . h Mpc − k max = 0 . h Mpc − FIG. 8. Here we show the modes that are required to achievea given k max for the given radial boundaries. The solid blackcurves are constant- k contours for r min “ . of the SFB decomposition, the coupling matrix, and thenumber of modes with the parameters of the SFB powerspectrum estimation.Typically, it is desirable to calculate all modes upto some k max . The total number of modes can beestimated in the same way as for a standard Fouriertransform by estimating the fundamental frequency fromthe volume that is being transformed, i.e., V SFB » π ` r ´ r ˘ and k F » π { V { . Then, thetotal number of modes is approximately N modes » π k { k F » π k ` r ´ r ˘ . These are themodes that need to be calculated for transform.However, as shown in Fig. 8, the boundary at r min changes the structure about which modes need to be cal-culated. In the figure, modes with k ď k max are in theshaded region when, for illustration, r min “ h ´ Mpc ,and all the modes below the solid black curves need tobe calculated if r min “ . The figure shows that fewerlow- (cid:96) modes are needed when r min is finite. However,most modes are at large (cid:96) , and, therefore, this is only asmall computational reduction.The algorithm now scales as follows. First, the sumin Eq. (60) is performed, and then for each p n, (cid:96) q -combination the spherical harmonic transform Eq. (64)is performed. Hence, the execution time of the transformscales as T „ O ” n max (cid:96) max ´ N gal ` N Healpix (cid:96) ¯ı , (118)where N Healpix (cid:96) „ N { is the number of operationsneeded for the spherical harmonic transform, and N pix is the number of HEALPix pixels.To estimate N pix “ n , we need to estimate n side ,which we do by considering the angular resolution. Rec-ommended is (cid:96) max » n side . However, our algorithmsdealing with the window function will need to go to L “ (cid:96) max . Hence, we estimate n side “ ceil p log p (cid:96) max ` qq , (119) https://healpix.jpl.nasa.gov/ where ceil p x q is the smallest integer greater than x .Eq. (119) guarantees that n side is a power of two.The angular resolution is determined by (cid:96) max , whichis determined by k max and r max by Limber’s relationEq. (36). However, we note that the actual numberneeded for (cid:96) max tends to be smaller by a few percent. IV. FUTURE APPLICATIONS
In this section we will test the SFB estimator pre-sented in this paper for several use cases. First, wewill consider simplistic simulations of surveys similar tothe High-Latitude Spectroscopic Survey (HLSS) of the
Nancy Grace Roman Space Telescope which will benefitfrom large radial modes, and we will consider
SPHEREx and
Euclid for wide-angle surveys.We bin into bandpowers by selecting ∆ (cid:96) » f sky , (120)and then we round to the nearest integer. Eq. (115) thensuggests ∆ n » f sky f vol . (121)For all cases in this paper, this results in ∆ n “ .To do the window deconvolution, it is important to en-sure that all modes are complete. For the angular powerspectrum, Alonso et al. (2019), Leistedt et al. (2013) sug-gest estimating up to (cid:96) max , and then discarding all themodes above (cid:96) max . Wang et al. (2020) argue that (inour notation) the sum Eq. (66) only converges with theinclusion of modes k n(cid:96) ą k max , and they do numericalexperiments to estimate the maximum k needed.We take a similar approach, which is demonstrated inFig. 9, where in the left panel the relative contribution ofwindow-convolved modes to a physical mode near k max is shown. That is, we plot the relative contributions ofall observed modes p (cid:96)nn q that contribute to the physi-cal mode p LN N q using Eq. (66). Assuming a flat powerspectrum and summing the absolute values of the cou-pling matrix M ´ , then, allows us to estimate the contri-bution from all modes above some k large . This is shownin the right panel of Fig. 9.Next, we iteratively increase k large until the contri-bution from k n(cid:96) ą k large to the most affected mode k NL ď k max is less than . The results for k large arethe dashed vertical lines in either panel of Fig. 9.That is, by including all modes up to k large ą k max inthe SFB power spectrum estimation, we get reasonableconfidence that all modes k ă k max can be fully decon-volved by the inversion of Eq. (66).8 .
04 0 .
05 0 . k n‘ in h/ Mpc0 . . . . ( M − ) ‘ nn L NN / P ‘ nn ( M − ) ‘ nn L NN Observed modes contributing to k NL ∼ . h/ Mpc
RomanSPHERExEuclid 0 .
040 0 .
045 0 .
050 0 .
055 0 .
060 0 . k NL in h/ Mpc10 − − − F r a c t i o n o f k > k l a r g e c o n tr i bu t i o n s k n‘ > k large contribution to k NL RomanSPHERExEuclid
FIG. 9. This figure shows the relative contributions from hi- k observed (window-convolved) modes to physical modes closeto k max “ . h Mpc ´ (vertical black line). The left panel shows a histogram of the contributions to a single physical mode.Since there are modes above k max that are contributing in the deconvolution of the window (the inversion of Eq. (66)) weinclude all the modes up to the dashed lines at k large , which are colored by survey. The specific physical mode chosen is theone that has the most contribution from modes above k large (the dashed line), and the dashed lines are chosen so that theircumulative contribution is less than . The right panel shows the cumulative contribution from observed modes k ą k large for each physical mode k NL . The vertical lines are the same as in the left plot, the horizontal line marks our target maximumcontribution of . .
50 0 .
75 1 .
00 1 .
25 1 .
50 1 .
75 2 . z r in Mpc /h N u m b e r d e n s i t y n ( z ) i n − h / M p c Input n ( r )Log-normal n ( r ) Log-Normal Angular Distribution
100 200
Number of galaxies
FIG. 10. The plots show the approximate Roman radial selection function and angular mask. In addition one realization of alog-normal simulation is shown.
A. Roman
In this section we apply the SFB power spectrum esti-mator to a log-normal simulation for the High-LatitudeSpectroscopy Survey (HLSS) of the
Nancy Grace RomanSpace Telescope (Spergel et al. „ of the sky. Ourmain objective here, however, is to exploit the large ra-dial selection that Roman will provide.The grism spectroscopy of the HLSS will yield observed wavelengths µ m to . µ m , which for the H α line at . Å results in a redshift range . ď z ď . ,and for the simulations we round this to the range h ´ Mpc to h ´ Mpc . The radial selection forour simulation is shown in the left panel of Fig. 10 (Eifler et al. https://roman.ipac.caltech.edu/sims/Param_db.html .
01 0 .
02 0 .
03 0 .
04 0 .
05 0 .
06 0 . k in h/ Mpc0200040006000800010000 N ‘ ( k , k ) i n [ h / M p c ] n Theory N shot Measured N shot .
01 0 .
02 0 .
03 0 .
04 0 .
05 0 .
06 0 . k in h/ Mpc − C ‘ ( k , k ) i n [ h / M p c ] Input P ( k )Theory C Measured C A .
01 0 .
02 0 .
03 0 .
04 0 .
05 0 .
06 0 . k in h/ Mpc0200040006000800010000 ¯ N ‘ ( k , k ) i n [ h / M p c ] n Theory N shot Measured ¯ N shot .
01 0 .
02 0 .
03 0 .
04 0 .
05 0 .
06 0 . k in h/ Mpc0250050007500100001250015000 ¯ C ‘ ( k , k ) i n [ h / M p c ] Input P ( k )Theory C Measured ¯ C A FIG. 11. SFB power spectrum measurement from uniform shot-noise-only (left column) and log-normal (right column) simula-tions with Roman window function. The top row is for a single simulation, the bottom shows the average over 50 simulations.In each panel, the grey-painted modes near k max are incomplete bandpower bins, and we show them for illustration only. Onthe right panels, the theory points take into account the bandpower binning with ∆ (cid:96) “ . The horizontal line in each plotis n . The local average effect suppresses the largest-scale mode in the simulations. Note that here we only estimate the H α sample of Roman . tion of our log-normal simulation (Agrawal et al. b “ . . The number of galaxies is „ . ˆ . We use aflat Λ CDM Planck cosmology.The results of the SFB analysis are shown in Fig. 11.The top left panel shows the shot noise from a simulationwith vanishing power spectrum as well as the theoreticalshot noise prediction from Eq. (80). The bottom leftshows the average over 50 simulations.For a simulation with signal, we show the SFB mea-surement from our log-normal simulation in the top rightof Fig. 11. Here, we have subtracted the theoretical shot noise. The theory points and the input P p k q differ dueto application of Eq. (88) to account for the bandpowerbinning with ∆ (cid:96) “ . The bottom right panel showsthe same as an average over our ensemble of simulations.We get good agreement between the ensemble measure-ment and theory power spectra if we restrict ourselves tothe blue modes in Fig. 11. The greyed-out modes cannotbe fully deconvolved due to the possibility of high- k con-tributions, as explained at the beginning of Section IV.Clearly, our estimations there were conservative, becausehigh- k contributions vary in sign and can cancel eachother.The black curves in the lower left panel of Fig. 11 con-nect theory points with the same (cid:96) . For a given (cid:96) , then,the limber ratio Eq. (36) suggests that higher k corre-sponds to lower redshift. Since the number density tendsto be higher at lower redshift, we expect the shot noise0to decrease with k given constant (cid:96) . This effect is muchmore pronounced for SPHEREx below.
B. SPHEREx
In this section we aim to show the feasibility of apply-ing our estimator to
SPHEREx (Doré et al. SPHEREx publicproducts , and for demonstration we limit ourselves tothe range ď r ď h ´ Mpc , corresponding to amaximum redshift of 0.83. We impose this limit due toour current lack of a large number of better simulations.The radial selection is shown in Fig. 12.
SPHEREx isable to go down to essentially z “ since it is not limitedby the detection of a single emission line but measuresgalaxy redshifts with 102 narrow photometric bands.For the mask we use the HFI GAL080 mask with noapodization from the Planck
Collaboration . This cutsout the galactic plane, as shown in Fig. 12. The binningstrategy in Eqs. (120) and (121) yields no binning, or ∆ (cid:96) “ ∆ n “ .Since our log-normal simulations do not take into ac-count redshift-evolution effects, we choose a fixed effec-tive redshift z eff “ . and galaxy bias b “ . .The estimation of the SFB shot noise and power spec-trum is shown in Fig. 13. The “dotted curves” of thetheory shot noise that rise quickly and then settle on anapproximately constant value are curves of constant n ,and each dot along a curve signifies the increase of (cid:96) byone. That is, lines of constant (cid:96) start on the first of thesecurves, and then decrease rapidly, as expected from theLimber ratio Eq. (36) in conjuction with a high numberdensity at low redshifts.Fig. 13 shows that we get good agreement between ourmeasured SFB power spectrum with the theory powerspectrum. It is only in the greyed-out modes that arenot fully deconvolved that a spurious oscillatory patternis introduced, and measuring those modes accurately issimply a matter of increasing k large .Furthermore, since in this paper we are primarily inter-ested in testing the SFB estimator, Fig. 13 shows everymode by itself. A more intuitive visualization of the con-straining power of the survey should bin the informationfrom neighboring modes, and this would bring the errorbars down significantly. We leave such visualization to afuture paper. For now, every mode for itself! http://spherex.caltech.edu https://github.com/SPHEREx/Public-products https://irsa.ipac.caltech.edu/data/Planck/release_2/ancillary-data/previews/HFI_Mask_GalPlane-apo0_2048_R2.00/index.html C. Euclid
As a final test for a realistic mask and selection, inthis section we apply the SFB estimator to make a fore-cast for the spectroscopic survey of
Euclid , which is anall-sky mission that covers approximately
40 % of thesky. For the number density, we adopt the reference casegiven in Amendola et al. (2018). We use the radial range ď rh ´ Mpc ď , shown in Fig. 14. Also shown inthe figure is the mask that cuts the galactic and eclipticplanes (Markovič 2020).Our simulation results are shown in Fig. 15, for bothshot noise-only and with a power spectrum signal withand galaxy bias b “ . .As for the other surveys, once we ignore the modesthat are not fully deconvolved, we get good agreementfor both cases. V. CONCLUSION
In this paper we present a new pseudo-SFB powerspectrum estimator,
SuperFaB . The estimator analyti-cally accounts for shot noise, mask, and selection effects.We also investigate the impact of the local average effectand the covariance matrix.
SuperFaB works by performing the radial transformbefore the angular transform, similar to Leistedt et al. (2012). In the radial direction the galaxies are treatedas point-particles so that no radial pixel window needscorrection. The angular transform is performed usingHEALPix (Górski et al. et al. r min and r max , as sug-gested by Samushia (2019). The boundary at r min ‰ eliminates the need for bandpower-binning in the radialdirection, the boundary at r max discretizes the measuredmodes k n(cid:96) for integer n and (cid:96) .We demonstrate that SuperFaB will be able to analyzeall large-scale modes of upcoming wide and deep galaxysurveys such as
Roman , SPHEREx , and
Euclid .We also review the SFB power spectrum theory andprovide intuition using the Limber approximation. No-tably, redshift-dependence of the power spectrum andbias factors primarily enter as the interplay between k and (cid:96) modes, such that r „ p (cid:96) ` . q{ k is an approxima-tion for the angular diameter distance. We leave a moredetailed analysis for the projection of the 3D power spec-trum to SFB space and the connection with cosmologicalparameters for a future paper.We also leave for a future paper the extension of theestimator to cross-correlations between samples with dif-fering selection functions.Since the SFB power spectrum is uniquely suited forall-sky surveys, we expect a particularly intriguing ap-plication of SuperFaB will be intensity mapping at highredshift, and we look forward to this possibility.1 . . . . . z r in Mpc /h N u m b e r d e n s i t y n ( z ) i n − h / M p c Input n ( r )Log-normal n ( r ) Log-Normal Angular Distribution
81 248
Number of galaxies
FIG. 12. The plots show the approximate radial selection function and angular mask for our
SPHEREx -like survey. .
00 0 .
01 0 .
02 0 .
03 0 .
04 0 .
05 0 . k in h/ Mpc05000100001500020000 ¯ N ‘ ( k , k ) i n [ h / M p c ] n Theory N shot Measured ¯ N shot .
00 0 .
01 0 .
02 0 .
03 0 .
04 0 .
05 0 . k in h/ Mpc − − ¯ C ‘ ( k , k ) i n [ h / M p c ] Input P ( k )Theory C Measured ¯ C A FIG. 13. The plot on the left shows the shot noise and the plot on the right the SFB power spectrum similar to Fig. 11, butnow for a full-sky mission like
SPHEREx , averaged over 50 simulations. Note that while the mask and selection are realistic,we only consider a small part of the full
SPHEREx volume due to limitations of our mocks.
For surveys covering a compact area on the sky we canintroduce boundary conditions to the angular basis func-tions as well, as pointed out by Samushia (2019). Thiswould significantly reduce the size of the computationalproblem in these cases.
ACKNOWLEDGMENTS ©2020. All rights reserved. The authors thank Kata-rina Markovič, Christopher Hirata, Lado Samushia, ChenHeinrich, and the Nancy Grace Roman Telescope
Cos-mology with the High Latitude Survey
Science Investi- gation Team (SIT) for discussions and providing selec-tion functions and masks. Part of this work was done atJet Propulsion Laboratory, California Institute of Tech-nology, under a contract with the National Aeronauticsand Space Administration. This work was supportedby NASA grant 15-WFIRST15-0008
Cosmology with theHigh Latitude Survey
Roman Science Investigation Team(SIT). Henry S. G. Gebhardt’s research was supportedby an appointment to the NASA Postdoctoral Programat the Jet Propulsion Laboratory, administered by Uni-versities Space Research Association under contract withNASA.
D. Spergel, N. Gehrels, C. Baltay, D. Bennett, J. Breckin-ridge, M. Donahue, A. Dressler, B. S. Gaudi, T. Greene, O. Guyon, et al. , arXiv e-prints , arXiv:1503.03757 (2015), .
75 1 .
00 1 .
25 1 .
50 1 .
75 2 . z r in Mpc /h N u m b e r d e n s i t y n ( z ) i n − h / M p c Input n ( r )Log-normal n ( r ) Log-Normal Angular Distribution
140 301
Number of galaxies
FIG. 14. The plots show the approximate radial selection function and angular mask for
Euclid with one log-normal simulation. .
00 0 .
01 0 .
02 0 .
03 0 .
04 0 .
05 0 . k in h/ Mpc − − ¯ N ‘ ( k , k ) i n [ h / M p c ] n Theory N shot Measured ¯ N shot .
00 0 .
01 0 .
02 0 .
03 0 .
04 0 .
05 0 . k in h/ Mpc − ¯ C ‘ ( k , k ) i n [ h / M p c ] Input P ( k )Theory C Measured ¯ C A FIG. 15. The plot on the left shows the shot noise and the plot on the right the SFB power spectrum similar to Fig. 11, butnow for a full-sky mission like
Euclid , averaged over 20 simulations.arXiv:1503.03757 [astro-ph.IM].O. Doré, J. Bock, M. Ashby, P. Capak, A. Cooray, R. dePutter, T. Eifler, N. Flagey, Y. Gong, S. Habib, et al. , arXive-prints , arXiv:1412.4872 (2014), arXiv:1412.4872 [astro-ph.CO].L. Amendola, S. Appleby, A. Avgoustidis, D. Bacon, T. Baker,M. Baldi, N. Bartolo, A. Blanchard, C. Bonvin, S. Bor-gani, et al. , Living Reviews in Relativity , 2 (2018),arXiv:1606.00180 [astro-ph.CO].DESI Collaboration, A. Aghamousa, J. Aguilar, S. Ahlen,S. Alam, L. E. Allen, C. Allende Prieto, J. Annis, S. Bailey,C. Balland, et al. , arXiv e-prints , arXiv:1611.00036 (2016),arXiv:1611.00036 [astro-ph.IM].M. Takada, R. S. Ellis, M. Chiba, J. E. Greene, H. Aihara,N. Arimoto, K. Bundy, J. Cohen, O. Doré, G. Graves, et al. ,PASJ , R1 (2014), arXiv:1206.0737 [astro-ph.CO]. D. Munshi, G. Pratten, P. Valageas, P. Coles, and P. Brax,MNRAS , 1627 (2016), arXiv:1508.00583 [astro-ph.CO].D. S. Salopek and J. R. Bond, Phys. Rev. D , 3936 (1990).E. Komatsu and D. N. Spergel, Phys. Rev. D , 063002(2001), arXiv:astro-ph/0005036 [astro-ph].N. Dalal, O. Doré, D. Huterer, and A. Shirokov, Phys. Rev. D , 123514 (2008), arXiv:0710.4560 [astro-ph].V. Desjacques, D. Jeong, and F. Schmidt, Phys. Rep. , 1(2018), arXiv:1611.09787 [astro-ph.CO].K. Yamamoto, M. Nakamichi, A. Kamino, B. A. Bassett,and H. Nishioka, PASJ , 93 (2006), arXiv:astro-ph/0505115[astro-ph].D. Bianchi, H. Gil-Marín, R. Ruggeri, and W. J. Percival,MNRAS , L11 (2015), arXiv:1505.05341 [astro-ph.CO].R. Scoccimarro, Phys. Rev. D , 083532 (2015),arXiv:1506.02729 [astro-ph.CO]. F. Beutler, H.-J. Seo, A. J. Ross, P. McDonald, S. Saito,A. S. Bolton, J. R. Brownstein, C.-H. Chuang, A. J.Cuesta, D. J. Eisenstein, et al. , MNRAS , 3409 (2017),arXiv:1607.03149 [astro-ph.CO].G. Pratten and D. Munshi, MNRAS , 3792 (2013),arXiv:1301.3673 [astro-ph.CO].P. Reimberg, F. Bernardeau, and C. Pitrou, J. CosmologyAstropart. Phys. , 048 (2016), arXiv:1506.06596 [astro-ph.CO].J. Binney and T. Quinn, MNRAS , 678 (1991).O. Lahav, in
Cosmic Velocity Fields , Vol. 9, edited byF. Bouchet and M. Lachieze-Rey (1993) p. 205, arXiv:astro-ph/9309030 [astro-ph].A. F. Heavens and A. N. Taylor, MNRAS , 483 (1995),arXiv:astro-ph/9409027 [astro-ph].H. Tadros, W. E. Ballinger, A. N. Taylor, A. F. Heav-ens, G. Efstathiou, W. Saunders, C. S. Frenk, O. Keeble,R. McMahon, S. J. Maddox, et al. , MNRAS , 527 (1999),arXiv:astro-ph/9901351 [astro-ph].W. J. Percival, D. Burkey, A. Heavens, A. Taylor, S. Cole,J. A. Peacock, C. M. Baugh, J. Bland-Hawthorn, T. Bridges,R. Cannon, et al. , MNRAS , 1201 (2004), arXiv:astro-ph/0406513 [astro-ph].B. Leistedt, A. Rassat, A. Réfrégier, and J. L. Starck, A&A , A60 (2012), arXiv:1111.3591 [astro-ph.CO].M. S. Wang, S. Avila, D. Bianchi, R. Crittenden, and W. J.Percival, J. Cosmology Astropart. Phys. , 022 (2020),arXiv:2007.14962 [astro-ph.CO].L. Samushia, arXiv e-prints , arXiv:1906.05866 (2019),arXiv:1906.05866 [astro-ph.CO].K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K.Hansen, M. Reinecke, and M. Bartelmann, ApJ , 759(2005), arXiv:astro-ph/0409513 [astro-ph].A. Zonca, L. Singer, D. Lenz, M. Reinecke, C. Rosset,E. Hivon, and K. Gorski, The Journal of Open Source Soft-ware , 1298 (2019).E. Hivon, K. M. Górski, C. B. Netterfield, B. P. Crill,S. Prunet, and F. Hansen, ApJ , 2 (2002), arXiv:astro-ph/0105302 [astro-ph].D. Alonso, J. Sanchez, A. Slosar, and LSST Dark En-ergy Science Collaboration, MNRAS , 4127 (2019),arXiv:1809.09603 [astro-ph.CO].S. Camera, J. Fonseca, R. Maartens, and M. G. Santos, MN-RAS , 1251 (2018), arXiv:1803.10773 [astro-ph.CO].A. Nicola, A. Refregier, A. Amara, and A. Paranjape,Phys. Rev. D , 063515 (2014), arXiv:1405.3660 [astro-ph.CO].F. Lanusse, A. Rassat, and J. L. Starck, A&A , A10(2015), arXiv:1406.5989 [astro-ph.CO].E. Castorina and M. White, MNRAS , 4403 (2018),arXiv:1709.09730 [astro-ph.CO].F. Beutler, E. Castorina, and P. Zhang, J. Cosmology As-tropart. Phys. , 040 (2019), arXiv:1810.05051 [astro-ph.CO].H. S. Grasshorn Gebhardt and D. Jeong, Phys. Rev. D ,083521 (2020), arXiv:2008.08706 [astro-ph.CO]. M. LoVerde and N. Afshordi, Phys. Rev. D , 123506 (2008),arXiv:0809.5112 [astro-ph].F. Bernardeau, S. Colombi, E. Gaztañaga, and R. Scocci-marro, Phys. Rep. , 1 (2002), arXiv:astro-ph/0112551[astro-ph].K. B. Fisher, O. Lahav, Y. Hoffman, D. Lynden-Bell, andS. Zaroubi, MNRAS , 885 (1995), arXiv:astro-ph/9406009[astro-ph].A. Taruya, T. Nishimichi, and D. Jeong, Phys. Rev. D ,023501 (2021), arXiv:2007.05504 [astro-ph.CO].H. A. Feldman, N. Kaiser, and J. A. Peacock, ApJ , 23(1994), arXiv:astro-ph/9304022 [astro-ph].B. D. Wandelt, E. Hivon, and K. M. Górski, Phys. Rev. D , 083003 (2001), arXiv:astro-ph/0008111 [astro-ph].A. de Mattia, V. Ruhlmann-Kleider, A. Raichoor, A. J. Ross,A. Tamone, C. Zhao, S. Alam, S. Avila, E. Burtin, J. Bautista, et al. , MNRAS , 5616 (2021), arXiv:2007.09008 [astro-ph.CO].T. Sunayama, M. Takada, M. Reinecke, R. Makiya,T. Nishimichi, E. Komatsu, S. Saito, N. Tamura, andK. Yabe, J. Cosmology Astropart. Phys. , 057 (2020),arXiv:1912.06583 [astro-ph.CO].P. J. E. Peebles, ApJ , 413 (1973).A. Dresden, Science , 393 (1920).R. Penrose, Proceedings of the Cambridge Philosophical So-ciety , 406 (1955).F. Beutler, S. Saito, H.-J. Seo, J. Brinkmann, K. S.Dawson, D. J. Eisenstein, A. Font-Ribera, S. Ho, C. K.McBride, F. Montesano, et al. , MNRAS , 1065 (2014),arXiv:1312.4611 [astro-ph.CO].A. de Mattia and V. Ruhlmann-Kleider, J. Cosmology As-tropart. Phys. , 036 (2019), arXiv:1904.08851 [astro-ph.CO].R. de Putter, C. Wagner, O. Mena, L. Verde, and W. J.Percival, J. Cosmology Astropart. Phys. , 019 (2012),arXiv:1111.6596 [astro-ph.CO].D. Wadekar, M. M. Ivanov, and R. Scoccimarro, Phys. Rev. D , 123521 (2020), arXiv:2009.00622 [astro-ph.CO].D. Wadekar and R. Scoccimarro, Phys. Rev. D , 123517(2020), arXiv:1910.02914 [astro-ph.CO].N. S. Sugiyama, S. Saito, F. Beutler, and H.-J. Seo, MNRAS , 1684 (2020), arXiv:1908.06234 [astro-ph.CO].F. Lacasa and J. Grain, A&A , A61 (2019),arXiv:1809.05437 [astro-ph.CO].Y. Li, M. Schmittfull, and U. Seljak, J. Cosmology Astropart.Phys. , 022 (2018), arXiv:1711.00018 [astro-ph.CO].B. Leistedt, H. V. Peiris, D. J. Mortlock, A. Benoit-Lévy,and A. Pontzen, MNRAS , 1857 (2013), arXiv:1306.0005[astro-ph.CO].T. Eifler, H. Miyatake, E. Krause, C. Heinrich, V. Miranda,C. Hirata, J. Xu, S. Hemmati, M. Simet, P. Capak, et al. ,arXiv e-prints , arXiv:2004.05271 (2020), arXiv:2004.05271[astro-ph.CO].A. Agrawal, R. Makiya, C.-T. Chiang, D. Jeong, S. Saito,and E. Komatsu, J. Cosmology Astropart. Phys. , 003(2017), arXiv:1706.09195 [astro-ph.CO]. K. Markovič, personal communication, paper in prep (2020).D. N. Limber, ApJ , 655 (1954).
Appendix A: Useful formulae
For any function f p k q ż k d k d ˆ k δ D p k ´ k q f p k q“ ż d k δ D p k ´ k q ż d ˆ k δ D p ˆ k ´ ˆ k q f p k q . (A1)Therefore, δ D p k ´ k q “ k ´ δ D p k ´ k q δ D p ˆ k ´ ˆ k q . (A2)Furthermore, r δ D ˆ r ´ r ˙ “ r δ D p r ´ r q . (A3)The first-order result from LoVerde and Afshordi(2008) can be written as J ν p kr q » δ D p kr ´ ν q , (A4)where J ν p x q is the Bessel function. Therefore, for a spher-ical Bessel function j (cid:96) p x q “ a π { x J (cid:96) ` p x q we get j (cid:96) p kr q » c π rk k δ D ˆ r ´ (cid:96) ` k ˙ (A5)to first order. This is called Limber’s approximation (Limber 1954).Spherical Bessel functions and spherical harmonics sat-isfy orthogonality relations δ D p k ´ k q “ kk π ż d r r j (cid:96) p kr q j (cid:96) p k r q , (A6) δ K(cid:96)(cid:96) δ Kmm “ ż dΩ ˆ r Y (cid:96)m p ˆ r q Y ˚ (cid:96) m p ˆ r q . (A7)Spherical harmonics can be expressed in terms of a com-plex exponential and real associated Legendre functions P m(cid:96) p x q as Y (cid:96)m p ˆ r q “ e imφ ˆ p (cid:96) ´ m q ! p (cid:96) ` q π p (cid:96) ` m q ! ˙ P m(cid:96) p cos θ q . (A8)The completeness relation is ÿ (cid:96)m Y (cid:96)m p ˆ r q Y ˚ (cid:96)m p ˆ r q “ δ D p ˆ r ´ ˆ r q . (A9)Rayleigh’s formula decomposes the plane waves intospherical Bessels and spherical harmonics, e i q ¨ r “ π ÿ (cid:96) ,m i (cid:96) j (cid:96) p qr q Y ˚ (cid:96) m p ˆ q q Y (cid:96) m p ˆ r q . (A10) Legendre polynomials can be expressed as a sum overspherical harmonics as P (cid:96) p ˆ k ¨ ˆ r q “ π (cid:96) ` ÿ m Y (cid:96)m p ˆ k q Y ˚ (cid:96)m p ˆ r q . (A11)Flipping the sign of the component angular momentumor the direction of the argument to spherical harmonicsgives Y ˚ (cid:96)m p ˆ r q “ p´ q m Y (cid:96), ´ m p ˆ r q , (A12) Y (cid:96)m p´ ˆ r q “ p´ q (cid:96) Y (cid:96),m p ˆ r q . (A13)The Gaunt factor is G (cid:96)LL mMM “ ż d ˆ r Y (cid:96)m p ˆ r q Y LM p ˆ r q Y L M p ˆ r q , (A14)and it can be expressed in terms of Wigner- j symbols, G (cid:96)LL mMM “ ˆ p (cid:96) ` qp L ` qp L ` q π ˙ ˆ (cid:96) L L ˙ ˆ ˆ (cid:96) L L m M M ˙ . (A15)The Wigner j symbols obey an orthogonality relation ÿ mM ˆ (cid:96) L L m M M ˙ ˆ (cid:96) L L m M M ˙ “ δ KL L δ KM M δ T p (cid:96), L, L q L ` , (A16)where δ T p (cid:96), L, L q enforces the triangle relation that isalso obeyed by the 3 j -symbols. That is, the Gaunt factoris only nonzero when m ` M ` M “ , (A17) | (cid:96) ´ L | ď L ď (cid:96) ` L . (A18)Assuming the triangle condition is satisfied, for even J “ (cid:96) ` L ` L we have ˆ (cid:96) L L ˙ “ p´ q J ˆ p J ´ (cid:96) q ! p J ´ L q ! p J ´ L q ! p J ` q ! ˙ ˆ ` J ˘ ! ` J ´ (cid:96) ˘ ! ` J ´ L ˘ ! ` J ´ L ˘ ! , (A19)for odd J “ (cid:96) ` L ` L , those j ’s vanish when m “ M “ M “ . Appendix B: The Laplacian in an expandinguniverse
We use the SFB decomposition because that correlatesradial and angular modes of the same scale. Here we showthat the Laplacian in a flat expanding universe takes theform Eq. (38).5The flat Robertson-Walker metric is ds “ ´ dt ` a “ dr ` r dθ ` r sin θ dφ ‰ , (B1)where r is the comoving coordinate, and the metric hasthe non-zero Christoffel symbols Γ “ a H , Γ “ a r H , Γ “ a r sin θ H , Γ “ Γ “ H , Γ “ ´ r , Γ “ ´ r sin θ , Γ “ Γ “ H , Γ “ Γ “ r , Γ “ ´ cos θ sin θ , Γ “ Γ “ H , Γ “ Γ “ r , Γ “ Γ “ cos θ sin θ . Therefore, g µν Γ µν “ H , (B2) g µν Γ µν “ ´ a ´ r ´ , (B3) g µν Γ µν “ ´ a ´ r ´ cos θ sin θ , (B4) g µν Γ µν “ . (B5)The d’Alembertian operator for a scalar f is l f “ g µν f : µν “ g µν ` B µ B ν ´ Γ σµν B σ ˘ f “ „ ´ B ` a ´ B ` a ´ r ´ B ` a ´ r ´ sin ´ θ B ´ H B ` a ´ r ´ B ` a ´ r ´ cos θ sin θ B f “ „ ´ B ´ H B ` a ´ ˆ B ` r ´ B ` r ´ sin ´ θ B ` r ´ B ` r ´ cos θ sin θ B ˙ f . (B6)Identifying the term in parentheses as given by Eq. (38),we get l f “ “ ´B ´ H B ` a ´ ∇ ‰ f (B7)in a flat expanding universe. We can also write Eq. (B7)as l f “ “ ´ a ´ B t ` a B t ˘ ` a ´ ∇ ‰ f . (B8)The eigenfunctions to the d’Alembertian are, therefore,separable. If we write an eigenfunction f p t, r, ˆ r q “ p p t q g p r q h p ˆ r q (B9)with ∇ r g p r q h p ˆ r qs “ ´ k g p r q h p ˆ r q , (B10)then “ a ´ B t ` a B t p ˘ ` “ a ´ k ´ λ ‰ p , (B11)where ´ λ is the eigenvalue of the d’Alembertian. Since ar is the angular diameter distance and r is the comovingdistance (also comoving angular diameter distance), wecan call k a comoving mode. Appendix C: Radial spherical Fourier-Bessel modeswith potential boundary conditions
In this appendix we derive the radial basis functionsof the Laplacian with potential boundary conditions at r min and r max .We first isolate the radial part of Eq. (38). Writing f p r q “ g p r q h p ˆ r q , (C1)we require ´ k gh “ hr BB r ˆ r B g B r ˙ ` g ∇ h . (C2)Given the spherical harmonic solution for the angularterm h , ∇ h “ ´ (cid:96) p (cid:96) ` q r h , (C3)we get “ dd r ˆ r d g (cid:96) p kr q d r ˙ ` “ p kr q ´ (cid:96) p (cid:96) ` q ‰ g (cid:96) p kr q , (C4)where we now added that the function g depends on (cid:96) .Our first aim is to derive the discrete spectrum of k modesfor a given (cid:96) . We then use that to derive the form of the g (cid:96) . Following Fisher et al. (1995), we demand that theorthogonality relation Eq. (42) is satisfied. However, wemodify the approach in Fisher et al. (1995) to integratefrom r min to r max . Eq. (C4) multiplied by g (cid:96) p kr q thenyields ż r max r min d r dd r ˆ r d g (cid:96) p kr q d r ˙ g (cid:96) p k r q“ ż r max r min d r “ (cid:96) p (cid:96) ` q ´ p kr q ‰ g (cid:96) p kr q g (cid:96) p k r q . (C5)Subtract from this equation the same equation with k and k interchanged, “ k ´ k ‰ ż r max r min d r r g (cid:96) p kr q g (cid:96) p k r q“ ż r max r min d r " dd r ˆ r d g (cid:96) p kr q d r ˙ g (cid:96) p k r q´ dd r ˆ r d g (cid:96) p k r q d r ˙ g (cid:96) p kr q * . (C6)Partial integration with the terms on the right hand side(r.h.s.) yields ż d r dd r ` kr g (cid:96) p kr q ˘ g (cid:96) p k r q“ kr g (cid:96) p kr q g (cid:96) p k r q ˇˇ r max r min ´ kk ż d r r g (cid:96) p kr q g (cid:96) p k r q . (C7)6Then, Eq. (C6) becomes “ k ´ k ‰ ż r max r min d r r g (cid:96) p kr q g (cid:96) p k r q“ kr g (cid:96) p kr q g (cid:96) p k r q ˇˇ r max r min ´ k r g (cid:96) p k r q g (cid:96) p kr q ˇˇ r max r min . (C8)The r.h.s. will vanish for any k whenever “ Akr g (cid:96) p kr max q ´ Br g (cid:96) p kr max q´ akr g (cid:96) p kr min q ` br g (cid:96) p kr min q . (C9) for any constants a , b , A , and B .To choose a , b , A , and B , we note that the repre-sentable field δ p r q is written as a sum of the solutions toEq. (C4) inside the SFB volume, and we have some free-dom to choose the desired behavior outside of it. Sincethe field inside the SFB volume satisfies the Poisson equa-tion, it is natural to have it satisfy Laplace’s equation outside it, and demand that the solution is continuousand smooth at the boundaries. That is, δ p r q “ $’’’’’&’’’’’%ř (cid:96)m „ a (cid:96)m ´ rr min ¯ (cid:96) ` b (cid:96)m ` r min r ˘ (cid:96) ` Y (cid:96)m p ˆ r q , for r ă r min , ř n(cid:96)m ” c n(cid:96) j (cid:96) p k n(cid:96) r q ` d n(cid:96) y (cid:96) p k n(cid:96) r q ı Y (cid:96)m p ˆ r q δ n(cid:96)m , for r min ď r ď r max , ř (cid:96)m „ A (cid:96)m ´ rr max ¯ (cid:96) ` B (cid:96)m ` r max r ˘ (cid:96) ` Y (cid:96)m p ˆ r q , for r ą r max , (C10)where we defined the constants a (cid:96)m , b (cid:96)m , c n(cid:96) , d n(cid:96) , A (cid:96)m ,and B (cid:96)m , and we explicitly wrote g (cid:96) p kr q “ c n(cid:96) j (cid:96) p kr q ` d n(cid:96) y (cid:96) p kr q , (C11)and we anticipate that the function g (cid:96) will also dependon n . Continuity at the boundaries requires a (cid:96)m ` b (cid:96)m “ ÿ n g (cid:96) p k n(cid:96) r min q δ n(cid:96)m , (C12) A (cid:96)m ` B (cid:96)m “ ÿ n g (cid:96) p k n(cid:96) r max q δ n(cid:96)m . (C13)Smoothness further requires (cid:96) a (cid:96)m r min ´ p (cid:96) ` q b (cid:96)m r min “ ÿ n k n(cid:96) g (cid:96) p k n(cid:96) r min q δ n(cid:96)m , (C14) (cid:96) A (cid:96)m r max ´ p (cid:96) ` q B (cid:96)m r max “ ÿ n k n(cid:96) g (cid:96) p k n(cid:96) r max q δ n(cid:96)m , (C15)Now requiring δ p r q to be finite at r “ and r “ 8 sets b (cid:96)m “ A (cid:96)m “ , and requiring continuity and smoothnessfor any δ n(cid:96)m , we get (cid:96) g (cid:96) p k n(cid:96) r min q “ k n(cid:96) r min g (cid:96) p k n(cid:96) r min q , (C16) ´p (cid:96) ` q g (cid:96) p k n(cid:96) r max q “ k n(cid:96) r max g (cid:96) p k n(cid:96) r max q . (C17)These choices lead to a “ , b “ (cid:96) { r min , A “ , and B “ ´p (cid:96) ` q{ r max in Eq. (C9), which shows that the Laplace’s equation is Poisson’s equation without a source term. conditions Eqs. (C16) and (C17) on k n(cid:96) lead to an or-thogonality relation for the g (cid:96) .Both j (cid:96) and y (cid:96) satisfy the two recurrence relations j (cid:96) p kr q “ ´ j (cid:96) ` p kr q ` (cid:96)kr j (cid:96) p kr q , (C18) j (cid:96) p kr q “ j (cid:96) ´ p kr q ´ (cid:96) ` kr j (cid:96) p kr q . (C19)Then, Eqs. (C16) and (C17) simplify to c n(cid:96) j (cid:96) ` p k n(cid:96) r min q ` d n(cid:96) y (cid:96) ` p k n(cid:96) r min q “ , (C20) c n(cid:96) j (cid:96) ´ p k n(cid:96) r max q ` d n(cid:96) y (cid:96) ´ p k n(cid:96) r max q “ . (C21)The normalization of g (cid:96) is obtained by dividing Eq. (C8)by k ´ k , and taking the limit k Ñ k “ k n(cid:96) , “ ż r max r min d r r g (cid:96) p kr q“ lim k Ñ k kr g (cid:96) p kr q g (cid:96) p k r q ˇˇ r max r min ´ k r g (cid:96) p k r q g (cid:96) p kr q ˇˇ r max r min k ´ k . (C22)When kr “ k n(cid:96) r min as in Eq. (C16), kr g (cid:96) p kr q g (cid:96) p k r q ´ k r g (cid:96) p k r q g (cid:96) p kr q“ k r g (cid:96) p kr q “ c n(cid:96) j (cid:96) ` p k r q ` d n(cid:96) y (cid:96) ` p k r q ‰ , (C23)and when kr “ k n(cid:96) r max as in Eq. (C17), kr g (cid:96) p kr q g (cid:96) p k r q ´ k r g (cid:96) p k r q g (cid:96) p kr q“ ´ k r g (cid:96) p kr q “ c n(cid:96) j (cid:96) ´ p k r q ` d n(cid:96) y (cid:96) ´ p k r q ‰ . (C24)7The terms in brackets vanish in the limit k Ñ k “ k n(cid:96) as per Eqs. (C20) and (C21). That is, we need limits lim q Ñ q c n(cid:96) j (cid:96) ` p q q ` d n(cid:96) y (cid:96) ` p q q q ´ q “ g (cid:96) p q q q , (C25) lim q Ñ q c n(cid:96) j (cid:96) ´ p q q ` d n(cid:96) y (cid:96) ´ p q q q ´ q “ ´ g (cid:96) p q q q , (C26)for q “ kr min and q “ kr max , respectively. Then,Eq. (C22) becomes “ r g (cid:96) p k n(cid:96) r max q ´ r g (cid:96) p k n(cid:96) r min q . (C27)Choosing k n(cid:96) , c n(cid:96) , and d n(cid:96) that satisfy Eqs. (C20),(C21) and (C27) guarantees the orthonormality of the g (cid:96) , ż r max r min d r r g (cid:96) p k n(cid:96) r q g (cid:96) p k n (cid:96) r q “ δ Knn . (C28)Note that the condition (cid:96) “ (cid:96) is not enforced by the g (cid:96) .Instead, (cid:96) “ (cid:96) comes from the spherical harmonics, i.e.,Eq. (A7).
1. Phase factor
We are free to introduce a phase factor for the g n(cid:96) p r q ,and we choose it so that the sign of g n(cid:96) p r min q alternateswith n , but stays constant with (cid:96) , g n(cid:96) p r q “ p´ qr ´ floor p (cid:96) ` qsr ´ floor p n qs ˜ g n(cid:96) p r q , (C29)where the tilde indicates that we have not included thephase factor. This flips the sign unless either (cid:96) ‰ or n ‰ . Thus, the basis functions in Fig. 4 are obtained.
2. Numerical concerns
To calculate k n(cid:96) , solve each of Eqs. (C20) and (C21)for the ratio d n(cid:96) { c n(cid:96) , and set them equal to get “ j (cid:96) ´ p k n(cid:96) r max q y (cid:96) ` p k n(cid:96) r min q´ y (cid:96) ´ p k n(cid:96) r max q j (cid:96) ` p k n(cid:96) r min q . (C30)Examples for the resulting zeros and the first few basisfunctions are shown in Fig. 4. The ratio d n(cid:96) { c n(cid:96) thenfollows from Eq. (C20) or Eq. (C21). Finally, the overallnormalization is fixed by Eq. (C27) up to a sign.When r max is large and r min is small, the k n(cid:96) may needto be computed using arbitrary precision floats, and the g n(cid:96) may need to be calculated with arbitrary precision aswell. Caching the result in double precision should thenprovide for sufficient speed for the actual transform. Appendix D: Covariance matrix simplification
In this appendix we simplify the covariance matrix Eq. (112). For simplicity we ignore the local average effect. Weexplicitly treat the shot noise Eq. (79) because it is inhomogeneous and anisotropic. We get A δ W,ANLM δ W,A, ˚ N L M E “ ÿ n (cid:96) m W n (cid:96) m NLM ÿ n (cid:96) m W n (cid:96) m , ˚ N L M „ δ K(cid:96) (cid:96) δ Km m C (cid:96) n n ` n ` W ´ ˘ n (cid:96) m n (cid:96) m (D1) “ ÿ (cid:96) n n C (cid:96) n n ÿ m W n (cid:96) m NLM W N L M n (cid:96) m ` n W N L M NLM , (D2)where we used Eq. (53). Similarly, when neither density contrast has the complex conjugate attached, A δ W,ANLM δ W,AN L M E “ p´ q M ÿ (cid:96) n n C (cid:96) n n ÿ m W n (cid:96) m NLM W N L , ´ M n (cid:96) m ` p´ q M ¯ n W N L , ´ M NLM , (D3)where we used Eq. (A12). Therefore, both terms in Eq. (112) are of the form p (cid:96) ` qp L ` q ÿ mM A δ W,An(cid:96)m δ W,A, ˚ NLM
E A δ W,An (cid:96)m δ W,A, ˚ N LM E “ p (cid:96) ` qp L ` q ÿ mM « ÿ (cid:96) n n C (cid:96) n n ÿ m W n (cid:96) m n(cid:96)m W NLMn (cid:96) m ` n W NLMn(cid:96)m ff ˆ « ÿ (cid:96) n n C (cid:96) n n ÿ m W n (cid:96) m , ˚ n (cid:96)m W N LM, ˚ n (cid:96) m ` n W N LM, ˚ n (cid:96)m ff (D4)8 “ p (cid:96) ` qp L ` q ÿ mM « ÿ (cid:96) n n C (cid:96) n n ÿ m W n (cid:96) m n(cid:96)m W NLMn (cid:96) m ÿ (cid:96) n n C (cid:96) n n ÿ m W n (cid:96) m , ˚ n (cid:96)m W N LM, ˚ n (cid:96) m ` n W NLMn(cid:96)m n W N LM, ˚ n (cid:96)m ` n W N LM, ˚ n (cid:96)m ÿ (cid:96) n n C (cid:96) n n ÿ m W n (cid:96) m n(cid:96)m W NLMn (cid:96) m ` @ n Ø n , N Ø N D ˚ ff (D5) “ A LNN (cid:96)nn ` A LNN (cid:96)nn ` A LN N(cid:96)n n ` A LNN (cid:96)nn . (D6)Using Eq. (53), we get A LNN (cid:96)nn “ p (cid:96) ` qp L ` q ÿ (cid:96) n n C (cid:96) n n ÿ (cid:96) n n C (cid:96) n n W ¨˝ (cid:96) L (cid:96) (cid:96)n N n n n N n n ˛‚ , (D7)and A LNN (cid:96)nn “ n p (cid:96) ` qp L ` q ÿ (cid:96) n n C (cid:96) n n W ¨˝ (cid:96) L (cid:96)n N n n N n ˛‚ , (D8)and A LNN (cid:96)nn “ n p (cid:96) ` qp L ` q W ¨˝ L (cid:96)N n N n ˛‚ “ n (cid:96) ` M LNN (cid:96)nn . (D9)The W k symbols are defined in Eq. (E1) and discussed in Appendix E. Then Eq. (112) becomes V W LNN (cid:96)nn “ A LNN (cid:96)nn ` A LNN (cid:96)nn ` A LN N(cid:96)n n ` A LNN (cid:96)nn ` @ N Ø N D . (D10)The A term dominates if the power spectrum is much larger than the shot noise, and A dominates if shot noise islarger. Appendix E: Chains of window functions
Throughout the paper, we find that traces of density-contrast window coupling matrices W appear with summationsover the azimuthal modes m . That is, we find tensors W k of the form W k ¨˝ (cid:96) (cid:96) . . . (cid:96) k n n . . . n k n n . . . n k ˛‚ “ ÿ m m ...m k W n (cid:96) m n k (cid:96) k m k W n (cid:96) m n (cid:96) m ¨ ¨ ¨ W n k (cid:96) k m k n k ´ (cid:96) k ´ m k ´ . (E1)This starts with a single window function Eq. (80) ( k “ ) for shot noise, two window functions Eqs. (67) and (D9)( k “ ) for the pseudo power spectrum mixing matrix and shot noise covariance, three window functions Eqs. (107)and (D8) for part of the local average effect and covariance calculations, four window functions Eq. (D7). If we wereto include the local average effect in the covariance, then k “ and k “ would also occur.The W k are real, which is trivially shown by substituting each window with its definition Eq. (51) and usingEq. (A11). Also, cyclical permutations of the argument columns leave W k invariant, and anticyclical permutationsleave it invariant if all n i and n i are switched as well. That is, W k ¨˝ (cid:96) (cid:96) . . . (cid:96) k n n . . . n k n n . . . n k ˛‚ “ W k ¨˝ (cid:96) k (cid:96) . . . (cid:96) k ´ n k n . . . n k ´ n k n . . . n k ´ ˛‚ “ W k ¨˝ (cid:96) k (cid:96) k ´ . . . (cid:96) n k n k ´ . . . n n k n k ´ . . . n ˛‚ . (E2)The first equality trivially follows from the definition Eq. (E1), and the second from applying Eq. (53) to all thewindows.9
1. Evaluation for separable window function
For a separable window W p r q “ φ p r q W p ˆ r q , Eq. (56) becomes W n(cid:96)mn (cid:96) m “ I n(cid:96)n (cid:96) W (cid:96)m(cid:96) m , (E3)where I n(cid:96)n (cid:96) “ ż r max r min d r r g n (cid:96) p r q g n(cid:96) p r q φ p r q , (E4) W (cid:96)m(cid:96) m “ p´ q m ÿ LM G (cid:96) (cid:96)L ´ m ,m,M W LM , (E5)where W LM was defined in Eq. (71). Thus, Eq. (E1) is W k ¨˝ (cid:96) (cid:96) . . . (cid:96) k n n . . . n k n n . . . n k ˛‚ “ I n (cid:96) n k (cid:96) k I n (cid:96) n (cid:96) ¨ ¨ ¨ I n k (cid:96) k n k ´ (cid:96) k ´ ÿ m m ...m k W (cid:96) m (cid:96) k m k W (cid:96) m (cid:96) m ¨ ¨ ¨ W (cid:96) k m k (cid:96) k ´ m k ´ . (E6)Eq. (E6) contains O ` (cid:96) k ˘˘