Beyond Yamamoto: Anisotropic Power Spectra and Correlation Functions with Pairwise Lines-of-Sight
MMNRAS , 1–27 (2021) Preprint 18 February 2021 Compiled using MNRAS L A TEX style file v3.0
Beyond Yamamoto: Anisotropic Power Spectra and CorrelationFunctions with Pairwise Lines-of-Sight
Oliver H. E. Philcox , ★ and Zachary Slepian , † Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA School of Natural Sciences, Institute for Advanced Study, 1 Einstein Drive, Princeton, NJ 08540, USA Department of Astronomy, University of Florida, 211 Bryant Space Science Center, Gainesville, FL 32611, USA Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA
18 February 2021
ABSTRACT
Conventional estimators of the anisotropic power spectrum and two-point correlation function(2PCF) adopt the ‘Yamamoto approximation’, fixing the line-of-sight of a pair of galaxies to thatof just one of its members. Whilst this is accurate only to first-order in the characteristic openingangle 𝜃 max , it allows for efficient implementation via Fast Fourier Transforms (FFTs). This workpresents practical algorithms for computing the power spectrum and 2PCF multipoles usingpairwise lines-of-sight, adopting either the galaxy midpoint or angle bisector definitions. Usingnewly derived infinite series expansions for spherical harmonics and Legendre polynomials,we construct estimators accurate to arbitrary order in 𝜃 max , though note that the midpoint andbisector formalisms themselves differ at fourth order. Each estimator can be straightforwardlyimplemented using FFTs, requiring only modest additional computational cost relative tothe Yamamoto approximation. We demonstrate the algorithms by applying them to a set ofrealistic mock galaxy catalogs, and find both procedures produce comparable results for the2PCF, with a slight preference for the bisector power spectrum algorithm, albeit at the costof greater memory usage. Such estimators provide a useful method to reduce wide-anglesystematics for future surveys. Key words: cosmology: large-scale structure of Universe, theory – methods: statistical, dataanalysis
We have now entered the epoch of ‘precision cosmology’. In the coming years, the volume of cosmological data available will increase at aprodigious rate, thanks to the advent of large spectroscopic surveys such as DESI (DESI Collaboration et al. 2016), Euclid (Laureijs et al.2011) and SPHEREx (Doré et al. 2014). As the number of observed galaxies grows, so too does the precision on fundamental parameters suchas the growth rate, energy densities and Hubble parameter. Given that we will soon be able to measure summary statistics at the sub-percentlevel, it is vital to understand also their systematics to this precision, else we risk biasing our inference or losing effective survey volume.Whilst there is growing interest in more complicated statistics (e.g., Gil-Marín et al. 2017; Slepian et al. 2017; Chudaykin & Ivanov2019; Philcox et al. 2020; Samushia et al. 2021), the information content of future surveys will be dominated by the two-point correlator,masquerading either as the configuration-space two-point correlation function (2PCF), 𝜉 ( r ) , or the Fourier-space power spectrum, 𝑃 ( k ) (e.g.,Beutler et al. 2017; Alam et al. 2017; eBOSS Collaboration et al. 2020). In a statistically isotropic universe, both of these will depend only thedistance between galaxies, be it 𝑟 = | r | or the momentum-space equivalent 𝑘 = | k | . In our Universe this is not the case, since redshift-spacedistortions (RSD) impart a preferred origin to the observer (Kaiser 1987), sourcing additional cosmological information (e.g., Lesgourgues& Pastor 2006; Weinberg et al. 2013).To fully encapsulate RSD, two-point statistics should depend on the position vectors to the two galaxies in question, rather than justa single length. Taking into account the various rotational symmetries, such a configuration can be specified by three degrees of freedom;options include the separation 𝑟 (or 𝑘 ) and the angles between the separation vector and the line-of-sight (LoS) to each galaxy (e.g., Pápai & ★ E-mail: [email protected] † E-mail: zslepian@ufl.edu © a r X i v : . [ a s t r o - ph . C O ] F e b Philcox & Slepian
Szapudi 2008; Yoo & Seljak 2015; Castorina & White 2018) or the separation, the mean distance to the galaxy pair, and a single angle (e.g.,Reimberg et al. 2016; Beutler et al. 2019). Unless our interest lies in the largest-possible scales (for instance in 𝑓 NL -analyses), it is usuallysufficient to parametrize the two-point correlators by just two variables; the inter-galaxy distance 𝑟 or 𝑘 , and the angle of the galaxy separationvector to a joint LoS, 𝜇 . In this case, the functions can be robustly expanded as a Legendre series in 𝜇 , and theory and observations simplycompared. Of course, any such approximation necessarily induces wide-angle effects on the largest scales, which are the subject of extensivediscussion in the literature (e.g., Hamilton 1992; Hamilton & Culhane 1996; Hamilton 1998; Zaroubi & Hoffman 1996; Szalay et al. 1998;Szapudi 2004; Datta et al. 2007; Pápai & Szapudi 2008; Shaw & Lewis 2008; Bonvin & Durrer 2011; Raccanelli et al. 2014; Yoo & Seljak2015; Slepian & Eisenstein 2015; Reimberg et al. 2016; Castorina & White 2018; Beutler et al. 2019). In general, the error in these approachesdepends on the characteristic opening angle 𝜃 max , defined either as the galaxy pair opening angle in the maximum radial bin (2PCF) or thesurvey opening angle (power spectrum). The dichotomy arises since the power spectrum depends on an integral over all galaxy pairs, whilstthe 2PCF only requires pairs separated by the scale of interest.If the two-parameter formalism is adopted (as has become commonplace), an important question must be asked: how should one choosethe joint LoS to the galaxy pair? The early literature adopted a single LoS for the whole survey (e.g., Kaiser 1987; Hamilton 1992), which,whilst simple to implement, incurs significant errors if the survey is wide. A more accurate prescription is to fix the LoS to the direction vectorof a single galaxy, in the ‘Yamamoto approximation’ (Yamamoto et al. 2006). Whilst this gives an error at O( 𝜃 ) for characteristic size 𝜃 max , which becomes important for DESI volumes (Sugiyama et al. 2019), it is straightforward to implement using Fast Fourier Transforms(FFTs), thus is the approach found in most recent analyses (e.g., Scoccimarro 2015; Bianchi et al. 2015; Hand et al. 2017; Beutler et al. 2017;Hand et al. 2018). Beyond this approximation, there are two appealing LoS choices: the angle between the galaxy midpoint and separationvector (cf. Yamamoto et al. 2006; Scoccimarro 2015; Samushia et al. 2015; Bianchi et al. 2015) and the angle bisector (cf. Szalay et al. 1998;Matsubara 2000; Szapudi 2004; Yoo & Seljak 2015). Both are consistent to third order in the opening angle (demonstrated in Slepian &Eisenstein 2015), but incur an error at fourth order, and, for 𝜃 max < ° , lose little information compared to the double LoS approach (Beutleret al. 2012; Samushia et al. 2012; Yoo & Seljak 2015). However, their naïve implementation scales as O( 𝑁 ) for 𝑁 galaxies, rather than the O( 𝑁 g log 𝑁 g ) dependence enjoyed by algorithms based on FFTs with 𝑁 g grid cells.In this work, we will demonstrate that the midpoint and bisector power spectrum and 2PCF estimators may be efficiently computed in O( 𝑁 g log 𝑁 g ) time using FFTs. In both cases, it is necessary to perform a series expansion of the angular dependence in the galaxy pairopening angle 𝜃 (with 𝜃 < 𝜃 max ); however, we give the explicit form of these corrections at arbitrary order. The different LoS definitionsrequire different mathematical treatments for greatest efficiency; for the bisector case, we implement the suggestion of Castorina & White(2018) and extend it to arbitrary order, whilst for the midpoint approach, we provide novel formulae based on a newly-derived sphericalharmonic shift theorem (similar to that of Garcia & Slepian 2020 for the three-point function). This paper is an extension also of Slepian &Eisenstein (2015), which gave the lowest-order corrections for the 2PCF, but did not consider the Fourier-space counterpart (which carriessomewhat more subtleties). In contrast, Samushia et al. (2015) considered the power-spectrum in the midpoint formalism, but only appliedto a simplified spherical cap geometry. Our work goes beyond the previous by giving a full catalog of arbitrary-order expressions for themidpoint and bisector formalism in real- and Fourier-space. We further consider their application to realistic data using the MultiDark-patchymock catalogs (Kitaura et al. 2016), and make the analysis code publicly available. The remainder of this work is structured as follows. We begin in §2 by recapitulating the basic two-point correlator estimators. §3 & §4present our implementations of the power spectrum and 2PCF algorithms in the midpoint formalism, before the same is done in the bisectorformalism in §5 & §6. §7 considers the application of the algorithms to data, before we conclude in §8. A list of useful mathematical identitiesis given in Appendix A, with Appendices B-F giving mathematical derivations of results central to this work, in particular, a shift theoremfor spherical harmonics and Legendre polynomials. For the reader less interested in mathematical derivations, we recommend skipping §3.3and the Appendices, and note that the key equations in this work are boxed.
We begin by stating the Fourier conventions used throughout this work. We define the Fourier and inverse Fourier transforms by 𝑋 ( k ) ≡ F [ 𝑋 ] ( k ) = ∫ 𝑑 x 𝑒 − 𝑖 k · x 𝑋 ( x ) , 𝑋 ( x ) ≡ F − [ 𝑋 ] ( x ) = ∫ 𝑑 k ( 𝜋 ) 𝑒 𝑖 k · x 𝑋 ( k ) , (2.1)leading to the definition of the Dirac delta function, 𝛿 D , as ( 𝜋 ) 𝛿 D ( k − k ) = ∫ 𝑑 x 𝑒 𝑖 ( k − k ) · x . (2.2)The correlation function and power spectrum of the density field, 𝛿 , are defined as 𝜉 ( r ) = (cid:104) 𝛿 ( x ) 𝛿 ( x + r )(cid:105) , ( 𝜋 ) 𝛿 D ( k + k (cid:48) ) 𝑃 ( k ) = (cid:10) 𝛿 ( k ) 𝛿 ( k (cid:48) ) (cid:11) , (2.3) github.com/oliverphilcox/BeyondYamamoto MNRAS , 1–27 (2021) eyond Yamamoto Correlators Fixed Yamamoto Midpoint Bisector
Figure 1.
Comparison of the various choices of lines-of-sight (LoS) for the galaxy pair. For each example, the observer is located at the based of the triangle,with the galaxies at r and r with separation vector 𝚫 = r − r . We define 𝜙 as the internal triangle angle with cos 𝜙 = ˆ r · ˆ r . In the small-angle limit, wehave 𝜙 ≈ 𝜃 , where 𝜃 ≡ Δ /| r + r | . The fixed LoS approximation (ˆ n = const . ) incurs an O ( 𝜃 ) error, whilst the Yamamoto approximation (ˆ n = ˆ r ) andmidpoint or bisector estimators (ˆ n ∝ r + r or ˆ r + ˆ r ) incur O ( 𝜃 ) and O ( 𝜃 ) errors respectively, averaging over 𝜃 < 𝜃 max . with the power spectrum as the Fourier transform of the correlation function. Additionally, we use the shorthand ∫ k ≡ ∫ 𝑑 k ( 𝜋 ) , ∫ Ω 𝑘 ≡ ∫ 𝑑 Ω 𝑘 𝜋 . (2.4) The conventional estimator for the galaxy power spectrum multipoles, 𝑃 ℓ ( 𝑘 ) , is defined as a Fourier transform of two density fields 𝛿 ( r ) and 𝛿 ( r ) : ˆ 𝑃 ℓ ( 𝑘 ) = ℓ + 𝑉 ∫ Ω 𝑘 ∫ 𝑑 r 𝑑 r 𝑒 − 𝑖 k ·( r − r ) 𝛿 ( r ) 𝛿 ( r ) 𝐿 ℓ ( ˆ k · ˆ n ) (2.5)(e.g., Hand et al. 2017), where ˆ n is the joint line-of-sight (LoS) to the pair of galaxies at ( r , r ) , 𝐿 ℓ is the Legendre polynomial of order ℓ , 𝑉 is the survey volume and hats denote unit vectors. Whilst we have not included 𝑘 -space binning, this is a straightforward addition, requiringan additional integral over | k | . In spectroscopic surveys, we do not have access to 𝛿 directly, only a set of galaxy and random particle positionswith associated density fields 𝑛 𝑔 ( r ) and 𝑛 𝑟 ( r ) respectively. In this context, we replace 𝛿 ( r ) → 𝑤 ( r )[ 𝑛 𝑔 ( r ) − 𝛼 𝑟 𝑛 𝑟 ( r )] 𝐼 / , 𝐼 ≡ ∫ 𝑑 r 𝑤 ( r ) ¯ 𝑛 ( r ) (2.6)(Yamamoto et al. 2006), where 𝑤 ( r ) are some weights (accounting for systematics, optimality and completeness), ¯ 𝑛 is the mean galaxy densityand 𝛼 𝑟 is the ratio of randoms to galaxies. (2.5) is strictly an estimator for the window-convolved power spectrum, and includes a shot-noiseterm. Since neither effect depends on the LoS choice, we will ignore both henceforth.In the simplest (‘plane-parallel’) approximation, we take ˆ n to be fixed across the survey, thus the estimator simplifies to the form familiarfrom 𝑁 -body estimators, involving the Fourier-space density field 𝛿 ( k ) : ˆ 𝑃 fix ℓ ( 𝑘 ) = ℓ + 𝑉 ∫ Ω 𝑘 𝐿 ℓ ( ˆ k · ˆ n ) | 𝛿 ( k )| . (2.7)(Kaiser 1987; Hamilton 1992; Feldman et al. 1994). This necessarily incurs an error at O( 𝜃 max ) where 𝜃 max is the survey opening angle, and is valid only for the smallest surveys.At the next order in approximation is the ‘Yamamoto formalism’ used by most current estimators; this approximates the LoS as theposition vector of a single galaxy, i.e. ˆ n = ˆ r or ˆ n = ˆ r , as shown in Fig. 1. In this case, one may expand 𝐿 ℓ ( ˆ k · ˆ n ) via the addition theorem One may remove the window function by judicious use of quadratic estimators (e.g., Philcox 2020). In full, 𝛿 ( k ) should be multiplied by a compensation function 𝑀 ( k ) to account for the mass assignment scheme window function. Note this depends on the survey opening angle, since the power spectrum is an integral over all galaxy pairs. Its importance increases as the wavenumberbecomes small.MNRAS000
Comparison of the various choices of lines-of-sight (LoS) for the galaxy pair. For each example, the observer is located at the based of the triangle,with the galaxies at r and r with separation vector 𝚫 = r − r . We define 𝜙 as the internal triangle angle with cos 𝜙 = ˆ r · ˆ r . In the small-angle limit, wehave 𝜙 ≈ 𝜃 , where 𝜃 ≡ Δ /| r + r | . The fixed LoS approximation (ˆ n = const . ) incurs an O ( 𝜃 ) error, whilst the Yamamoto approximation (ˆ n = ˆ r ) andmidpoint or bisector estimators (ˆ n ∝ r + r or ˆ r + ˆ r ) incur O ( 𝜃 ) and O ( 𝜃 ) errors respectively, averaging over 𝜃 < 𝜃 max . with the power spectrum as the Fourier transform of the correlation function. Additionally, we use the shorthand ∫ k ≡ ∫ 𝑑 k ( 𝜋 ) , ∫ Ω 𝑘 ≡ ∫ 𝑑 Ω 𝑘 𝜋 . (2.4) The conventional estimator for the galaxy power spectrum multipoles, 𝑃 ℓ ( 𝑘 ) , is defined as a Fourier transform of two density fields 𝛿 ( r ) and 𝛿 ( r ) : ˆ 𝑃 ℓ ( 𝑘 ) = ℓ + 𝑉 ∫ Ω 𝑘 ∫ 𝑑 r 𝑑 r 𝑒 − 𝑖 k ·( r − r ) 𝛿 ( r ) 𝛿 ( r ) 𝐿 ℓ ( ˆ k · ˆ n ) (2.5)(e.g., Hand et al. 2017), where ˆ n is the joint line-of-sight (LoS) to the pair of galaxies at ( r , r ) , 𝐿 ℓ is the Legendre polynomial of order ℓ , 𝑉 is the survey volume and hats denote unit vectors. Whilst we have not included 𝑘 -space binning, this is a straightforward addition, requiringan additional integral over | k | . In spectroscopic surveys, we do not have access to 𝛿 directly, only a set of galaxy and random particle positionswith associated density fields 𝑛 𝑔 ( r ) and 𝑛 𝑟 ( r ) respectively. In this context, we replace 𝛿 ( r ) → 𝑤 ( r )[ 𝑛 𝑔 ( r ) − 𝛼 𝑟 𝑛 𝑟 ( r )] 𝐼 / , 𝐼 ≡ ∫ 𝑑 r 𝑤 ( r ) ¯ 𝑛 ( r ) (2.6)(Yamamoto et al. 2006), where 𝑤 ( r ) are some weights (accounting for systematics, optimality and completeness), ¯ 𝑛 is the mean galaxy densityand 𝛼 𝑟 is the ratio of randoms to galaxies. (2.5) is strictly an estimator for the window-convolved power spectrum, and includes a shot-noiseterm. Since neither effect depends on the LoS choice, we will ignore both henceforth.In the simplest (‘plane-parallel’) approximation, we take ˆ n to be fixed across the survey, thus the estimator simplifies to the form familiarfrom 𝑁 -body estimators, involving the Fourier-space density field 𝛿 ( k ) : ˆ 𝑃 fix ℓ ( 𝑘 ) = ℓ + 𝑉 ∫ Ω 𝑘 𝐿 ℓ ( ˆ k · ˆ n ) | 𝛿 ( k )| . (2.7)(Kaiser 1987; Hamilton 1992; Feldman et al. 1994). This necessarily incurs an error at O( 𝜃 max ) where 𝜃 max is the survey opening angle, and is valid only for the smallest surveys.At the next order in approximation is the ‘Yamamoto formalism’ used by most current estimators; this approximates the LoS as theposition vector of a single galaxy, i.e. ˆ n = ˆ r or ˆ n = ˆ r , as shown in Fig. 1. In this case, one may expand 𝐿 ℓ ( ˆ k · ˆ n ) via the addition theorem One may remove the window function by judicious use of quadratic estimators (e.g., Philcox 2020). In full, 𝛿 ( k ) should be multiplied by a compensation function 𝑀 ( k ) to account for the mass assignment scheme window function. Note this depends on the survey opening angle, since the power spectrum is an integral over all galaxy pairs. Its importance increases as the wavenumberbecomes small.MNRAS000 , 1–27 (2021)
Philcox & Slepian for Legendre polynomials (A5), and arrive at the estimatorˆ 𝑃 Yama ℓ ( 𝑘 ) = 𝜋𝑉 ∫ Ω 𝑘 ℓ ∑︁ 𝑚 = − ℓ 𝑌 𝑚 ∗ ℓ ( ˆ k ) ∫ 𝑑 r 𝑒 𝑖 k · r 𝛿 ( r ) ∫ 𝑑 r 𝑒 − 𝑖 k · r 𝑌 𝑚ℓ ( ˆ r ) 𝛿 ( r ) (2.8) = 𝜋𝑉 ∫ Ω 𝑘 (cid:34) ℓ ∑︁ 𝑚 = − ℓ 𝑌 𝑚 ∗ ℓ ( ˆ k )F (cid:2) 𝑌 𝑚ℓ 𝛿 (cid:3) ( k ) (cid:35) 𝛿 ∗ ( k ) (Yamamoto et al. 2006; Scoccimarro 2015; Bianchi et al. 2015; Hand et al. 2017, 2018). This incurs an error O( 𝜃 ) in the ℓ > O( 𝜃 max ) part vanishes upon r ↔ r symmetrization), and has been employed in almost all recent analyses (e.g., Beutler et al.2017; Gil-Marín et al. 2017; eBOSS Collaboration et al. 2020). The error incurred is not insignificant however; the BAO scale (much belowthe characteristic survey size) has opening angle 𝑟 𝑑 / 𝑑 𝐴 ( 𝑧 eff ) ∼ . 𝑑 𝐴 ( 𝑧 eff ) is the angular diameter distance to the meansurvey redshift and 𝑟 𝑑 is the sound horizon scale at decoupling. At low redshifts, this scale is at the percent level; of importance for futuresurveys such as Euclid and DESI. Considering the whole survey, 𝜃 max is significantly larger, which will affect measurements particularly onlarge scales.Beyond the Yamamoto approximation, there are multiple options for how to proceed. Clearly, setting the LoS to the direction of justone galaxy not an optimal strategy on large scales, and a full treatment would include the positions of both galaxies, since their pairwisevelocity cannot be simply described by a single LoS. This results in a higher-dimensional data-vector however, thus is generally disfavored.Two primary options exist for defining a single LoS correct to O( 𝜃 ) , both of which are shown in Fig. 1; (a) using the position vector of the midpoint of the two galaxies, ˆ n ∝ ( r + r ) (cf. Scoccimarro 2015; Samushia et al. 2015; Bianchi et al. 2015), or (b) using the bisector of thegalaxy-observer-galaxy triangle, ˆ n ∝ ( ˆ r + ˆ r ) (cf. Szalay et al. 1998; Matsubara 2000; Szapudi 2004; Yoo & Seljak 2015). These differ onlyat higher-order, and we will consider both in this work. In full, these are given byˆ 𝑃 midpoint ℓ ( 𝑘 ) = ℓ + 𝑉 ∫ Ω 𝑘 ∫ 𝑑 r 𝑑 r 𝑒 − 𝑖 k ·( r − r ) 𝛿 ( r ) 𝛿 ( r ) 𝐿 ℓ ( ˆ k · (cid:155) r + r ) (2.9)ˆ 𝑃 bisector ℓ ( 𝑘 ) = ℓ + 𝑉 ∫ Ω 𝑘 ∫ 𝑑 r 𝑑 r 𝑒 − 𝑖 k ·( r − r ) 𝛿 ( r ) 𝛿 ( r ) 𝐿 ℓ ( ˆ k · (cid:155) ˆ r + ˆ r ) . Notably, neither straightforwardly factorizes into pieces depending only on r and r , making its computation more involved than that ofthe Yamamoto estimator. A naïve implementation would involve counting all pairs of galaxies individually; this results in an estimator with O( 𝑁 ) complexity for 𝑁 galaxies. In most cases, this is prohibitively slow (though shown to be remarkably efficient on small scales in Philcox& Eisenstein 2020 and Philcox 2021), but, as shown below, separable forms can be obtained using convergent series expansions. Similar estimators may be derived for the multipoles of the two-point correlation function, 𝜉 ℓ ( 𝑟 ) . Analogous to (2.5), the general form isgiven by ˆ 𝜉 ℓ ( 𝑟 ) = ℓ + 𝑉 ∫ 𝑑 r 𝑑 r 𝛿 ( r ) 𝛿 ( r ) 𝐿 ℓ ( ˆ 𝚫 · ˆ n ) (cid:20) 𝛿 D ( 𝑟 − Δ ) 𝜋𝑟 (cid:21) , (2.10)where ˆ n is again the LoS and 𝚫 = r − r is the separation vector. Here, the square brackets pick out a particular value of the galaxy separation Δ = 𝑟 , and, if we substitute (cid:20) 𝛿 D ( 𝑟 − Δ ) 𝜋𝑟 (cid:21) → 𝑣 𝑎 Θ 𝑎 ( 𝑟 ) , (2.11)where Θ 𝑎 ( 𝑟 ) is some binning function with volume 𝑣 𝑎 , (2.10) becomes the estimator for the 2PCF in a finite bin 𝑎 . Just as for the powerspectrum, we cannot access the overdensity field 𝛿 directly, and must instead work with galaxies and random particle catalogs. Conventionally,this leads to the 2PCF being computed via the Landy-Szalay estimator (Landy & Szalay 1993), using 𝐷𝐷 , 𝐷𝑅 and 𝑅𝑅 counts, each ofwhich is estimated via (2.10). Since such complexities do not depend on our choice of ˆ n , we ignore them here, alongside the intricacies ofedge-correction.Two types of correlation function algorithms abound in the literature. Firstly, they are often computed by exhaustive pair counting,scaling as 𝑁 (e.g., Sinha & Garrison 2020). Since we explicitly consider each pair of galaxies, any pairwise LoS can be simply included,rendering moot the analysis of this work. However, direct pair-counting can be exceedingly slow for large-datasets, thus it is commonplace to Specifically, due to the symmetry under r ↔ r permutations for even ℓ , any odd 𝜃 contribution must vanish, thus the difference between midpoint andbisector formalisms starts at O ( 𝜃 ) . MNRAS , 1–27 (2021) eyond Yamamoto Correlators compute the 2PCF via FFT-based approaches. In the plane-parallel approximation of fixed ˆ n , the estimator may be writtenˆ 𝜉 fix ℓ ( 𝑟 ) = ℓ + 𝑉 ∫ 𝑑 𝚫 (cid:20) 𝛿 D ( 𝑟 − Δ ) 𝜋𝑟 (cid:21) 𝐿 ℓ ( ˆ 𝚫 · ˆ n ) ∫ 𝑑 r 𝛿 ( r ) 𝛿 ( r + 𝚫 ) (2.12) = ℓ + 𝑉 ∫ 𝑑 𝚫 (cid:20) 𝛿 D ( 𝑟 − Δ ) 𝜋𝑟 (cid:21) 𝐿 ℓ ( ˆ 𝚫 · ˆ n ) F − (cid:104) | 𝛿 ( k )| (cid:105) ( 𝚫 ) , where we switch variables to r = r , 𝚫 = r − r in the first line, and write the r integral as a convolution in the second. Following the inverseFourier transform, (2.12) is a weighted real-space summation, which is easy to compute. This again incurs an error at O( 𝜃 ) , where 𝜃 max is now the characteristic opening of galaxy pairs separated by the maximum value of 𝑟 considered.Similarly to (2.8), the 2PCF may be estimated in the Yamamoto formalism by expanding the Legendre polynomial in spherical harmonicsand writing the result as a convolution integral:ˆ 𝜉 Yama ℓ ( 𝑟 ) = 𝜋𝑉 ∫ 𝑑 𝚫 (cid:20) 𝛿 D ( 𝑟 − Δ ) 𝜋𝑟 (cid:21) ℓ ∑︁ 𝑚 = − ℓ 𝑌 𝑚 ∗ ℓ ( ˆ 𝚫 ) ∫ 𝑑 r 𝛿 ( r ) 𝛿 ( r + 𝚫 ) 𝑌 𝑚ℓ ( (cid:154) r + 𝚫 ) (2.13) = 𝜋𝑉 ∫ 𝑑 𝚫 (cid:20) 𝛿 D ( 𝑟 − Δ ) 𝜋𝑟 (cid:21) ℓ ∑︁ 𝑚 = − ℓ 𝑌 𝑚 ∗ ℓ ( ˆ 𝚫 )F − (cid:2) 𝛿 ∗ ( k )F (cid:2) 𝑌 𝑚ℓ 𝛿 (cid:3) ( k ) (cid:3) ( 𝚫 ) , (e.g., Slepian & Eisenstein 2016), which is computable in a similar manner to the plane-parallel estimator, now requiring 2 ( ℓ + ) FFTs fora given ℓ . It again suffers an O( 𝜃 ) error.The O( 𝜃 ) pairwise 2PCF estimates can be written in an analogous form to (2.10):ˆ 𝜉 midpoint ℓ ( 𝑟 ) = ℓ + 𝑉 ∫ 𝑑 r 𝑑 r 𝛿 ( r ) 𝛿 ( r ) 𝐿 ℓ ( ˆ 𝚫 · (cid:155) r + r ) (cid:20) 𝛿 D ( 𝑟 − Δ ) 𝜋𝑟 (cid:21) (2.14)ˆ 𝜉 bisector ℓ ( 𝑟 ) = ℓ + 𝑉 ∫ 𝑑 r 𝑑 r 𝛿 ( r ) 𝛿 ( r ) 𝐿 ℓ ( ˆ 𝚫 · (cid:155) ˆ r + ˆ r ) (cid:20) 𝛿 D ( 𝑟 − Δ ) 𝜋𝑟 (cid:21) . Unlike the Yamamoto estimator, neither of these can be straightforwardly recast as a convolution, and thus evaluated via FFTs. It is howeverpossible via series expansions, as will be discussed below.
To obtain an efficient power spectrum algorithm within the midpoint formalism, it is necessary to perform a series expansion on the angulardependence, 𝑌 𝑚ℓ ( (cid:155) r + r ) , such that the estimator can be recast in a form conducive to FFT application. To motivate this, we start with thedefinition (2.9) in terms of the new variables 𝚫 ≡ r − r and R = r + r , using the addition theorem (A5) to write the Legendre polynomialin terms of spherical harmonics: ˆ 𝑃 pair ℓ ( 𝑘 ) = 𝜋𝑉 ℓ ∑︁ 𝑚 = − ℓ ∫ Ω 𝑘 𝑌 𝑚 ∗ ℓ ( ˆ k ) ∫ 𝑑 r 𝑑 𝚫 𝑒 − 𝑖 k · 𝚫 𝛿 ( r ) 𝛿 ( r + 𝚫 ) 𝑌 𝑚ℓ ( ˆ R ) . (3.1)If the spherical harmonic factor 𝑌 𝑚ℓ ( ˆ R ) can be written in a form separable in 𝚫 and r , the above expression can be evaluated asa convolution, allowing for acceleration by way of Fourier transforms (analogous to the 2PCF manipulations in (2.12) & (2.13)). Such anexpansion is indeed possible, via the spherical harmonic shift theorem , which states 𝑌 𝑚ℓ ( (cid:154) a + A ) = ∞ ∑︁ 𝛼 = 𝛼 ∑︁ 𝐽 = ℓ + 𝛼 ∑︁ 𝐽 = max ( ,ℓ − 𝛼 ) 𝐽 ∑︁ 𝑀 = − 𝐽 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 (cid:16) 𝑎𝐴 (cid:17) 𝛼 𝑌 𝑀𝐽 ( ˆ a ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ A ) (3.2)(B7), for arbitrary vectors a , A with 𝑎 (cid:28) 𝐴 . This is proved in Appendix B and is a major new result of this work. Essentially, (3.2) is aninfinite expansion in terms of two spherical harmonics and the (small) ratio of | a | and | A | , giving an arbitrarily accurate approximation of 𝑌 𝑚ℓ ( (cid:154) a + A ) if truncated at sufficiently large 𝛼 . This uses the numerical coefficients 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 given in (B8), which may be pre-computed andobey the relations 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 = 𝐽 + 𝛼 or 𝐽 + 𝛼 + ℓ is odd , (3.3) 𝜑 ,ℓ𝑚𝐽 𝐽 𝑀 = √ 𝜋𝛿 K 𝐽 𝛿 K 𝐽 ℓ 𝛿 K 𝑀 , 𝜑 𝛼, 𝐽 𝐽 𝑀 = √ 𝜋𝛿 K 𝛼 𝛿 K 𝐽 𝛿 K 𝐽 𝛿 K 𝑀 , (where 𝛿 𝐾𝑖 𝑗 is the Kronecker delta, equal to unity if 𝑖 = 𝑗 and zero else), as proved in Appendix B2. We note that it is formally possible to first perform the Ω 𝑘 integral analytically, which leads to the replacement ∫ Ω 𝑘 𝑒 − 𝑖 k · 𝚫 𝐿 ℓ ( ˆ k · ˆ R ) → (− 𝑖 ) ℓ 𝑗 ℓ ( 𝑘 Δ ) 𝐿 ℓ ( 𝚫 · R ) using (A11) & (A12). However, this is not exact when the finite nature of the 𝑘 -space grid is considered, thus will not be adopted here.MNRAS000
To obtain an efficient power spectrum algorithm within the midpoint formalism, it is necessary to perform a series expansion on the angulardependence, 𝑌 𝑚ℓ ( (cid:155) r + r ) , such that the estimator can be recast in a form conducive to FFT application. To motivate this, we start with thedefinition (2.9) in terms of the new variables 𝚫 ≡ r − r and R = r + r , using the addition theorem (A5) to write the Legendre polynomialin terms of spherical harmonics: ˆ 𝑃 pair ℓ ( 𝑘 ) = 𝜋𝑉 ℓ ∑︁ 𝑚 = − ℓ ∫ Ω 𝑘 𝑌 𝑚 ∗ ℓ ( ˆ k ) ∫ 𝑑 r 𝑑 𝚫 𝑒 − 𝑖 k · 𝚫 𝛿 ( r ) 𝛿 ( r + 𝚫 ) 𝑌 𝑚ℓ ( ˆ R ) . (3.1)If the spherical harmonic factor 𝑌 𝑚ℓ ( ˆ R ) can be written in a form separable in 𝚫 and r , the above expression can be evaluated asa convolution, allowing for acceleration by way of Fourier transforms (analogous to the 2PCF manipulations in (2.12) & (2.13)). Such anexpansion is indeed possible, via the spherical harmonic shift theorem , which states 𝑌 𝑚ℓ ( (cid:154) a + A ) = ∞ ∑︁ 𝛼 = 𝛼 ∑︁ 𝐽 = ℓ + 𝛼 ∑︁ 𝐽 = max ( ,ℓ − 𝛼 ) 𝐽 ∑︁ 𝑀 = − 𝐽 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 (cid:16) 𝑎𝐴 (cid:17) 𝛼 𝑌 𝑀𝐽 ( ˆ a ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ A ) (3.2)(B7), for arbitrary vectors a , A with 𝑎 (cid:28) 𝐴 . This is proved in Appendix B and is a major new result of this work. Essentially, (3.2) is aninfinite expansion in terms of two spherical harmonics and the (small) ratio of | a | and | A | , giving an arbitrarily accurate approximation of 𝑌 𝑚ℓ ( (cid:154) a + A ) if truncated at sufficiently large 𝛼 . This uses the numerical coefficients 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 given in (B8), which may be pre-computed andobey the relations 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 = 𝐽 + 𝛼 or 𝐽 + 𝛼 + ℓ is odd , (3.3) 𝜑 ,ℓ𝑚𝐽 𝐽 𝑀 = √ 𝜋𝛿 K 𝐽 𝛿 K 𝐽 ℓ 𝛿 K 𝑀 , 𝜑 𝛼, 𝐽 𝐽 𝑀 = √ 𝜋𝛿 K 𝛼 𝛿 K 𝐽 𝛿 K 𝐽 𝛿 K 𝑀 , (where 𝛿 𝐾𝑖 𝑗 is the Kronecker delta, equal to unity if 𝑖 = 𝑗 and zero else), as proved in Appendix B2. We note that it is formally possible to first perform the Ω 𝑘 integral analytically, which leads to the replacement ∫ Ω 𝑘 𝑒 − 𝑖 k · 𝚫 𝐿 ℓ ( ˆ k · ˆ R ) → (− 𝑖 ) ℓ 𝑗 ℓ ( 𝑘 Δ ) 𝐿 ℓ ( 𝚫 · R ) using (A11) & (A12). However, this is not exact when the finite nature of the 𝑘 -space grid is considered, thus will not be adopted here.MNRAS000 , 1–27 (2021) Philcox & Slepian
In our context, we may use (3.2) to expand 𝑌 𝑚ℓ ( ˆ R ) using R = r + 𝚫 or R = r − 𝚫 . Suppressing summation limits for clarity, theselead to 𝑌 𝑚ℓ ( ˆ R ) = ∑︁ 𝛼𝐽 𝐽 𝑀 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 (cid:18) Δ 𝑟 (cid:19) 𝛼 𝑌 𝑀𝐽 ( ˆ 𝚫 ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) = ∑︁ 𝛼𝐽 𝐽 𝑀 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 (− ) 𝐽 (cid:18) Δ 𝑟 (cid:19) 𝛼 𝑌 𝑀𝐽 ( ˆ 𝚫 ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) , (3.4)using 𝑌 𝑀𝐽 (− ˆ 𝚫 ) = (− ) 𝐽 𝑌 𝑀𝐽 ( ˆ 𝚫 ) . The two may be combined to give the symmetrized form: 𝑌 𝑚ℓ ( ˆ R ) = ∑︁ 𝛼𝐽 𝐽 𝑀 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 𝑌 𝑀𝐽 ( ˆ 𝚫 ) (cid:104) 𝜖 𝛼 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) + (− ) 𝐽 𝜖 𝛼 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) (cid:105) , (3.5)where 𝜖 𝑖 ≡ Δ /( 𝑟 𝑖 ) . Each term is fully separable in ˆ 𝚫 and ˆ r 𝑖 , and the expansion is simply a power series in 𝜖 𝑖 , which is closely related to thepair opening angle 𝜃 ≡ Δ / 𝑅 . Indeed, including all terms up to 𝛼 = 𝐾 gives an approximation incurring an error only at O( 𝜃 𝐾 + ) . At lowestorder ( 𝛼 = 𝜑 ,ℓ𝑚𝐽 𝐽 𝑀 = 𝛿 K 𝐽 𝛿 K 𝐽 ℓ 𝛿 K 𝑀 × √ 𝜋 (3.3), thus 𝑌 𝑚ℓ ( ˆ R ) → (cid:104) 𝑌 𝑚ℓ ( ˆ r ) + (− ) ℓ 𝑌 𝑚ℓ ( ˆ r ) (cid:105) ; (3.6)this simply yields the Yamamoto estimator in symmetrized form.As an example, we consider the expansion of 𝑌 ( ˆ R ) , including all terms up to 𝛼 =
2. From the above expressions, we obtain2 𝑌 ( ˆ R ) = (cid:110) 𝑌 ( ˆ r ) (3.7) + √︂ 𝜋 𝜖 (cid:104) √ 𝑌 ( ˆ 𝚫 ) 𝑌 ( ˆ r ) + √ 𝑌 ( ˆ 𝚫 ) 𝑌 ( ˆ r ) + √ 𝑌 ( ˆ 𝚫 ) 𝑌 ( ˆ r ) − √ 𝑌 ( ˆ 𝚫 ) 𝑌 ( ˆ r ) + √ 𝑌 − ( ˆ 𝚫 ) 𝑌 ( ˆ r ) (cid:105) + √ 𝜋 𝜖 (cid:104) 𝑌 ( ˆ 𝚫 ) 𝑌 ( ˆ r ) − √ 𝑌 ( ˆ 𝚫 ) 𝑌 ( ˆ r ) − 𝑌 ( ˆ 𝚫 ) 𝑌 ( ˆ r ) − 𝑌 ( ˆ 𝚫 ) 𝑌 ( ˆ r ) − √ 𝑌 ( ˆ 𝚫 ) 𝑌 ( ˆ r )+ √ 𝑌 ( ˆ 𝚫 ) 𝑌 , − ( ˆ r ) + √ 𝑌 − ( ˆ 𝚫 ) 𝑌 ( ˆ r ) + √ 𝑌 ( ˆ 𝚫 ) 𝑌 − ( ˆ r ) + √ 𝑌 ( ˆ 𝚫 ) 𝑌 ( ˆ r ) − √ 𝑌 − ( ˆ 𝚫 ) 𝑌 ( ˆ r )+ √ 𝑌 − ( ˆ 𝚫 ) 𝑌 ( ˆ r ) (cid:105) +O( 𝜖 ) (cid:111) + ( r ↔ r ) , with zeroth-, first- and second-order pieces shown in red, orange and green. Whilst this is lengthy, it is nonetheless computationally tractable. The above series expansion may be used to write the pairwise power spectrum estimator as a series of convolutions. Inserting (3.5) into (3.1)gives ˆ 𝑃 midpoint ℓ ( 𝑘 ) = 𝜋𝑉 ℓ ∑︁ 𝑚 = − ℓ ∫ Ω 𝑘 𝑌 𝑚 ∗ ℓ ( ˆ k ) ∑︁ 𝛼𝐽 𝐽 𝑀 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 ∫ 𝑑 𝚫 𝑒 − 𝑖 k · 𝚫 𝑌 𝑀𝐽 ( ˆ 𝚫 ) Δ 𝛼 (3.8) × (cid:34)∫ 𝑑 r 𝛿 ( r ) 𝛿 ( r + 𝚫 ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ r )( 𝑟 ) 𝛼 + (− ) 𝐽 ∫ 𝑑 r 𝛿 ( r − 𝚫 ) 𝛿 ( r ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ r )( 𝑟 ) 𝛼 (cid:35) . Relabelling 𝚫 → − 𝚫 , r → r and k → − k in the second term shows that the two are equivalent up to a factor (− ) ℓ . Here and henceforthwe will assume even ℓ , allowing us to absorb the symmetrization: ˆ 𝑃 midpoint ℓ ( 𝑘 ) = 𝜋𝑉 ℓ ∑︁ 𝑚 = − ℓ ∫ Ω 𝑘 𝑌 𝑚 ∗ ℓ ( ˆ k ) ∑︁ 𝛼𝐽 𝐽 𝑀 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 (− ) 𝐽 ∫ 𝑑 𝚫 𝑒 − 𝑖 k · 𝚫 𝑌 𝑀𝐽 ( ˆ 𝚫 ) Δ 𝛼 ∫ 𝑑 r 𝛿 ( r ) 𝛿 ( r − 𝚫 ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ r )( 𝑟 ) 𝛼 . (3.9)Next, we note that the r integral can be written as a convolution: 𝑔 𝛼𝐽 ( 𝑚 − 𝑀 ) ( 𝚫 ) ≡ (− ) 𝐽 ∫ 𝑑 r (cid:104) 𝛿 ( r )( 𝑟 ) − 𝛼 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) (cid:105) 𝛿 ( r − 𝚫 ) (3.10) = (− ) 𝐽 F − (cid:104) 𝛿 ∗ ( k )F (cid:104) ( 𝑟 ) − 𝛼 𝑌 𝑚 − 𝑀𝐽 𝛿 (cid:105) ( k ) (cid:105) ( 𝚫 ) , An alternative approach would be to find an expansion of 𝑌 𝑚ℓ ( ˆ R ) that is separable in r and r rather than r 𝑖 and 𝚫 . In this case, the power spectrum couldbe computed in the same manner as in the Yamamoto approximation, however, this is more difficult to obtain since 𝑟 / 𝑟 is order unity, so one cannot simplyapply (3.2). An approach similar to this will prove useful for the bisector formalism however (§5). Note the distinction between 𝜃 , the opening angle for a particular pair, and 𝜃 max , the characteristic opening angle, which is fixed for a particular analysis. We replace (− ) 𝐽 with (− ) 𝐽 for later convenience; for even ℓ , these must have the same sign, due to the parity-rules on 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 .MNRAS , 1–27 (2021) eyond Yamamoto Correlators and the 𝚫 integral is then just a Fourier transform:ˆ 𝑃 midpoint ℓ ( 𝑘 ) = 𝜋𝑉 ℓ ∑︁ 𝑚 = − ℓ ∫ Ω 𝑘 𝑌 𝑚 ∗ ℓ ( ˆ k ) ∑︁ 𝛼𝐽 𝐽 𝑀 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 F (cid:104) 𝑌 𝑀𝐽 ( ˆ 𝚫 ) Δ 𝛼 𝑔 𝛼𝐽 ( 𝑚 − 𝑀 ) ( 𝚫 ) (cid:105) ( k ) . (3.11)Calculation of the spectra is thus reduced to computing a convolution for each { 𝛼, 𝐽 , ( 𝑚 − 𝑀 )} triplet, then performing a Fourier transformfor each ( ℓ𝑚 ) pair. We note that the summation over { 𝛼, 𝐽 , 𝐽 } can be moved inside the k -space Fourier transform; we separate it here forclarity. For 𝛼 =
0, we require 𝐽 = 𝐽 = ℓ , 𝑀 = 𝑔 ℓ𝑚 ( 𝚫 ) = (− ) ℓ F − (cid:2) 𝛿 ∗ ( k )F (cid:2) 𝑌 𝑚ℓ 𝛿 (cid:3) ( k ) (cid:3) ( 𝚫 ) , (3.12)and, since 𝑌 𝑀𝐽 ( ˆ 𝚫 ) → ( 𝜋 ) − / , the estimator is equal to that of Hand et al. (2017), as expected.When implementing (3.11), we must be aware of a certain subtlety. Our formalism requires the Fourier transform of a function of 𝚫 weighted by Δ 𝛼 where 𝛼 ≥
0. If 𝑔 𝛼𝐽 ( 𝑚 − 𝑀 ) ( 𝚫 ) decays slower than Δ − 𝛼 , we will obtain an integrand that, in the limit of infinite surveyvolume, is not square integrable, i.e. it diverges at large Δ . This is particularly clear when the mean survey distance, ¯ 𝑟 is large; in this case Δ /( 𝑟 𝑖 ) ≈ Δ /( 𝑟 ) , which is a monotonically increasing function of Δ . Whilst the infinite volume limit is somewhat academic (since ourexpansion parameter 𝜃 max cannot be assumed small), for finite volumes, we note that the magnitude of Δ 𝛼 𝑔 𝛼𝐽 ( 𝑚 − 𝑀 ) increases towards thesurvey edges for 𝛼 >
0. As such, it is important that we consider the full extent of the 𝚫 -space function. For a survey of characteristic width 𝐿 ,the convolved width is 2 𝐿 , thus this requires us to use a grid at least twice the survey width when painting particles. This restriction reducesthe efficiency of the algorithm, since we must double the number of grid cells per dimension to obtain the same Nyquist frequency.The selection rules on { 𝐽 , 𝐽 , 𝑀 } allow us to quantify the method’s complexity. In general, if one wishes to compute all even powerspectrum multipoles up to ℓ max using even (odd) 𝛼 = 𝛼 , we must compute all 𝑔 𝛼𝐽 (cid:48) 𝑀 (cid:48) functions with even (odd) 𝐽 (cid:48) up to 𝐽 (cid:48) = ℓ max + 𝛼 ; a totalof [ + ( ℓ max + 𝛼 )/ ] , each of which requires two Fourier transforms. The Fourier transform over 𝚫 can then be performed just once per ( ℓ, 𝑚 ) pair ( i.e. [ + ℓ max / ] times). As a concrete example, computing the spectra up to ℓ max = 𝑔 𝛼𝐽 (cid:48) 𝑀 (cid:48) coefficients for 𝛼 = 𝛼 = Whilst the above estimator is mathematically valid, closer inspection reveals a curious property; it contains terms both odd and even inthe small angles 𝜖 𝑖 . Since (for even ℓ ) the power-spectrum definition is symmetric under permutation of r and r , we would expect anycontributions of O( 𝜃 𝐾 ) to vanish for odd 𝐾 (recalling 𝜃 ≡ Δ / 𝑅 ). In fact, this is the case, and the above expansion can be recast in a mannerto make this manifest. Below, we consider a straightforward way to achieve this, based on iterated infinite sequences. An alternative method,which is less obvious a priori , but simpler to implement, is described in Appendix D.In order to demonstrate that the odd terms in 𝜃 vanish, we first consider the term (cid:104) 𝜖 𝛼 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) + (− ) 𝐽 𝜖 𝛼 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) (cid:105) (3.13)appearing in (3.5). At lowest order in 𝜃 , 𝜖 ≈ 𝜖 ≈ 𝜃 and r ≈ r ≈ r , giving (cid:104) 𝜖 𝛼 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) + (− ) 𝐽 𝜖 𝛼 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) (cid:105) = 𝜃 𝛼 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) (cid:104) + (− ) 𝐽 (cid:105) + O( 𝜃 𝛼 + ) . (3.14)For even 𝐽 , the contribution is O( 𝜃 𝛼 ) , yet for odd 𝐽 , the term starts only at O( 𝜃 𝛼 + ) . Since 𝐽 + 𝛼 must be even for 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 to be non-zero,we find that even 𝛼 terms in (3.5) ( i.e. those with even powers of 𝜖 𝑖 ) begin to contribute at O( 𝜃 𝛼 ) , whilst the leading-order piece of odd 𝛼 terms vanishes, and their contribution starts at O( 𝜃 𝛼 + ) . This is not sufficient to demonstrate that there are no terms odd in 𝜃 however, since,the term containing 𝜖 𝑖 could include a non-cancelling 𝜖 𝑖 contribution for example.To obtain a manifestly parity-even expansion, we first split the 𝛼 summation of (3.5) into even and odd pieces, noting that the lattercontributions must start at O( 𝜃 𝛼 + ) , viz. the above discussion: 𝑌 𝑚ℓ ( ˆ R ) = ∑︁ even 𝛼 ∑︁ 𝐽 𝐽 𝑀 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 𝑌 𝑀𝐽 ( ˆ 𝚫 ) (cid:104) 𝜖 𝛼 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) + 𝜖 𝛼 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) (cid:105) (3.15) + ∑︁ odd 𝛼 ∑︁ 𝐽 𝐽 𝑀 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 𝑌 𝑀𝐽 ( ˆ 𝚫 ) (cid:104) 𝜖 𝛼 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) − 𝜖 𝛼 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) (cid:105) . Next, we rewrite the odd-parity piece by expanding the r term around 𝚫 = , which allows us to explicitly cancel the lowest-order piece. Todo so, we require a generalized version of the spherical harmonic shift theorem, proved in appendix C: (cid:18) | a || A + a | (cid:19) 𝛼 𝑌 𝑀𝐽 ( ˆ a ) 𝑌 𝑚 − 𝑀𝐽 ( (cid:154) a + A ) = ∑︁ 𝛽𝑆 𝑆 𝑇 Ω 𝛽,𝛼𝐽 𝐽 𝑚𝑀𝑆 𝑆 𝑇 (cid:16) 𝑎𝐴 (cid:17) 𝛼 + 𝛽 𝑌 𝑇𝑆 ( ˆ a ) 𝑌 𝑚 − 𝑇𝑆 ( ˆ A ) , (3.16)where the Ω coefficients are given in (C3), and the 𝑆 𝑖 summation is limited to the range [ max ( , 𝐽 𝑖 − 𝛽 ) , 𝐽 𝑖 + 𝛽 ] . Notably, Ω ,𝛼𝐽 𝐽 𝑚𝑀𝑆 𝑆 𝑇 = 𝛿 K 𝑆 𝐽 𝛿 K 𝑆 𝐽 𝛿 K 𝑇 𝑀 from (C4). Applying this to 𝜖 𝛼 𝑌 𝑀𝐽 ( ˆ 𝚫 ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) with a = 𝚫 , A = r yields 𝜖 𝛼 𝑌 𝑀𝐽 ( ˆ 𝚫 ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) = 𝜖 𝛼 𝑌 𝑀𝐽 ( ˆ 𝚫 ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) + ∑︁ 𝛽> ∑︁ 𝑆 𝑆 𝑇 𝛽 Ω 𝛽,𝛼𝐽 𝐽 𝑚𝑀𝑆 𝑆 𝑇 𝜖 𝛼 + 𝛽 𝑌 𝑇𝑆 ( ˆ 𝚫 ) 𝑌 𝑚 − 𝑇𝑆 ( ˆ r ) , (3.17) MNRAS , 1–27 (2021)
Philcox & Slepian separating out the 𝛽 = Δ / 𝑟 𝑖 ≈ 𝜃 , thus the radius of convergence is somewhat diminished compared tothe all-parity form. Inserting this in the parity-odd part of (3.15) and symmetrizing over r ↔ r gives 𝑌 𝑚ℓ ( ˆ R ) (cid:12)(cid:12) odd = − ∑︁ odd 𝛼 ∑︁ 𝐽 𝐽 𝑀 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 ∑︁ 𝛽> ∑︁ 𝑆 𝑆 𝑇 𝛽 Ω 𝛽,𝛼𝐽 𝐽 𝑚𝑀𝑆 𝑆 𝑇 𝑌 𝑇𝑆 ( ˆ 𝚫 ) (cid:104) 𝜖 𝛼 + 𝛽 𝑌 𝑚 − 𝑇𝑆 ( ˆ r ) − (− ) 𝐽 + 𝑆 𝜖 𝛼 + 𝛽 𝑌 𝑚 − 𝑇𝑆 ( ˆ r ) (cid:105) . (3.18)Importantly the 𝛽 = 𝛼 contains terms only starting at O( 𝜖 𝛼 + 𝑖 ) . We have therefore succeeded inremoving the asymmetric piece at lowest order. In practice, this means that we have removed all terms proportional to 𝜖 𝑖 , reducing the totalof convolutions that need to be performed, since there is no longer a requirement to compute 𝛼 = 𝜖 𝑖 . This is possible via a similar prescription,first noting that the lowest-order piece of the summand in (3.18) ( i.e. that with r = r ) is non-zero only for odd 𝛽 ; a consequence of the parityrules given in Appendix C, restricting 𝛽 + 𝐽 + 𝑆 to be even. For even 𝛽 , we again have a cancellation between terms involving r and r .As before, we may expand the relevant r piece in terms of r , and cancel the lowest-order contribution. This procedure of “split accordingto parity, expand, symmetrize” may be iterated, and leads to the following form: 𝑌 𝑚ℓ ( ˆ R ) = ∑︁ even 𝛼 ∑︁ 𝐽 𝐽 𝑀 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 𝑌 𝑀𝐽 ( ˆ 𝚫 ) (cid:104) 𝜖 𝛼 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) + 𝜖 𝛼 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) (cid:105) (3.19) − ∑︁ odd 𝛼 ∑︁ 𝐽 𝐽 𝑀 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 ∑︁ odd 𝛽 ∑︁ 𝑆 𝑆 𝑇 𝛽 Ω 𝛽,𝛼𝐽 𝐽 𝑚𝑀𝑆 𝑆 𝑇 𝑌 𝑇𝑆 ( ˆ 𝚫 ) (cid:104) 𝜖 𝛼 + 𝛽 𝑌 𝑚 − 𝑇𝑆 ( ˆ r ) + 𝜖 𝛼 + 𝛽 𝑌 𝑚 − 𝑇𝑆 ( ˆ r ) (cid:105) + ∑︁ odd 𝛼 ∑︁ 𝐽 𝐽 𝑀 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 ∑︁ even 𝛽> ∑︁ 𝑆 𝑆 𝑇 𝛽 Ω 𝛽,𝛼𝐽 𝐽 𝑚𝑀𝑆 𝑆 𝑇 ∑︁ odd 𝛾 ∑︁ 𝑄 𝑄 𝑅 𝛾 Ω 𝛾, ( 𝛼 + 𝛽 ) 𝑆 𝑆 𝑚𝑇𝑄 𝑄 𝑅 𝑌 𝑅𝑄 ( ˆ 𝚫 )× (cid:104) 𝜖 𝛼 + 𝛽 + 𝛾 𝑌 𝑚 − 𝑅𝑄 ( ˆ r ) + 𝜖 𝛼 + 𝛽 + 𝛾 𝑌 𝑚 − 𝑅𝑄 ( ˆ r ) (cid:105) + ... , which contains only even powers of 𝜖 𝑖 . Whilst each successive line requires more work to evaluate, we note that the only terms up to thesecond (third) line are required for an expansion correct to third (fifth) order in 𝜃 . Furthermore, we may simplify the above by introducingredefined coefficients Φ 𝛼,ℓ𝑚𝐽 𝐽 𝑀 , such that 𝑌 𝑚ℓ ( ˆ R ) = ∑︁ even 𝛼 ∑︁ 𝐽 𝐽 𝑀 Φ 𝛼,ℓ𝑚𝐽 𝐽 𝑀 𝑌 𝑀𝐽 ( ˆ 𝚫 ) (cid:104) 𝜖 𝛼 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) + 𝜖 𝛼 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) (cid:105) , (3.20)for even 𝐽 , 𝐽 , 𝛼 , using Φ 𝛼,ℓ𝑚𝐽 𝐽 𝑀 = 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 (3.21) − ∑︁ odd 𝛽<𝛼 ∑︁ 𝐽 (cid:48) 𝐽 (cid:48) 𝑀 (cid:48) 𝛼 − 𝛽 𝜑 𝛽,ℓ𝑚𝐽 (cid:48) 𝐽 (cid:48) 𝑀 (cid:48) Ω ( 𝛼 − 𝛽 ) ,𝛽𝐽 (cid:48) 𝐽 (cid:48) 𝑚𝑀 (cid:48) 𝐽 𝐽 𝑀 + ∑︁ odd 𝛽<𝛼 ∑︁ 𝐽 (cid:48) 𝐽 (cid:48) 𝑀 (cid:48) 𝜑 𝛽,ℓ𝑚𝐽 (cid:48) 𝐽 (cid:48) 𝑀 (cid:48) 𝛼 − 𝛽 ∑︁ even 𝛾<𝛼 − 𝛽,𝛾> ∑︁ 𝐽 (cid:48)(cid:48) 𝐽 (cid:48)(cid:48) 𝑀 (cid:48)(cid:48) Ω 𝛾,𝛽𝐽 (cid:48) 𝐽 (cid:48) 𝑚𝑀 (cid:48) 𝐽 (cid:48)(cid:48) 𝐽 (cid:48)(cid:48) 𝑀 (cid:48)(cid:48) Ω ( 𝛼 − 𝛽 − 𝛾 ) , ( 𝛽 + 𝛾 ) 𝐽 (cid:48)(cid:48) 𝐽 (cid:48)(cid:48) 𝑚𝑀 (cid:48)(cid:48) 𝐽 𝐽 𝑀 + ... . As an example, for ℓ = 𝑚 =
1, the parity-even expansion truncating at 𝛼 = 𝑌 ( ˆ R ) = (cid:110) 𝑌 ( ˆ r ) (3.22) − √ 𝜋 𝜖 (cid:104) 𝑌 ( ˆ 𝚫 ) 𝑌 ( ˆ r ) − √ 𝑌 ( ˆ 𝚫 ) 𝑌 ( ˆ r ) − 𝑌 ( ˆ 𝚫 ) 𝑌 ( ˆ r ) − 𝑌 ( ˆ 𝚫 ) 𝑌 ( ˆ r ) − √ 𝑌 ( ˆ 𝚫 ) 𝑌 ( ˆ r )+ √ 𝑌 ( ˆ 𝚫 ) 𝑌 , − ( ˆ r ) + √ 𝑌 − ( ˆ 𝚫 ) 𝑌 ( ˆ r ) + √ 𝑌 ( ˆ 𝚫 ) 𝑌 − ( ˆ r ) + √ 𝑌 ( ˆ 𝚫 ) 𝑌 ( ˆ r ) − √ 𝑌 − ( ˆ 𝚫 ) 𝑌 ( ˆ r )+ √ 𝑌 − ( ˆ 𝚫 ) 𝑌 ( ˆ r ) (cid:105) + O( 𝜖 ) (cid:111) + ( r ↔ r ) , marking zeroth- (second-)order terms in red (green). At fixed maximum order in 𝜖 𝑖 , this contains significantly fewer terms than (3.7), sinceall odd powers vanish. We note that Φ ,ℓ𝑚𝐽 𝐽 𝑀 = − 𝜑 ,ℓ𝑚𝐽 𝐽 𝑀 here. Adopting this notation, the full parity-even power spectrum estimator is given byˆ 𝑃 midpoint ℓ ( 𝑘 ) = 𝜋𝑉 ℓ ∑︁ 𝑚 = − ℓ ∫ Ω 𝑘 𝑌 𝑚 ∗ ℓ ( ˆ k ) ∑︁ even 𝛼 ∑︁ 𝐽 𝐽 𝑀 Φ 𝛼,ℓ𝑚𝐽 𝐽 𝑀 ∫ 𝑑 𝚫 𝑒 − 𝑖 k · 𝚫 𝑌 𝑀𝐽 ( ˆ 𝚫 ) Δ 𝛼 ∫ 𝑑 r 𝛿 ( r ) 𝛿 ( r + 𝚫 ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ r )( 𝑟 ) 𝛼 This arises naturally in the alternate derivation given in Appendix D. MNRAS , 1–27 (2021) eyond Yamamoto Correlators ⇒ ˆ 𝑃 midpoint ℓ ( 𝑘 ) = 𝜋𝑉 ℓ ∑︁ 𝑚 = − ℓ ∫ Ω 𝑘 𝑌 𝑚 ∗ ℓ ( ˆ k ) ∑︁ even 𝛼 ∑︁ 𝐽 𝐽 𝑀 Φ 𝛼,ℓ𝑚𝐽 𝐽 𝑀 F (cid:104) 𝑌 𝑀𝐽 ( ˆ 𝚫 ) Δ 𝛼 𝑔 𝛼𝐽 ( 𝑚 − 𝑀 ) ( 𝚫 ) (cid:105) ( k ) , (3.23)analogous to (3.9), but now requiring only even 𝛼 (and thus even 𝐽 , 𝐽 ), significantly reducing the necessary number of 𝑔 𝛼𝐽 (cid:48) 𝑀 (cid:48) fields (3.10).For ℓ max =
2, we require 28 (45) functions for 𝛼 = 𝛼 = 𝛼 ( i.e. those that contribute to even 𝛼 < 𝛼 max ), and the expansion formally requires 2 𝜖 𝑖, max ≈ 𝜃 max (cid:28)
1, rather than 𝜖 𝑖, max ≈ 𝜃 max (cid:28)
1. If the latter conditions are met, bothforms are fully convergent.
Just as for the power spectrum, we may construct an efficient 2PCF estimator via series expansions coupled with FFTs. In this instance, theangular dependence appears through 𝐿 ℓ ( ˆ 𝚫 · ˆ R ) ; our goal therefore is to expand this in a form separable in 𝚫 and r . To this end, we use the Legendre polynomial shift theorem , which states 𝐿 ℓ ( ˆ a · (cid:154) a + A ) = ∞ ∑︁ 𝛼 = ℓ + 𝛼 ∑︁ 𝐽 = max ( ℓ − 𝛼, ) 𝐽 + ℓ + 𝑓 𝛼,ℓ𝐽 (cid:16) 𝑎𝐴 (cid:17) 𝛼 𝐿 𝐽 ( ˆ a · ˆ A ) (4.1)(cf. E4), for 𝑎 (cid:28) 𝐴 . This is proved in Appendix E1, and uses the coefficients 𝑓 𝛼,ℓ𝐽 presented in (E11). In brief, the derivation proceeds bynoting that 𝐿 ℓ ( ˆ a · (cid:154) a + A ) can be written as a sum of spherical harmonics in ˆ a and (cid:154) a + A using the addition theorem (A5), then expanding 𝑌 𝑚ℓ ( (cid:154) a + A ) using the spherical harmonic shift theorem of Appendix B, and simplifying the resulting coefficients. An alternative approachwould be to write the Legendre polynomial as a power series then perform a Taylor expansion; this gives the same results.In our context, we set a = 𝚫 , A = r or a = − 𝚫 , A = r in (4.1) to yield the symmetrized form 𝐿 ℓ ( ˆ 𝚫 · ˆ R ) = ∞ ∑︁ 𝛼 = ℓ + 𝛼 ∑︁ 𝐽 = max ( ℓ − 𝛼, ) 𝐽 + ℓ + 𝑓 𝛼.ℓ𝐽 (cid:104) 𝜖 𝛼 𝐿 𝐽 ( 𝜇 ) + (− ) ℓ + 𝐽 𝜖 𝛼 𝐿 𝐽 ( 𝜇 ) (cid:105) , (4.2)where 𝜇 𝑖 = ˆ 𝚫 · ˆ r 𝑖 and 𝜖 𝑖 = Δ /( 𝑟 𝑖 ) . This is markedly simpler than the result for the spherical harmonic shift theorem, and can equivalentlybe derived by an expansion of 𝐿 ℓ ( ˆ 𝚫 · ˆ R ) in powers of ˆ 𝚫 · ˆ R , which is then expanded via the binomial theorem. Key properties of thisexpansion are discussed in Appendix E1; here we note that the coefficients are non-zero only for even 𝛼 + ℓ + 𝐽 and 𝑓 ,ℓ𝐽 = 𝛿 K 𝐽ℓ , such that 𝐿 ℓ ( ˆ 𝚫 · ˆ R ) = 𝐿 ℓ ( 𝜇 ) = 𝐿 ℓ (− 𝜇 ) at lowest order, i.e. the Yamamoto approximation.As an example, we consider the expansion of the ℓ = 𝜖 𝑖 (and hence the opening angle 𝜃 ):2 𝐿 ( ˆ 𝚫 · ˆ R ) = 𝐿 ( 𝜇 ) + 𝜖 [ 𝐿 ( 𝜇 ) − 𝐿 ( 𝜇 )] + 𝜖 [ 𝐿 ( 𝜇 ) − 𝐿 ( 𝜇 ) + 𝐿 ( 𝜇 )] (4.3) − 𝜖 [ 𝐿 ( 𝜇 ) − 𝐿 ( 𝜇 ) + 𝐿 ( 𝜇 )] + 𝜖 [ 𝐿 ( 𝜇 ) + 𝐿 ( 𝜇 ) − 𝐿 ( 𝜇 ) + 𝐿 ( 𝜇 )]+ O (cid:16) 𝜖 (cid:17) + ( r ↔ − r ) , where colors separate the different orders, as before. For even (odd) 𝛼 , the expansion simply consists of Legendre polynomials of even (odd)order up to ℓ + 𝛼 , weighted by powers of 𝜖 𝑖 . Inserting (4.2) into the 2PCF estimator (2.14) yieldsˆ 𝜉 midpoint ℓ ( 𝑟 ) = 𝑉 ∑︁ 𝛼𝐽 ( 𝐽 + ) 𝑓 𝛼,ℓ𝐽 ∫ 𝑑 r 𝑑 r 𝛿 ( r ) 𝛿 ( r ) (cid:18) Δ 𝑟 (cid:19) 𝛼 𝐿 𝐽 (− ˆ 𝚫 · ˆ r ) (cid:20) 𝛿 D ( 𝑟 − Δ ) 𝜋𝑟 (cid:21) , (4.4)suppressing summation indices for clarity and absorbing the r ↔ − r symmetrization (for even ℓ ). Expanding the Legendre polynomial viathe addition theorem (A5), we may express the integrand in separable form:ˆ 𝜉 midpoint ℓ ( 𝑟 ) = 𝜋𝑉 ∑︁ 𝛼𝐽 𝐽 ∑︁ 𝑀 = − 𝐽 𝑓 𝛼,ℓ𝐽 ∫ 𝑑 𝚫 𝑌 𝑀 ∗ 𝐽 ( ˆ 𝚫 ) Δ 𝛼 (cid:20) 𝛿 D ( 𝑟 − Δ ) 𝜋𝑟 (cid:21) ∫ 𝑑 r 𝛿 ( r ) 𝛿 ( r − 𝚫 ) 𝑌 𝑀𝐽 (− ˆ r )( 𝑟 ) 𝛼 . (4.5)As for the power spectrum, the r integral is simply a convolution: 𝑔 𝛼𝐽 𝑀 ( 𝚫 ) = (− ) 𝐽 ∫ 𝑑 r (cid:104) 𝛿 ( r )( 𝑟 ) − 𝛼 𝑌 𝑀𝐽 ( ˆ r ) (cid:105) 𝛿 ( r − 𝚫 ) (4.6) = (− ) 𝐽 F − (cid:104) 𝛿 ∗ ( k )F (cid:104) ( 𝑟 ) − 𝛼 𝑌 𝑀𝐽 𝛿 (cid:105) ( k ) (cid:105) ( 𝚫 ) , MNRAS000
Just as for the power spectrum, we may construct an efficient 2PCF estimator via series expansions coupled with FFTs. In this instance, theangular dependence appears through 𝐿 ℓ ( ˆ 𝚫 · ˆ R ) ; our goal therefore is to expand this in a form separable in 𝚫 and r . To this end, we use the Legendre polynomial shift theorem , which states 𝐿 ℓ ( ˆ a · (cid:154) a + A ) = ∞ ∑︁ 𝛼 = ℓ + 𝛼 ∑︁ 𝐽 = max ( ℓ − 𝛼, ) 𝐽 + ℓ + 𝑓 𝛼,ℓ𝐽 (cid:16) 𝑎𝐴 (cid:17) 𝛼 𝐿 𝐽 ( ˆ a · ˆ A ) (4.1)(cf. E4), for 𝑎 (cid:28) 𝐴 . This is proved in Appendix E1, and uses the coefficients 𝑓 𝛼,ℓ𝐽 presented in (E11). In brief, the derivation proceeds bynoting that 𝐿 ℓ ( ˆ a · (cid:154) a + A ) can be written as a sum of spherical harmonics in ˆ a and (cid:154) a + A using the addition theorem (A5), then expanding 𝑌 𝑚ℓ ( (cid:154) a + A ) using the spherical harmonic shift theorem of Appendix B, and simplifying the resulting coefficients. An alternative approachwould be to write the Legendre polynomial as a power series then perform a Taylor expansion; this gives the same results.In our context, we set a = 𝚫 , A = r or a = − 𝚫 , A = r in (4.1) to yield the symmetrized form 𝐿 ℓ ( ˆ 𝚫 · ˆ R ) = ∞ ∑︁ 𝛼 = ℓ + 𝛼 ∑︁ 𝐽 = max ( ℓ − 𝛼, ) 𝐽 + ℓ + 𝑓 𝛼.ℓ𝐽 (cid:104) 𝜖 𝛼 𝐿 𝐽 ( 𝜇 ) + (− ) ℓ + 𝐽 𝜖 𝛼 𝐿 𝐽 ( 𝜇 ) (cid:105) , (4.2)where 𝜇 𝑖 = ˆ 𝚫 · ˆ r 𝑖 and 𝜖 𝑖 = Δ /( 𝑟 𝑖 ) . This is markedly simpler than the result for the spherical harmonic shift theorem, and can equivalentlybe derived by an expansion of 𝐿 ℓ ( ˆ 𝚫 · ˆ R ) in powers of ˆ 𝚫 · ˆ R , which is then expanded via the binomial theorem. Key properties of thisexpansion are discussed in Appendix E1; here we note that the coefficients are non-zero only for even 𝛼 + ℓ + 𝐽 and 𝑓 ,ℓ𝐽 = 𝛿 K 𝐽ℓ , such that 𝐿 ℓ ( ˆ 𝚫 · ˆ R ) = 𝐿 ℓ ( 𝜇 ) = 𝐿 ℓ (− 𝜇 ) at lowest order, i.e. the Yamamoto approximation.As an example, we consider the expansion of the ℓ = 𝜖 𝑖 (and hence the opening angle 𝜃 ):2 𝐿 ( ˆ 𝚫 · ˆ R ) = 𝐿 ( 𝜇 ) + 𝜖 [ 𝐿 ( 𝜇 ) − 𝐿 ( 𝜇 )] + 𝜖 [ 𝐿 ( 𝜇 ) − 𝐿 ( 𝜇 ) + 𝐿 ( 𝜇 )] (4.3) − 𝜖 [ 𝐿 ( 𝜇 ) − 𝐿 ( 𝜇 ) + 𝐿 ( 𝜇 )] + 𝜖 [ 𝐿 ( 𝜇 ) + 𝐿 ( 𝜇 ) − 𝐿 ( 𝜇 ) + 𝐿 ( 𝜇 )]+ O (cid:16) 𝜖 (cid:17) + ( r ↔ − r ) , where colors separate the different orders, as before. For even (odd) 𝛼 , the expansion simply consists of Legendre polynomials of even (odd)order up to ℓ + 𝛼 , weighted by powers of 𝜖 𝑖 . Inserting (4.2) into the 2PCF estimator (2.14) yieldsˆ 𝜉 midpoint ℓ ( 𝑟 ) = 𝑉 ∑︁ 𝛼𝐽 ( 𝐽 + ) 𝑓 𝛼,ℓ𝐽 ∫ 𝑑 r 𝑑 r 𝛿 ( r ) 𝛿 ( r ) (cid:18) Δ 𝑟 (cid:19) 𝛼 𝐿 𝐽 (− ˆ 𝚫 · ˆ r ) (cid:20) 𝛿 D ( 𝑟 − Δ ) 𝜋𝑟 (cid:21) , (4.4)suppressing summation indices for clarity and absorbing the r ↔ − r symmetrization (for even ℓ ). Expanding the Legendre polynomial viathe addition theorem (A5), we may express the integrand in separable form:ˆ 𝜉 midpoint ℓ ( 𝑟 ) = 𝜋𝑉 ∑︁ 𝛼𝐽 𝐽 ∑︁ 𝑀 = − 𝐽 𝑓 𝛼,ℓ𝐽 ∫ 𝑑 𝚫 𝑌 𝑀 ∗ 𝐽 ( ˆ 𝚫 ) Δ 𝛼 (cid:20) 𝛿 D ( 𝑟 − Δ ) 𝜋𝑟 (cid:21) ∫ 𝑑 r 𝛿 ( r ) 𝛿 ( r − 𝚫 ) 𝑌 𝑀𝐽 (− ˆ r )( 𝑟 ) 𝛼 . (4.5)As for the power spectrum, the r integral is simply a convolution: 𝑔 𝛼𝐽 𝑀 ( 𝚫 ) = (− ) 𝐽 ∫ 𝑑 r (cid:104) 𝛿 ( r )( 𝑟 ) − 𝛼 𝑌 𝑀𝐽 ( ˆ r ) (cid:105) 𝛿 ( r − 𝚫 ) (4.6) = (− ) 𝐽 F − (cid:104) 𝛿 ∗ ( k )F (cid:104) ( 𝑟 ) − 𝛼 𝑌 𝑀𝐽 𝛿 (cid:105) ( k ) (cid:105) ( 𝚫 ) , MNRAS000 , 1–27 (2021) Philcox & Slepian where the 𝑔 𝛼𝐽 𝑀 ( 𝚫 ) functions are the same as those in (3.10). In this notation, the 2PCF estimator is given byˆ 𝜉 midpoint ℓ ( 𝑟 ) = 𝜋𝑉 ∑︁ 𝛼𝐽 𝑀 𝑓 𝛼,ℓ𝐽 ∫ 𝑑 𝚫 𝑌 𝑀 ∗ 𝐽 ( ˆ 𝚫 ) Δ 𝛼 𝑔 𝛼𝐽 𝑀 ( 𝚫 ) (cid:20) 𝛿 D ( 𝑟 − Δ ) 𝜋𝑟 (cid:21) , (4.7)which may be computed via a simple summation over real-space pixels, given some 𝑟 = | 𝚫 | binning scheme. This is analogous to the estimatorobtained within the Yamamoto scheme:ˆ 𝜉 Yama ℓ ( 𝑟 ) = 𝜋𝑉 𝑀 = ℓ ∑︁ 𝑀 = − ℓ ∫ 𝑑 𝚫 𝑌 𝑀 ∗ ℓ ( ˆ 𝚫 ) 𝑔 ℓ𝑀 ( 𝚫 ) (cid:20) 𝛿 D ( 𝑟 − Δ ) 𝜋𝑟 (cid:21) , (4.8)and simply requires additional 2PCF multipoles to be computed, weighted by powers of Δ /( 𝑟 ) . Since algorithm this does not require aFourier transform in 𝚫 space, we do not need to use an increased boxsize, contrary to the power spectrum case (cf. §3.2). Instead, we requireonly that the width of the box is at least the maximum galaxy dimension plus the maximum separation of interest, to avoid particle overlapon periodic wrapping. To obtain a 2PCF estimate in multipoles up to ℓ max including O( 𝜃 𝐾 max ) corrections, we need to compute the 𝑔 𝛼𝐽 𝑀 functions with 𝛼 ≤ 𝐾 , 𝐽 ≤ ℓ max + 𝛼 , | 𝑀 | ≤ 𝐽 , each of which requires a forward and inverse Fourier transform. For ℓ max = 𝑔 𝛼𝐽 𝑀 coefficients for 𝛼 = 𝛼 =
2) just as for ˆ 𝑃 ℓ ( 𝑘 ) , each of which is computed in O( 𝑁 g log 𝑁 g ) time for 𝑁 g grid points. We may further simplify the estimator by writing it as a sum over only even 𝛼 (and thus even multipoles 𝐽 ). This is possible via recasting(4.2) in a manifestly parity-even form: 𝐿 ℓ ( ˆ 𝚫 · ˆ R ) = ∞ ∑︁ even 𝛼 = ℓ + 𝛼 ∑︁ 𝐽 = max ( ℓ − 𝛼, ) 𝐹 𝛼.ℓ𝐽 (cid:104) 𝜖 𝛼 𝐿 𝐽 ( 𝜇 ) + (− ) ℓ + 𝐽 𝜖 𝛼 𝐿 𝐽 ( 𝜇 ) (cid:105) (4.9)(E14), where the corresponding coefficients, 𝐹 𝛼,ℓ𝐽 , are given in (E15). The corresponding estimator is thusˆ 𝜉 midpoint ℓ ( 𝑟 ) = 𝜋𝑉 ∑︁ even 𝛼 ∑︁ 𝐽 𝑀 𝐹 𝛼,ℓ𝐽 ∫ 𝑑 𝚫 𝑌 𝑀 ∗ 𝐽 ( ˆ 𝚫 ) Δ 𝛼 𝑔 𝛼𝐽 𝑀 ( 𝚫 ) (cid:20) 𝛿 D ( 𝑟 − Δ ) 𝜋𝑟 (cid:21) . (4.10)The derivation of (4.9) is analogous to that of §3.3 and sketched in Appendix E2. As an example, the even-parity expansions of ℓ = ℓ = 𝐿 ( ˆ 𝚫 · ˆ R ) = 𝐿 ( 𝜇 ) + 𝜖 [− 𝐿 ( 𝜇 ) + 𝐿 ( 𝜇 ) − 𝐿 ( 𝜇 )] (4.11) + 𝜖 [ 𝐿 ( 𝜇 ) + 𝐿 ( 𝜇 ) − 𝐿 ( 𝜇 ) + 𝐿 ( 𝜇 )] + O (cid:16) 𝜖 (cid:17) + ( r ↔ − r ) 𝐿 ( ˆ 𝚫 · ˆ R ) = 𝐿 ( 𝜇 ) − 𝜖 [ 𝐿 ( 𝜇 ) − 𝐿 ( 𝜇 ) + 𝐿 ( 𝜇 )]− 𝜖 [ 𝐿 ( 𝜇 ) + 𝐿 ( 𝜇 ) − 𝐿 ( 𝜇 ) + 𝐿 ( 𝜇 ) − 𝐿 ( 𝜇 )] + O (cid:16) 𝜖 (cid:17) + ( r ↔ − r ) , correct to fifth-order in 𝜃 . These results are consistent with the second-order calculation of Slepian & Eisenstein (2015, Eqs. 19 & 20), andmay be extended to arbitrarily high order in ℓ and 𝛼 , with each order just requiring computation of more 𝑔 𝛼𝐽 𝑀 ( 𝚫 ) functions. We now consider the estimators for the power spectrum in the bisector formalism. The discussion in this section closely follows that ofCastorina & White (2018, Appendix E), but is extended to give closed-form estimators at arbitrary order in 𝜃 max . To derive the midpoint power spectrum estimator (§3), we began by expanding the spherical harmonic 𝑌 𝑚ℓ ( ˆ n ) as a separable function in ˆ r andˆ 𝚫 . This was both tractable, since one can write ˆ n ∝ ( r + 𝚫 ) , and favorable, since it led to an explicit expansion in powers of 𝜖 = Δ /( 𝑟 ) ≈ 𝜃 .When adopting the bisector LoS, we have ˆ n ∝ ˆ r + ˆ r , which does not facilitate as straightforward separation into 𝚫 and r pieces due to theadditional normalization by | r | and | r | . For this reason, we take a different approach, proposed by Castorina & White (2018), whereuponwe first expand instead in r and r , then in the small angle 𝜃 .Whilst the explicit details of this can be found in Appendix F, for now we note that the derivation proceeds by writing the bisector vector MNRAS , 1–27 (2021) eyond Yamamoto Correlators d = 𝑡 r + ( − 𝑡 ) r , where 𝑡 = 𝑟 /( 𝑟 + 𝑟 ) . Expanding the spherical harmonic 𝑌 𝑚ℓ ( ˆ d ) in terms of 𝑡 r and ( − 𝑡 ) r via the solid harmonicaddition theorem (A10) gives the finite series 𝑌 𝑚ℓ ( ˆ d ) = ℓ ∑︁ 𝜆 = 𝜆 ∑︁ 𝜇 = − 𝜆 𝐴 𝜆𝜇ℓ𝑚 (cid:16) 𝑡𝑟 𝑑 (cid:17) ℓ 𝑌 𝜇𝜆 ( ˆ r ) 𝑌 𝑚 − 𝜇ℓ − 𝜆 ( ˆ r ) , (5.1)with 𝑑 = | d | and coupling coefficient 𝐴 𝜆𝜇ℓ𝑚 defined in (B1). Due to the ( 𝑡𝑟 / 𝑑 ) − ℓ factor, this is not directly separable in r , r ; however, wecan write 𝑡𝑟 𝑑 = [ ( + cos 𝜙 )] − / , (5.2)and expand in the small parameter ( − cos 𝜙 ) , where cos 𝜙 = ˆ r · ˆ r . From the geometry of Fig. 1, it is clear that 𝜃 ≈ 𝜙 / Δ (cid:28) 𝑟 , 𝑟 , thusthis expansion is effectively one in 𝜃 (since 1 − cos 𝜙 ≈ 𝜃 ). Omitting lengthy algebra, this leads to the series expansion for 𝑌 𝑚ℓ ( ˆ d ) : 𝑌 𝑚ℓ ( ˆ d ) = ∞ ∑︁ 𝑘 = ℓ + 𝑘 ∑︁ 𝐽 = ℓ + 𝑘 ∑︁ 𝐽 = 𝐽 ∑︁ 𝑀 = − 𝐽 𝐵 𝑘,ℓ𝑚𝐽 𝐽 𝑀 𝑌 𝑀𝐽 ( ˆ r ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) (5.3)(cf. F7), where the coupling coefficients 𝐵 𝑘,ℓ𝑚𝐽 𝐽 𝑀 are defined in (F8). These obey the properties 𝐵 𝑘,ℓ𝑚𝐽 𝐽 𝑀 = ℓ + 𝐽 + 𝐽 is odd (5.4) 𝐵 ,ℓ𝑚𝐽 𝐽 𝑀 = − ℓ 𝐴 𝐽 𝑀ℓ𝑚 𝛿 K 𝐽 ( ℓ − 𝐽 ) , 𝐵 𝑘, 𝐽 𝐽 𝑀 = √ 𝜋𝛿 K 𝐽 𝛿 K 𝐽 𝛿 K 𝑀 (cf. Appendix F), with the latter indicating that 𝑌 ( ˆ d ) = /√ 𝜋 , as expected.Notably, (5.3) is an infinite expansion involving only separable functions of ˆ r and ˆ r (and no powers of | r 𝑖 | ), which will allow for anefficient power spectrum estimator to be wrought. The approximation order is controlled by the value of 𝑘 ; fixing to 𝑘 ≤ 𝐾 expands theinternal angle up to ( − cos 𝜙 ) 𝐾 , which scales as (√ 𝜃 ) 𝐾 . Unlike in the midpoint formalism, the expansion is automatically symmetric in 𝜙 and hence 𝜃 , thus there is no need for the parity-even manipulations of §3.3. In this case, the symmetries of 𝐵 𝑘,ℓ𝑚𝐽 𝐽 𝑀 require both odd andeven 𝐽 , 𝐽 (albeit with the restriction 𝐽 + 𝐽 = even for even ℓ ), unlike for the midpoint estimator.As an example, we consider the expansion of 𝑌 ( ˆ d ) at 𝑘 = 𝑘 = i.e. the O( 𝜃 ) and O( 𝜃 ) contributions. This yields 𝑌 ( ˆ d ) (cid:12)(cid:12)(cid:12) 𝑘 = = √ 𝜋 (cid:104) 𝑌 ( ˆ r ) 𝑌 ( ˆ r ) + √ 𝑌 ( ˆ r ) 𝑌 ( ˆ r ) (cid:105) + ( r ↔ r ) (5.5) 𝑌 ( ˆ d ) (cid:12)(cid:12)(cid:12) 𝑘 = = √ 𝜋 (cid:104) 𝑌 ( ˆ r ) 𝑌 ( ˆ r ) + √ 𝑌 ( ˆ r ) 𝑌 − ( ˆ r ) + √ 𝑌 ( ˆ r ) 𝑌 ( ˆ r ) − √ 𝑌 ( ˆ r ) 𝑌 ( ˆ r )+ √ 𝑌 ( ˆ r ) 𝑌 ( ˆ r ) + √ 𝑌 ( ˆ r ) 𝑌 − ( ˆ r ) − √ 𝑌 ( ˆ r ) 𝑌 ( ˆ r ) (cid:105) + ( r ↔ r ) . Several points are of note. Firstly, in this formalism we have the same terms appearing at 𝑘 = 𝑘 =
1; for instance there is a 𝑌 ( ˆ r ) 𝑌 ( ˆ r ) term sourced both at 𝑘 = O( 𝜃 𝐾 ) part contains a factor 𝜖 𝐾𝑖 . Secondly, setting 𝑘 = 𝑌 ( ˆ d ) = 𝑌 ( ˆ r ) + 𝑌 ( ˆ r ) ). This is clear from thenon-trivial 𝑘 = 𝛼 = 𝑘 =
1, even though both are O( 𝜃 ) (the order to which the two LoS definitions agree). The difference is due todifferent higher-order terms (and exists also between the midpoint expansion and its even-parity equivalent).Given the simple nature of (5.3) (containing only spherical harmonics in ˆ r and ˆ r ), it is interesting to consider whether the aboveapproach may be applied also to the midpoint estimator. In this case, we use R = r + r , with (5.1) becoming 𝑌 𝑚ℓ ( ˆ R ) = ℓ ∑︁ 𝜆 = 𝜆 ∑︁ 𝜇 = − 𝜆 𝐴 𝜆𝜇ℓ𝑚 (cid:16) 𝑟 𝑅 (cid:17) ℓ (cid:18) 𝑟 𝑟 (cid:19) 𝜆 − ℓ 𝑌 𝜇𝜆 ( ˆ r ) 𝑌 𝑚 − 𝜇ℓ − 𝜆 ( ˆ r ) (5.6) = (cid:104) + 𝑠 + 𝑠 cos 𝜙 (cid:105) − ℓ / ℓ ∑︁ 𝜆 = 𝜆 ∑︁ 𝜇 = − 𝜆 𝐴 𝜆𝜇ℓ𝑚 𝑠 ℓ − 𝜆 𝑌 𝜇𝜆 ( ˆ r ) 𝑌 𝑚 − 𝜇ℓ − 𝜆 ( ˆ r ) , where 𝑠 = 𝑟 / 𝑟 , using the definition 𝑅 = 𝑟 + 𝑟 + 𝑟 𝑟 cos 𝜙 . As before, this can be expanded around the point cos 𝜙 =
1, via (cid:104) + 𝑠 + 𝑠 cos 𝜙 (cid:105) − ℓ / = ∞ ∑︁ 𝑘 = (cid:18) − ℓ / 𝑘 (cid:19) ( + 𝑠 ) − ℓ − 𝑘 ( 𝑠 ) 𝑘 ( cos 𝜙 − ) 𝑘 , (5.7)provided that ( + 𝑠 ) > 𝑠 , i.e. 𝑠 >
0. The difficulty arises from the additional factor of 𝑠 ; we require a separable perturbative expansionof ( + 𝑠 ) − ℓ − 𝑘 yet 𝑠 is order unity. It is this complexity (not found in the bisector approach) that yields the ( 𝚫 , r ) separation of §3.1 moreuseful for the midpoint LoS definition. MNRAS000
0. The difficulty arises from the additional factor of 𝑠 ; we require a separable perturbative expansionof ( + 𝑠 ) − ℓ − 𝑘 yet 𝑠 is order unity. It is this complexity (not found in the bisector approach) that yields the ( 𝚫 , r ) separation of §3.1 moreuseful for the midpoint LoS definition. MNRAS000 , 1–27 (2021) Philcox & Slepian
By inserting the expansion of (5.3) into (3.1), we obtain an estimator for the power spectrum in the bisector formalism:ˆ 𝑃 bisector ℓ ( 𝑘 ) = 𝜋𝑉 ℓ ∑︁ 𝑚 = − ℓ ∫ Ω 𝑘 𝑌 𝑚 ∗ ℓ ( ˆ k ) ∫ 𝑑 r 𝑑 r 𝑒 − 𝑖 k ·( r − r ) 𝛿 ( r ) 𝛿 ( r ) ∑︁ 𝑘𝐽 𝐽 𝑀 𝐵 𝑘,ℓ𝑚𝐽 𝐽 𝑀 𝑌 𝑀𝐽 ( ˆ r ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) (5.8) = 𝜋𝑉 ℓ ∑︁ 𝑚 = − ℓ ∫ Ω 𝑘 𝑌 𝑚 ∗ ℓ ( ˆ k ) ∑︁ 𝑘𝐽 𝐽 𝑀 𝐵 𝑘,ℓ𝑚𝐽 𝐽 𝑀 (cid:18)∫ 𝑑 r 𝑒 𝑖 k · r 𝛿 ( r ) 𝑌 𝑀𝐽 ( ˆ r ) (cid:19) (cid:18)∫ 𝑑 r 𝛿 ( r ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) (cid:19) ⇒ ˆ 𝑃 bisector ℓ ( 𝑘 ) = 𝜋𝑉 ℓ ∑︁ 𝑚 = − ℓ ∫ Ω 𝑘 𝑌 𝑚 ∗ ℓ ( ˆ k ) ∑︁ 𝑘𝐽 𝐽 𝑀 𝐵 𝑘,ℓ𝑚𝐽 𝐽 𝑀 F (cid:104) 𝑌 𝑀𝐽 𝛿 (cid:105) (− k )F (cid:104) 𝑌 𝑚 − 𝑀𝐽 𝛿 (cid:105) ( k ) . In the last line, we have written the estimator in terms of the Fourier transforms of 𝑌 𝑀𝐽 𝛿 , resulting in a simple-to-implement estimator. Thisbears strong similarities to the Yamamoto approximation (but now involves two density-weighted spherical harmonics), and does not includeany factors of 𝜖 𝑖 = Δ /( 𝑟 𝑖 ) . In comparison to the midpoint estimator, the lack of 𝜖 𝑖 terms has the advantage that we do not need to use a widerbox-size relative to the standard estimators (cf. §3.2).Practically, one estimates the power spectra by first computing F (cid:2) 𝑌 𝑀𝐽 𝛿 (cid:3) ( k ) for all (odd and even) 𝐽 up to ℓ max + 𝑘 max , giving a totalof ( + ℓ max + 𝑘 max )( + ℓ max + 𝑘 max )/ F (cid:2) 𝑌 𝑀𝐽 𝛿 (cid:3) (− k ) = (− ) 𝑀 F (cid:2) 𝑌 − 𝑀𝐽 𝛿 (cid:3) ∗ ( k ) . For ℓ max =
4, this requires 15(21) terms for 𝑘 max = 𝑘 max = ) . For modest computational resources, it may be impractical to store all possible F (cid:2) 𝑌 𝑀𝐽 𝛿 (cid:3) grids, sinceeach requires substantial memory. If these are computed separately, we note that the total number of relevant 𝑌 𝑚 ∗ ℓ 𝑌 𝑀𝐽 𝑌 𝑚 − 𝑀𝐽 combinations issignificant; 20 for 𝑘 = 𝑘 = ℓ =
2. We caution that this approach may thus be quite computationally intensive.
We finally turn to the two-point correlation function in the bisector formalism. In this case, the angular dependence of the estimator (2.14) isgiven by 𝐿 ℓ ( ˆ 𝚫 · ˆ d ) where d is the bisector vector as before. For this, we do not require a specialized perturbative expansion, instead expandingthe Legendre polynomial via the addition theorem (A5) and utilizing the expansion of the previous section (5.3): 𝐿 ℓ ( ˆ 𝚫 · ˆ d ) = 𝜋 ℓ + ℓ ∑︁ 𝑚 = − ℓ 𝑌 𝑚 ∗ ℓ ( ˆ 𝚫 ) ∞ ∑︁ 𝑘 = ℓ + 𝑘 ∑︁ 𝐽 = ℓ + 𝑘 ∑︁ 𝐽 = 𝐽 ∑︁ 𝑀 = − 𝐽 𝐵 𝑘,ℓ𝑚𝐽 𝐽 𝑀 𝑌 𝑀𝐽 ( ˆ r ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) . (6.1)Since we expand 𝑌 𝑚ℓ ( ˆ d ) in terms of ˆ r and ˆ r , contracting with 𝑌 𝑚ℓ ( ˆ 𝚫 ) does not simplify the coefficients (unlike that seen for the midpointcase in §4.1). Whilst this form is appealing since it involves no new derivations or powers of Δ /( 𝑟 𝑖 ) , we note that, at least at leading order, itis possible to obtain a series expansion of 𝐿 ℓ ( ˆ 𝚫 · ˆ d ) in powers of 𝜖 𝑖 and 𝐿 ℓ ( ˆ 𝚫 · ˆ r 𝑖 ) as in the midpoint formalism. This can be done by writingthe Legendre polynomial as a power series then performing a Taylor expansion (as in Slepian & Eisenstein 2015).To obtain an efficient 2PCF estimator, we insert (6.1) into (2.14), givingˆ 𝜉 bisector ℓ ( 𝑟 ) = 𝜋𝑉 ℓ ∑︁ 𝑚 = − ℓ ∫ 𝑑 𝚫 𝑌 𝑚 ∗ ℓ ( ˆ 𝚫 ) (cid:20) 𝛿 D ( 𝑟 − Δ ) 𝜋𝑟 (cid:21) ∑︁ 𝑘𝐽 𝐽 𝑀 𝐵 𝑘,ℓ𝑚𝐽 𝐽 𝑀 ∫ 𝑑 r (cid:104) 𝑌 𝑀𝐽 ( ˆ r ) 𝛿 ( r ) (cid:105) (cid:104) 𝑌 𝑚 − 𝑀𝐽 ( (cid:155) r + 𝚫 ) 𝛿 ( r + 𝚫 ) (cid:105) , (6.2)writing r = r + 𝚫 . As in §2, the r integral is a convolution and can thus be computed via Fourier methods; specifically ℎ 𝑚𝐽 𝐽 𝑀 ( 𝚫 ) ≡ ∫ 𝑑 r (cid:104) 𝑌 𝑀𝐽 ( ˆ r ) 𝛿 ( r ) (cid:105) (cid:104) 𝑌 𝑚 − 𝑀𝐽 ( (cid:155) r + 𝚫 ) 𝛿 ( r + 𝚫 ) (cid:105) (6.3) = F − (cid:104) F (cid:104) 𝑌 𝑀𝐽 𝛿 (cid:105) (− k )F (cid:104) 𝑌 𝑚 − 𝑀𝐽 𝛿 (cid:105) ( k ) (cid:105) ( 𝚫 ) . This gives the full form: ˆ 𝜉 bisector ℓ ( 𝑟 ) = 𝜋𝑉 ℓ ∑︁ 𝑚 = − ℓ ∫ 𝑑 𝚫 𝑌 𝑚 ∗ ℓ ( ˆ 𝚫 ) (cid:20) 𝛿 D ( 𝑟 − Δ ) 𝜋𝑟 (cid:21) ∑︁ 𝑘𝐽 𝐽 𝑀 𝐵 𝑘,ℓ𝑚𝐽 𝐽 𝑀 ℎ 𝑚𝐽 𝐽 𝑀 ( 𝚫 ) . (6.4)Whilst this does not require powers of 𝜖 𝑖 (as in §4), it is quite computationally expensive to implement since we require separate convolutionsfor each { 𝐽 , 𝐽 , 𝑚, 𝑀 } quartet; a total of 20 (68) for 𝑘 = 𝑘 =
1) at ℓ =
2, as before. This may be contrasted from the midpoint case (4.7),which requires estimation of the 𝑔 𝛼𝐽 𝑀 functions, depending only on a single total angular momentum. However, note that the summation over 𝑘, 𝐽 , 𝐽 , 𝑀 can be performed before the inverse transform, leading to significant expedition. MNRAS , 1–27 (2021) eyond Yamamoto Correlators Having established the existence of efficient pairwise power spectrum and 2PCF estimators, we are now ready to implement them. Below,we will consider their application on realistic survey, including a comparison between the two LoS schemes, before first discussing theconvergence of the aforementioned series expansions.
Three series expansions have been introduced in this work; one for each of 𝑌 𝑚ℓ ( ˆ R ) and 𝐿 ℓ ( ˆ 𝚫 · ˆ R ) (§3.1 & §4.1), relevant to the midpointestimators, and one for 𝑌 𝑚ℓ ( ˆ d ) (§5.1) for the bisector estimators. To test these, we consider a simple scenario where we fix R and vary 𝚫 such that we scan over the convergence parameter 𝜃 = Δ / 𝑅 . We compute the three angular statistics using the associated r , r vectors as afunction of 𝜃 , storing both the true value and the approximation at a given value of 𝛼 or 𝑘 . For the midpoint estimators, we consider both thestandard expansions and those after the even-parity transformations of §3.3 & §4.3.The results are shown in Fig. 2 for ℓ = 𝜃 increases (and thus the survey becomes wider), and, for 𝜃 = . 𝜃 (cid:28) 𝜃 (cid:28)
1) and (b) each odd-parity termcontributes to an infinite number of even-parity terms; our formalism only captures their contributions to even terms up to the given 𝛼 max ,and will thus incur a larger error. That being said, we find the even parity expansions (which are cheaper to compute) to converge fairly wellfor 𝜃 (cid:46) .
1. As a point of comparison, the mock data-sets used in §7.2 have 𝜃 ∼ . 𝛼 max = ℓ = .
5% (2 . 𝜃 = . 𝜃 = .
2) or 1 .
8% (7 . ℓ =
4; not an insignificant error!
To provide a practical demonstration of above algorithms, we apply them to a set of realistic mock galaxy catalogs. For this purpose, we use 24MultiDark-patchy (hereafter ‘patchy’) simulations (Rodríguez-Torres et al. 2016; Kitaura et al. 2016), created for the analysis of the twelfthdata-release (DR12) (Alam et al. 2017) of the Baryon Oscillation Spectroscopic Survey (BOSS), part of SDSS-III (Eisenstein et al. 2011;DESI Collaboration et al. 2016). The data is split according to the criteria discussed in (Beutler et al. 2017); here we specialize to the patchwith the largest number density (and volume); the north Galactic cap in the redshift range 0 . < 𝑧 < .
5. This has total volume 1 . ℎ − Gpc and mean redshift 𝑧 = .
38. The simulations are generated with the cosmology { Ω 𝑚 = . , Ω 𝑏 = . , 𝜎 = . , ℎ = . } ,and each contains ∼ . × simulated galaxies, alongside a random catalog 50 × larger.For each simulation, the mock galaxy positions are painted to a cuboidal grid using nbodykit (Hand et al. 2018), using FKP weights(Feldman et al. 1994) and triangle-shaped-cell interpolation, with a fiducial value Ω 𝑚 = .
31 used to convert redshifts and angles into comovingcoordinates. We use the same gridding parameters as in the final BOSS data release but double the box-size (as discussed in §3.2), giving aNyquist frequency 𝑘 Nyq = . ℎ Mpc − . All further computations are carried out in Python, making use of the pyfftw library to implementthe algorithms of §3.2, §4.2, §5.2 & §6. When binning spectra, we adopt the parameters { 𝑘 min = . ℎ Mpc − , 𝑘 max = . ℎ Mpc − , Δ 𝑘 = . ℎ Mpc − } , { 𝑟 min = ℎ − Mpc , 𝑟 max = ℎ − Mpc , Δ 𝑟 = ℎ − Mpc } in Fourier- and configuration-space respectively. In both cases,we use a maximum Legendre multipole of ℓ max = Contributions to the power spectrum multipoles in the midpoint and bisector formalisms are shown in Fig. 3. The figure displays thecontributions as a function of their order in the (assumed small) parameter 𝜃 max , with, for example, the 𝑘 = O( 𝜃 ) . For the midpoint, the same order is obtained by summing the 𝛼 = 𝛼 = 𝛼 = O( 𝜃 ) piece is equal to the Yamamoto spectrum inthe midpoint formalism, but includes higher-order corrections for the bisector LoS definition. Our first note is that, in all cases, the leadingorder contribution is subdominant, indicating that the Yamamoto approximation is, in general, a fair approximation for BOSS. Relative tothe errorbars, any wide-angle corrections represent a few-percent correction at best for this survey. Whilst they seem somewhat smaller for Publicly available at data.sdss.org/sas/dr12/boss/lss. Jupyter notebooks containing our analysis pipeline are publicly available at github.com/oliverphilcox/BeyondYamamoto. Our implementation takes ∼ 𝑃 ℓ ( 𝑘 ) of a single BOSS-like simulation on 4 Intel Skylake processors using 𝛼 max = 𝑘 max =
2, and ℓ max =
4. The analysisrequires ∼
150 GB of memory, though this can be significantly reduced at the expense of longer computation time.MNRAS000
150 GB of memory, though this can be significantly reduced at the expense of longer computation time.MNRAS000 , 1–27 (2021) Philcox & Slepian max F r a c t i o n a l E rr o r Expansion of Y ( R ) max Expansion of L ( R ) k max Expansion of Y ( d ) Midpoint: All ParityMidpoint: Even ParityBisector = 0.20= 0.10= 0.05 max F r a c t i o n a l E rr o r Expansion of Y ( R ) max Expansion of L ( R ) k max Expansion of Y ( d ) Midpoint: All ParityMidpoint: Even ParityBisector = 0.20= 0.10= 0.05
Figure 2.
Convergence plots for the series expansions used in this work. The left, center and right panels show the expansions of 𝑌 𝑚ℓ ( ˆ R ) , 𝐿 ℓ ( ˆ 𝚫 · ˆ R ) and 𝑌 𝑚ℓ ( ˆ d ) respectively (§3.1, §4.1 & §5.1), where ˆ R is the direction vector of the galaxy midpoint, ˆ d is the angle bisector and ˆ 𝚫 is the separation vector. Resultsare shown for both ℓ = ℓ = | 𝐴 true − 𝐴 approx |/| 𝐴 true | for (possibly complex) statistic 𝐴 . Whilst we pick a single value of 𝑚 for the spherical harmonicplots, we have checked that this is representative of all 𝑌 𝑚ℓ functions for the stated ℓ . The horizontal axis shows the approximation order, with 𝛼 = 𝐾 ( 𝑘 = 𝐾 )corresponding to 𝜃 𝐾 ( 𝜃 𝐾 ) corrections, where 𝜃 = Δ / 𝑅 ; the galaxy pair opening angle (with 𝜃 ∼ . 𝜃 max becomes large, the convergence is reduced, and divergences occur. This is more prominent for the even-parity expansions, whichare formally valid only for small 2 𝜃 . For the bisector cases, the 𝛼 = 𝜃 . the bisector LoS, this is primarily due to the additional O( 𝜃 ) terms absorbed into the 𝑘 = 𝛼 = 𝛼 =
2. This may be rationalized by noting that convergence condition in §3.3 isstricter than that of §4.1 (2 𝜃 max (cid:28) 𝜃 max (cid:28) ℎ − Mpc. This increases themean galaxy distance by a factor ∼
3, thus reducing 𝜃 max by the same factor. Resulting power spectrum contributions are shown in Fig. 4, andmatch our expectations; we see a greater difference between the O( 𝜃 ) and O( 𝜃 ) contributions than before, and significantly improvedconvergence for all estimators. Clearly then, our estimators will perform better in regimes where the characteristic angle 𝜃 max is smaller. Evenin this example, the systematic errors are at the percent level (and again artificially smaller for the bisector formalism, due to the resummationof higher orders into the zeroth-order term), which may become important for surveys such as DESI. MNRAS , 1–27 (2021) eyond Yamamoto Correlators k [ h Mpc ] | P ( k ) | [ h M p c ] Midpoint: All Parity k [ h Mpc ]Midpoint: Even Parity k [ h Mpc ]Bisector = 2= 4 Term Term Term k [ h Mpc ] | P |/ P ( k ) Midpoint: All Parity k [ h Mpc ]Midpoint: Even Parity k [ h Mpc ]Bisector = 2= 4 Term Term
Figure 3.
Contributions to the pairwise power spectra from the midpoint and bisector formalisms for a BOSS-like sample. The left (middle) panel shows theall parity (even parity) midpoint estimator of §3.2 (§3.3), whilst the right gives the bisector approach of §5.2, based on Castorina & White (2018). We show theabsolute magnitudes of the relevant terms in the top panel, whilst the bottom gives the ratio to the standard deviation of the Yamamoto spectra. In all cases, weplot the mean of 24 patchy mocks and display the quadrupole (full lines) and hexadecapole (dashed lines), noting that the monopole receives no wide-anglecorrections. The red, green and blue lines correspond to corrections of order 𝜃 , 𝜃 and 𝜃 in the characteristic angle 𝜃 max ( i.e. 𝜃 corresponds to 𝛼 = 𝛼 = 𝛼 = 𝑘 = 𝜃 term is the Yamamoto approximation, whilst for the bisector it is the sum of this and higher-order corrections, as discussed in §5.1. In generalthe corrections are a small fraction of the errorbars, and we find some convergence issues for the midpoint approximations, indicating that the convergencecriterion 𝜃 max (cid:28) Fig. 5 displays the analogous results for the 2PCF multipoles of the (unshifted) patchy mocks. We reiterate that these do not include window-function corrections for simplicity. The 2PCF presents a very different story to the power spectrum. We see a clear demarcation of theexpansion orders, with a high degree of convergence seen for the quadrupole and hexadecapole on all scales tested for both midpoint andbisector estimates. At larger radii the fractional contribution of the higher-order terms increases, and, extrapolating by eye, one would expectthe expansions to break down around 300 − ℎ − Mpc. These results clearly demonstrate that our expansions are performing as expected.As for the power spectrum, the post-Yamamoto corrections are small; ∼
1% of the error bars on all scales, with a slight reduction at highradius (whereupon cosmic variance dominates).The marked difference between the 2PCF and power spectrum begs the question: why does our power spectrum algorithm fare somuch worse? Given that we observe similar behavior in the midpoint and bisector formalisms (with a very different expansion scheme), it isunlikely to stem from some algebraic fault. Our justification for this relies on the following observations: (1) the fractional contribution of thehigher-order terms to the 2PCF increases as a function of 𝑟 , becoming unity on scales comparable with the survey width, and (2) the powerspectrum is an (Bessel function weighted) integral over the 2PCF. The scales on which the small- 𝜃 approximation is most tenuous thus havean impact on all the 𝑘 -modes of the power spectrum. This is clearly seen from the (more convergent) bisector formalism; whilst the ratio ofquartic to quadratic contribution strongly increases as a function of scale for the 2PCF, it is roughly constant for 𝑃 ℓ ( 𝑘 ) . As previously noted,for a convergent series expansion, the power spectrum requires the survey opening angle be small, rather than just that of the maximum radial The method of Slepian & Eisenstein (2016) gives a straightforward approach by which to include these.MNRAS000
1% of the error bars on all scales, with a slight reduction at highradius (whereupon cosmic variance dominates).The marked difference between the 2PCF and power spectrum begs the question: why does our power spectrum algorithm fare somuch worse? Given that we observe similar behavior in the midpoint and bisector formalisms (with a very different expansion scheme), it isunlikely to stem from some algebraic fault. Our justification for this relies on the following observations: (1) the fractional contribution of thehigher-order terms to the 2PCF increases as a function of 𝑟 , becoming unity on scales comparable with the survey width, and (2) the powerspectrum is an (Bessel function weighted) integral over the 2PCF. The scales on which the small- 𝜃 approximation is most tenuous thus havean impact on all the 𝑘 -modes of the power spectrum. This is clearly seen from the (more convergent) bisector formalism; whilst the ratio ofquartic to quadratic contribution strongly increases as a function of scale for the 2PCF, it is roughly constant for 𝑃 ℓ ( 𝑘 ) . As previously noted,for a convergent series expansion, the power spectrum requires the survey opening angle be small, rather than just that of the maximum radial The method of Slepian & Eisenstein (2016) gives a straightforward approach by which to include these.MNRAS000 , 1–27 (2021) Philcox & Slepian k [ h Mpc ] | P ( k ) | [ h M p c ] Midpoint: All Parity k [ h Mpc ]Midpoint: Even Parity k [ h Mpc ]Bisector = 2= 4 Term Term Term k [ h Mpc ] | P |/ P ( k ) Midpoint: All Parity k [ h Mpc ]Midpoint: Even Parity k [ h Mpc ]Bisector = 2= 4 Term Term
Figure 4.
As Fig. 3 but for a BOSS-like galaxy sample shifted radially outwards by an additional 2000 ℎ − Mpc to reduce the 𝜃 max parameter and enforcegreater convergence. In this instance, we use the mean of 5 patchy catalogs. As expected, the higher-order contributions are suppressed relative to Fig. 3, andthe expansion is significantly more convergent (in the sense that the O ( 𝜃 ) term is subdominant to the O ( 𝜃 ) piece). scale considered. This being said, we conclude that the algorithms developed herein are of particular use for the 2PCF, and may be appliedalso to the power spectrum, though the BOSS width lies in the limit of their convergence. Given that the bisector and midpoint are both valid LoS definitions at O( 𝜃 ) it is important to compare their utility in the context ofthe algorithms presented in this work. Below we list several differences, both in terms of efficiency and accuracy of their associated powerspectrum and 2PCFs. • Memory Requirements : A naïve implementation of the bisector formalism requires holding in memory all possible F (cid:2) 𝑌 𝑀𝐽 𝛿 (cid:3) ( k ) fieldsup to a maximum multipole ℓ max + 𝑘 max (such that we can take their outer product); a total of 21 for the O( 𝜃 ) correction with ℓ max = in double precision, this requires 60 GB of memory, which is infeasible on many computer architectures.Whilst one could compute each function ‘on-the-fly’, this will be significantly slower. For the same order of approximation, we require 28 𝑔 𝛼𝐽 𝑀 functions (using even multipoles up to ℓ max + 𝛼 ), but, since no outer product is required, only one needs to be held in memory at anytime. The most efficient CPU-wise implementation of the midpoint algorithm thus requires significantly less memory. • Computation Time : An efficient metric with which to judge the runtime of the two approaches is the number of FFTs required at agiven order. Assuming ℓ max = 𝑃 ℓ ( 𝑘 ) algorithm requires 70 FFTs (two for each 𝑔 𝛼𝐽 𝑀 function, onefor 𝛿 ∗ ( k ) and an addition 13 for transforming into k -space), whilst the bisector approach needs just 28 (one per F (cid:2) 𝑌 𝑀𝐽 𝛿 (cid:3) ( k ) evaluation).The bisector method is thus somewhat faster, though this holds only if high-memory computational resources are available. For the 2PCF, theconclusion is similar, with the midpoint approach requiring 56 FFTs, and the bisector 29 (if memory is no concern). • Convergence : As demonstrated in previous sections, the bisector 𝑃 ℓ ( 𝑘 ) algorithms exhibit stronger convergence than those using themidpoint LoS. This is particularly true for the even-parity approach, which breaks down on BOSS survey scales (also indicating that theYamamoto approximation works poorly there). For the 2PCF, all approaches converge well on scales of interest. • Grid Size : As discussed in §3.2, the midpoint power spectrum algorithm requires the particles to be placed on an FFT grid at least
MNRAS , 1–27 (2021) eyond Yamamoto Correlators
50 100 150 r [ h Mpc] | ( r ) | [ h M p c ] Midpoint: All Parity
50 100 150 r [ h Mpc]Midpoint: Even Parity
50 100 150 r [ h Mpc]Bisector = 2= 4 Term Term Term
50 100 150 r [ h Mpc] ||/ ( r ) Midpoint: All Parity
50 100 150 r [ h Mpc]Midpoint: Even Parity
50 100 150 r [ h Mpc]Bisector = 2= 4 Term Term
Figure 5.
As Fig. 3 but for the two-point correlation function multipoles. In this case, the hierarchy of terms is much clearer, with higher-order effects becomingimportant only at large 𝑟 . For BOSS volumes, the quadratic-order correction is ∼
1% of the error bar for both the monopole and quadrupole, though this willincrease with the survey volume. We find good convergence of both the midpoint and bisector formalism on all scales tested. Note that we do not perform anywindow-function correction on these 2PCF measurements. twice as wide as the survey itself to avoid errors. This requires a finer cell-size for the same Nyquist frequency, and thus slower (and moreexpensive) computaiton. This is not a concern for the bisector formalism or the 2PCF algorithms. • Zeroth-Order Contribution : The midpoint formalisms have the useful property that the zeroth-order ( 𝛼 =
0) contribution is equal tothe well-known Yamamoto approximation. In contrast, the 𝑘 = 𝑌 𝑚 ∗ ℓ ( ˆ k ) function rather than one. • Legendre Expansion : The midpoint 2PCF algorithm can be simply expressed in terms of Legendre multipoles (4.4) unlike the bisectorapproach, which must use spherical harmonics. This is far simpler to interpret, and requires significantly fewer coupling coefficients.Overall, it is clear that both approaches can be implemented in an efficient manner, and give sensible results if the characteristic size isnot too large. The optimal choice of LoS is left to the user, following the above considerations. In general, both approaches agree at O( 𝜃 ) thus this choice is not of particular importance, particularly since any single LoS definition necessarily incurs an O( 𝜃 ) error. This work has presented efficient and practical algorithms for computing the two-point correlators using pairwise lines-of-sight; an extensionto the usual single-particle Yamamoto approximation. Considering both the bisector and midpoint angle definitions, we have shown howconvergent series expansions may be used to write the 2PCF and power spectrum estimators in a form allowing for implementation via FFTs.Any such correction is a function of the characteristic size, 𝜃 max , which we have treated as a perturbation variable. Utilizing newly derivedspherical harmonic and Legendre polynomial shift theorems, we have computed the midpoint corrections at arbitrary order in 𝜃 max , payingclose attention to the existence (and excision) of odd-parity terms. For the bisector LoS definition, the work of Castorina & White (2018) hasbeen extended, including higher-order corrections and an efficient 2PCF algorithm. To demonstrate our approach, the methodology has been An analogous expression for the bisector is possible (and discussed in Slepian & Eisenstein 2015), but requires a different series expansion.MNRAS000
0) contribution is equal tothe well-known Yamamoto approximation. In contrast, the 𝑘 = 𝑌 𝑚 ∗ ℓ ( ˆ k ) function rather than one. • Legendre Expansion : The midpoint 2PCF algorithm can be simply expressed in terms of Legendre multipoles (4.4) unlike the bisectorapproach, which must use spherical harmonics. This is far simpler to interpret, and requires significantly fewer coupling coefficients.Overall, it is clear that both approaches can be implemented in an efficient manner, and give sensible results if the characteristic size isnot too large. The optimal choice of LoS is left to the user, following the above considerations. In general, both approaches agree at O( 𝜃 ) thus this choice is not of particular importance, particularly since any single LoS definition necessarily incurs an O( 𝜃 ) error. This work has presented efficient and practical algorithms for computing the two-point correlators using pairwise lines-of-sight; an extensionto the usual single-particle Yamamoto approximation. Considering both the bisector and midpoint angle definitions, we have shown howconvergent series expansions may be used to write the 2PCF and power spectrum estimators in a form allowing for implementation via FFTs.Any such correction is a function of the characteristic size, 𝜃 max , which we have treated as a perturbation variable. Utilizing newly derivedspherical harmonic and Legendre polynomial shift theorems, we have computed the midpoint corrections at arbitrary order in 𝜃 max , payingclose attention to the existence (and excision) of odd-parity terms. For the bisector LoS definition, the work of Castorina & White (2018) hasbeen extended, including higher-order corrections and an efficient 2PCF algorithm. To demonstrate our approach, the methodology has been An analogous expression for the bisector is possible (and discussed in Slepian & Eisenstein 2015), but requires a different series expansion.MNRAS000 , 1–27 (2021) Philcox & Slepian applied to a set of realistic galaxy catalogs, and shown to give good results for the 2PCF. For the power spectrum, the bisector approach issimilarly effective, though the midpoint algorithms suffer somewhat with convergence issues as the expansion variable becomes large. Suchseries expansions perform best when the angular survey size is relatively small. In the opposing limit, convergence is difficult to achieve, yetalso the Yamamoto approximation itself is poor.For BOSS, the size of the post-Yamamoto effects is small; ∼
1% of the statistical error. For future surveys such as DESI, the errorbarsshrink, though the mean redshift also increases. Following Castorina & White (2018, Fig. 6), we forecast that the estimator-induced wide angleeffects discussed herein can be marginally important, especially on larger scales than considered herein. Since much information regardingprimordial non-Gaussianity is found at low- 𝑘 , implementing estimators beyond the Yamamoto approximation should be seriously considered.As shown above, there is freedom in how best to do this, since both the midpoint and bisector approaches fix the O( 𝜃 ) systematic, butdiffer at O( 𝜃 ) . Here, we find little preference for one over the other; for efficient computation they require somewhat different algorithms,and the bisector approach fares a little better for the power spectrum, yet its memory consumption is significant.Pairwise lines-of-sight are not the perfect solution however. To fully extract information from the large-scale modes, we ought toparametrize the 2PCF by two lines-of-sight (e.g., Pápai & Szapudi 2008; Yoo & Seljak 2015), though this presents difficulties since the powerspectrum becomes ill-defined and the dimensionality increases. For surveys below ∼ ° , the pairwise approximations are valid (Samushiaet al. 2015), yet we should bear in mind their sub-optimality for future surveys such as SPHEREx (Doré et al. 2014). Even with improvedestimators, it is necessary to model the wide-angle effects arising from other sources such as the galaxy selection function, and much workhas been done to achieve this goal. Indeed, it may prove more straightforward to forward-model also the post-Yamamoto effects, and use asimpler estimator; given the results of this work, there is little justification for this on efficiency grounds. On a more philosophical note, it isinteresting to see how even on the largest scales, where the Universe is most linear, the two-point function is still highly non-trivial. Despitebeing the most basic cosmological observable, the two-point function still has secrets up its sleeve. ACKNOWLEDGEMENTS
We thank Emanuele Castorina, Karolina Garcia, and David Spergel for insightful discussions. OHEP acknowledges funding from the WFIRSTprogram through NNG26PJ30C and NNN12AA01C. The authors are pleased to acknowledge that the work reported on in this paper wassubstantially performed using the Princeton Research Computing resources at Princeton University which is consortium of groups led by thePrinceton Institute for Computational Science and Engineering (PICSciE) and Office of Information Technology’s Research Computing.
APPENDIX A: USEFUL MATHEMATICAL RELATIONS
We present a number of mathematical results used in this work. Firstly, we give the explicit form for 𝐿 ℓ ( 𝜇 ) as a finite polynomial in 𝜇 (with | 𝜇 | ≤ 𝐿 ℓ ( 𝜇 ) = (cid:98) ℓ / (cid:99) ∑︁ 𝑛 = 𝑐 𝑛ℓ × 𝜇 ℓ − 𝑛 , 𝑐 𝑛ℓ = (− ) 𝑛 ℓ (cid:18) ℓ𝑛 (cid:19) (cid:18) ℓ − 𝑛ℓ (cid:19) ≡ (− ) 𝑛 ( ℓ − 𝑛 ) !2 ℓ 𝑛 ! ( ℓ − 𝑛 ) ! ( ℓ − 𝑛 ) ! (A1)(Abramowitz & Stegun 1964, Eq. 22.3.8), where the 2 × (cid:98) ℓ / (cid:99) indicates the largest integer ≤ ℓ /
2. Theinverse relation also proves useful: 𝜇 𝑛 = ∑︁ ℓ = 𝑛 ↓ ˜ 𝑐 ℓ𝑛 × 𝐿 ℓ ( 𝜇 ) , ˜ 𝑐 ℓ𝑛 = ( ℓ + ) 𝑛 + √ 𝜋 Γ ( + 𝑛 ) Γ (( 𝑛 − ℓ )/ + ) Γ (( 𝑛 + ℓ + )/ ) ≡ ( ℓ + ) 𝑛 ! ( ℓ + 𝑛 + ) !! ( ℓ + 𝑛 + ) ! ( 𝑛 − ℓ ) !! (A2)(derived from Gradshteyn & Ryzhik 1994, Eq. 7.126.1 invoking Legendre polynomial orthogonality), where the summation is over all ℓ downwards from 𝑛 in steps of two, and we use the double factorial 𝑛 !! ≡ 𝑛 ( 𝑛 − )( 𝑛 − ) ... for positive integer 𝑛 . Note that Legendrepolynomials of even (odd) order depend only on even (odd) powers of 𝜇 (and vice versa ). To derive the simplified coefficients on theright-hand-side we have assumed 𝑛 to be integral, and noted that the summation rules imply that ℓ , 𝑛 have the same sign. Inserting (A2) into(A1) and using orthogonality gives the relation 𝐿 ℓ ( 𝜇 ) = (cid:98) ℓ / (cid:99) ∑︁ 𝑛 = ∑︁ 𝐿 = ( ℓ − 𝑛 )↓ 𝑐 𝑛ℓ ˜ 𝑐 𝐿ℓ − 𝑛 𝐿 𝐿 ( 𝜇 ) ⇒ (cid:98) ℓ / (cid:99) ∑︁ 𝑛 = 𝑐 𝑛ℓ ˜ 𝑐 𝐿ℓ − 𝑛 = 𝛿 K ℓ𝐿 , (A3)where 𝛿 K is the Kronecker delta.Legendre polynomials have the generating function (cid:104) − 𝑠𝑡 + 𝑡 (cid:105) − / = ∞ ∑︁ 𝑁 = 𝐿 𝐿 ( 𝑠 ) 𝑡 𝑁 (A4)(Gradshteyn & Ryzhik 1994, Eq. 8.921), and may related to spherical harmonics via the addition theorem : 𝐿 ℓ ( ˆ x · ˆ y ) = 𝜋 ℓ + ℓ ∑︁ 𝑚 = − ℓ 𝑌 𝑚ℓ ( ˆ x ) 𝑌 𝑚 ∗ ℓ ( ˆ y ) (A5) MNRAS , 1–27 (2021) eyond Yamamoto Correlators (NIST DLMF, Eq. 14.30.9), where ∗ indicates a complex conjugate and 𝑌 𝑚ℓ is the spherical harmonic of order ( ℓ, 𝑚 ) , which obeys thesymmetries 𝑌 𝑚ℓ (− ˆ x ) = (− ) ℓ 𝑌 𝑚ℓ ( ˆ x ) , 𝑌 𝑚 ∗ ℓ ( ˆ x ) = (− ) 𝑚 𝑌 − 𝑚ℓ ( ˆ x ) and 𝑌 = ( 𝜋 ) − / . Furthermore, the spherical harmonics are orthonormal: ∫ 𝑑 Ω 𝑥 𝑌 𝑚ℓ ( ˆ x ) 𝑌 𝑚 (cid:48) ∗ ℓ (cid:48) ( ˆ x ) = 𝛿 K ℓℓ (cid:48) 𝛿 K 𝑚𝑚 (cid:48) (A6)(NIST DLMF, Eq. 14.30.8). An important relation is the product-to-sum rule for spherical harmonics: 𝑌 𝑀 𝐿 ( ˆ x ) 𝑌 𝑀 𝐿 ( ˆ x ) = 𝐿 + 𝐿 ∑︁ 𝐿 = | 𝐿 − 𝐿 | 𝐿 ∑︁ 𝑀 = − 𝐿 G 𝑀 𝑀 (− 𝑀 ) 𝐿 𝐿 𝐿 (− ) 𝑀 𝑌 𝑀 𝐿 ( ˆ x ) , (A7)where G is the Gaunt integral obeying several selection rules including 𝑀 + 𝑀 + 𝑀 =
0. This can be derived from spherical harmonicorthogonality (A6) and the definition of G as the integral over three spherical harmonics (NIST DLMF, Eq. 34.3.22). The Gaunt factor maybe written explicitly in terms of Wigner 3- 𝑗 symbols as G ℓ ℓ ℓ 𝑚 𝑚 𝑚 = (cid:18) ( ℓ + )( ℓ + )( ℓ + ) 𝜋 (cid:19) / (cid:18) ℓ ℓ ℓ 𝑚 𝑚 𝑚 (cid:19) (cid:18) ℓ ℓ ℓ (cid:19) , (A8)(NIST DLMF, Eq. 34.3.22).Closely related to spherical harmonics are the regular solid harmonics , defined by 𝑅 𝑚ℓ ( r ) = √︂ 𝜋 ℓ + 𝑟 ℓ 𝑌 𝑚ℓ ( ˆ r ) . (A9)These obey an addition theorem: 𝑅 𝑚ℓ ( r + a ) = ℓ ∑︁ 𝜆 = 𝜆 ∑︁ 𝜇 = − 𝜆 (cid:18) ℓ + 𝑚𝜆 + 𝜇 (cid:19) / (cid:18) ℓ − 𝑚𝜆 − 𝜇 (cid:19) / 𝑅 𝜇𝜆 ( r ) 𝑅 𝑚 − 𝜇ℓ − 𝜆 ( a ) (A10)(Tough & Stone 1977), which is a finite series.The Rayleigh plane wave expansion gives 𝑒 𝑖 x · y = ∞ ∑︁ ℓ = 𝑖 ℓ ( ℓ + ) 𝑗 ℓ ( 𝑥𝑦 ) 𝐿 ℓ ( ˆ x · ˆ y ) = 𝜋 ∞ ∑︁ ℓ = ℓ ∑︁ 𝑚 = − ℓ 𝑖 ℓ 𝑗 ℓ ( 𝑥𝑦 ) 𝑌 𝑚ℓ ( ˆ x ) 𝑌 𝑚 ∗ ℓ ( ˆ y ) (A11)(Arfken et al. 2013, Eq. 16.63), where we have used (A5) to arrive at the second equality. Using this, we can prove the following identity: ∫ 𝑑 Ω 𝑥 𝑒 𝑖 x · y 𝐿 ℓ ( x · z ) = 𝜋 ∑︁ ℓ (cid:48) = ℓ (cid:48) ∑︁ 𝑚 (cid:48) = − ℓ (cid:48) 𝑖 ℓ (cid:48) 𝑗 ℓ (cid:48) ( 𝑥𝑦 ) ∫ 𝑑 Ω 𝑥 𝑌 𝑚 (cid:48) ℓ (cid:48) ( ˆ x ) 𝑌 𝑚 (cid:48) ∗ ℓ (cid:48) ( ˆ y ) 𝜋 ℓ + ℓ ∑︁ 𝑚 = − ℓ 𝑌 𝑚ℓ ( ˆ z ) 𝑌 𝑚 ∗ ℓ ( ˆ x ) (A12) = 𝜋 𝑖 ℓ 𝑗 ℓ ( 𝑥𝑦 ) 𝜋 ℓ + ℓ ∑︁ 𝑚 = − ℓ 𝑌 𝑚 ∗ ℓ ( ˆ y ) 𝑌 𝑚 ∗ ℓ ( ˆ z ) ≡ 𝜋 𝑖 ℓ 𝑗 ℓ ( 𝑥𝑦 ) 𝐿 ℓ ( ˆ y · ˆ z ) , using (A5) & (A11) in the first line, and (A6) to obtain the second.An additional class of orthogonal polynomials are Gegenbauer (or ultra-spherical) polynomials, with the generating function (cid:104) − 𝑠𝑡 + 𝑡 (cid:105) − 𝛼 = ∞ ∑︁ 𝑁 = 𝐶 ( 𝛼 ) 𝑁 ( 𝑠 ) 𝑡 𝑁 (A13)(Gradshteyn & Ryzhik 1994, Eq. 8.930), where 𝐶 ( 𝛼 ) 𝑁 is the Gegenbauer polynomial (for integer order 𝑁 ), and 𝛼 = / 𝐶 ( 𝛼 ) 𝑁 ( 𝑠 ) = (cid:98) 𝑁 / (cid:99) ∑︁ 𝑘 = 𝐶 𝛼𝑁 𝑘 × ( 𝑠 ) 𝑁 − 𝑘 , 𝐶 𝛼𝑁 𝑘 = (− ) 𝑘 Γ ( 𝛼 + 𝑁 − 𝑘 ) Γ ( 𝛼 ) 𝑘 ! ( 𝑁 − 𝑘 ) ! (A14)(Abramowitz & Stegun 1964, Eq. 22.3.4), with the special case of 𝐶 𝑁 𝑘 = 𝛿 K 𝑁 . Note that this is a sum of even (odd) powers of 𝑠 for even(odd) 𝑁 , as for the Legendre polynomials. MNRAS000
0. This can be derived from spherical harmonicorthogonality (A6) and the definition of G as the integral over three spherical harmonics (NIST DLMF, Eq. 34.3.22). The Gaunt factor maybe written explicitly in terms of Wigner 3- 𝑗 symbols as G ℓ ℓ ℓ 𝑚 𝑚 𝑚 = (cid:18) ( ℓ + )( ℓ + )( ℓ + ) 𝜋 (cid:19) / (cid:18) ℓ ℓ ℓ 𝑚 𝑚 𝑚 (cid:19) (cid:18) ℓ ℓ ℓ (cid:19) , (A8)(NIST DLMF, Eq. 34.3.22).Closely related to spherical harmonics are the regular solid harmonics , defined by 𝑅 𝑚ℓ ( r ) = √︂ 𝜋 ℓ + 𝑟 ℓ 𝑌 𝑚ℓ ( ˆ r ) . (A9)These obey an addition theorem: 𝑅 𝑚ℓ ( r + a ) = ℓ ∑︁ 𝜆 = 𝜆 ∑︁ 𝜇 = − 𝜆 (cid:18) ℓ + 𝑚𝜆 + 𝜇 (cid:19) / (cid:18) ℓ − 𝑚𝜆 − 𝜇 (cid:19) / 𝑅 𝜇𝜆 ( r ) 𝑅 𝑚 − 𝜇ℓ − 𝜆 ( a ) (A10)(Tough & Stone 1977), which is a finite series.The Rayleigh plane wave expansion gives 𝑒 𝑖 x · y = ∞ ∑︁ ℓ = 𝑖 ℓ ( ℓ + ) 𝑗 ℓ ( 𝑥𝑦 ) 𝐿 ℓ ( ˆ x · ˆ y ) = 𝜋 ∞ ∑︁ ℓ = ℓ ∑︁ 𝑚 = − ℓ 𝑖 ℓ 𝑗 ℓ ( 𝑥𝑦 ) 𝑌 𝑚ℓ ( ˆ x ) 𝑌 𝑚 ∗ ℓ ( ˆ y ) (A11)(Arfken et al. 2013, Eq. 16.63), where we have used (A5) to arrive at the second equality. Using this, we can prove the following identity: ∫ 𝑑 Ω 𝑥 𝑒 𝑖 x · y 𝐿 ℓ ( x · z ) = 𝜋 ∑︁ ℓ (cid:48) = ℓ (cid:48) ∑︁ 𝑚 (cid:48) = − ℓ (cid:48) 𝑖 ℓ (cid:48) 𝑗 ℓ (cid:48) ( 𝑥𝑦 ) ∫ 𝑑 Ω 𝑥 𝑌 𝑚 (cid:48) ℓ (cid:48) ( ˆ x ) 𝑌 𝑚 (cid:48) ∗ ℓ (cid:48) ( ˆ y ) 𝜋 ℓ + ℓ ∑︁ 𝑚 = − ℓ 𝑌 𝑚ℓ ( ˆ z ) 𝑌 𝑚 ∗ ℓ ( ˆ x ) (A12) = 𝜋 𝑖 ℓ 𝑗 ℓ ( 𝑥𝑦 ) 𝜋 ℓ + ℓ ∑︁ 𝑚 = − ℓ 𝑌 𝑚 ∗ ℓ ( ˆ y ) 𝑌 𝑚 ∗ ℓ ( ˆ z ) ≡ 𝜋 𝑖 ℓ 𝑗 ℓ ( 𝑥𝑦 ) 𝐿 ℓ ( ˆ y · ˆ z ) , using (A5) & (A11) in the first line, and (A6) to obtain the second.An additional class of orthogonal polynomials are Gegenbauer (or ultra-spherical) polynomials, with the generating function (cid:104) − 𝑠𝑡 + 𝑡 (cid:105) − 𝛼 = ∞ ∑︁ 𝑁 = 𝐶 ( 𝛼 ) 𝑁 ( 𝑠 ) 𝑡 𝑁 (A13)(Gradshteyn & Ryzhik 1994, Eq. 8.930), where 𝐶 ( 𝛼 ) 𝑁 is the Gegenbauer polynomial (for integer order 𝑁 ), and 𝛼 = / 𝐶 ( 𝛼 ) 𝑁 ( 𝑠 ) = (cid:98) 𝑁 / (cid:99) ∑︁ 𝑘 = 𝐶 𝛼𝑁 𝑘 × ( 𝑠 ) 𝑁 − 𝑘 , 𝐶 𝛼𝑁 𝑘 = (− ) 𝑘 Γ ( 𝛼 + 𝑁 − 𝑘 ) Γ ( 𝛼 ) 𝑘 ! ( 𝑁 − 𝑘 ) ! (A14)(Abramowitz & Stegun 1964, Eq. 22.3.4), with the special case of 𝐶 𝑁 𝑘 = 𝛿 K 𝑁 . Note that this is a sum of even (odd) powers of 𝑠 for even(odd) 𝑁 , as for the Legendre polynomials. MNRAS000 , 1–27 (2021) Philcox & Slepian
APPENDIX B: SPHERICAL HARMONIC SHIFT THEOREMB1 Derivation
We present the derivation of a useful series expansion for spherical harmonics. First, we consider the function 𝑌 𝑚ℓ ( (cid:154) a + A ) for 𝜖 ≡ 𝑎 / 𝐴 (cid:28) 𝑅 𝑚ℓ ( a + A ) then using the addition theorem (A10), we obtain 𝑌 𝑚ℓ ( (cid:154) a + A ) = ℓ ∑︁ 𝜆 = 𝜆 ∑︁ 𝜇 = − 𝜆 √︄ 𝜋 ( ℓ + )( 𝜆 + )( ℓ − 𝜆 + ) (cid:18) ℓ + 𝑚𝜆 + 𝜇 (cid:19) / (cid:18) ℓ − 𝑚𝜆 − 𝜇 (cid:19) / 𝑎 𝜆 𝐴 ℓ − 𝜆 | a + A | ℓ 𝑌 𝜇𝜆 ( ˆ a ) 𝑌 𝑚 − 𝜇ℓ − 𝜆 ( ˆ A ) (B1) ≡ ℓ ∑︁ 𝜆 = 𝜆 ∑︁ 𝜇 = − 𝜆 𝐴 𝜆𝜇ℓ𝑚 𝜖 𝜆 | ˆ A + 𝜖 ˆ a | ℓ 𝑌 𝜇𝜆 ( ˆ a ) 𝑌 𝑚 − 𝜇ℓ − 𝜆 ( ˆ A ) , defining the coefficients 𝐴 𝜆𝜇ℓ𝑚 in the second line. This is a finite expansion, allowing the composite spherical harmonic to be expressed interms of the harmonics of ˆ a and ˆ A . However, due to the angular dependence of | ˆ A + 𝜖 ˆ a | ℓ , we require an infinite expansion to separate a and A fully. To proceed, we recognize that the denominator is the generating function for a Gegenbauer series (A13) with 𝛼 = ℓ / 𝑡 = 𝜖 , 𝑠 = − ˆ a · ˆ A :1 | ˆ A + 𝜖 ˆ a | ℓ = (cid:16) + 𝜖 ˆ a · ˆ A + 𝜖 (cid:17) ℓ / = ∞ ∑︁ 𝑁 = 𝐶 ( ℓ / ) 𝑁 (− ˆ a · ˆ A ) 𝜖 𝑁 = ∞ ∑︁ 𝑁 = (cid:98) 𝑁 / (cid:99) ∑︁ 𝑘 = 𝐶 ℓ / 𝑁 𝑘 × (− a · ˆ A ) 𝑁 − 𝑘 𝜖 𝑁 , (B2)where the coefficients 𝐶 ℓ / 𝑁 𝑘 are defined in (A14). Inserting into (B1) gives: 𝑌 𝑚ℓ ( (cid:154) a + A ) = ℓ ∑︁ 𝜆 = 𝜆 ∑︁ 𝜇 = − 𝜆 𝐴 𝜆𝜇ℓ𝑚 ∞ ∑︁ 𝑁 = (cid:98) 𝑁 / (cid:99) ∑︁ 𝑘 = 𝐶 ℓ / 𝑁 𝑘 (− ) 𝑁 − 𝑘 × 𝜖 𝜆 + 𝑁 𝑌 𝜇𝜆 ( ˆ a ) 𝑌 𝑚 − 𝜇ℓ − 𝜆 ( ˆ A )( ˆ a · ˆ A ) 𝑁 − 𝑘 . (B3)To simplify this further, we express ( ˆ a · ˆ A ) 𝑁 − 𝑘 in spherical harmonics, using (A2): ( ˆ a · ˆ A ) 𝑁 − 𝑘 = ∑︁ 𝐿 (cid:48) = ( 𝑁 − 𝑘 )↓ ˜ 𝑐 𝐿 (cid:48) 𝑁 − 𝑘 𝐿 𝐿 (cid:48) ( ˆ a · ˆ A ) = ∑︁ 𝐿 (cid:48) = ( 𝑁 − 𝑘 )↓ 𝜋 𝐿 (cid:48) + 𝐿 (cid:48) ∑︁ 𝑀 (cid:48) = − 𝐿 (cid:48) ˜ 𝑐 𝐿 (cid:48) 𝑁 − 𝑘 𝑌 𝑀 (cid:48) 𝐿 (cid:48) ( ˆ a ) 𝑌 𝑀 (cid:48) ∗ 𝐿 (cid:48) ( ˆ A ) , (B4)expanding the Legendre polynomial into spherical harmonics via the addition theorem (A5). Finally, we can combine the two sphericalharmonics in ˆ a and ˆ A using the product-to-sum relation (A7): 𝑌 𝜇𝜆 ( ˆ a ) 𝑌 𝑀 (cid:48) 𝐿 (cid:48) ( ˆ a ) = 𝜆 + 𝐿 (cid:48) ∑︁ 𝐽 = | 𝜆 − 𝐿 (cid:48) | 𝐽 ∑︁ 𝑀 = − 𝐽 G 𝜇𝑀 (cid:48) (− 𝑀 ) 𝜆𝐿 (cid:48) 𝐽 (− ) 𝑀 𝑌 𝑀 𝐽 ( ˆ a ) (B5) 𝑌 𝑚 − 𝜇ℓ − 𝜆 ( ˆ A ) 𝑌 𝑀 (cid:48) ∗ 𝐿 (cid:48) ( ˆ A ) = ℓ − 𝜆 + 𝐿 (cid:48) ∑︁ 𝐽 = | ℓ − 𝜆 − 𝐿 (cid:48) | 𝐽 ∑︁ 𝑀 = − 𝐽 G ( 𝑚 − 𝜇 ) (− 𝑀 (cid:48) ) (− 𝑀 )( ℓ − 𝜆 ) 𝐿 (cid:48) 𝐽 (− ) 𝑀 (cid:48) + 𝑀 𝑌 𝑀 𝐽 ( ˆ A ) . Combining results, we obtain the shift theorem for spherical harmonics: 𝑌 𝑚ℓ ( (cid:154) a + A ) = ℓ ∑︁ 𝜆 = 𝜆 ∑︁ 𝜇 = − 𝜆 𝐴 𝜆𝜇ℓ𝑚 ∞ ∑︁ 𝑁 = (cid:98) 𝑁 / (cid:99) ∑︁ 𝑘 = 𝐶 ℓ / 𝑁 𝑘 (− ) 𝑁 − 𝑘 ∑︁ 𝐿 (cid:48) = ( 𝑁 − 𝑘 )↓ ˜ 𝑐 𝐿 (cid:48) 𝑁 − 𝑘 𝜋 𝐿 (cid:48) + × 𝜆 + 𝐿 (cid:48) ∑︁ 𝐽 = | 𝜆 − 𝐿 (cid:48) | 𝐽 ∑︁ 𝑀 = − 𝐽 G 𝜇 ( 𝑀 − 𝜇 ) (− 𝑀 ) 𝜆𝐿 (cid:48) 𝐽 ℓ − 𝜆 + 𝐿 (cid:48) ∑︁ 𝐽 = | ℓ − 𝜆 − 𝐿 (cid:48) | G ( 𝑚 − 𝜇 ) ( 𝜇 − 𝑀 ) ( 𝑀 − 𝑚 )( ℓ − 𝜆 ) 𝐿 (cid:48) 𝐽 (− ) 𝑚 + 𝑀 − 𝜇 × 𝜖 𝜆 + 𝑁 𝑌 𝑀𝐽 ( ˆ a ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ A ) , where we have noted that the Gaunt integrals imply 𝑀 = 𝜇 + 𝑀 (cid:48) , 𝑀 = 𝑚 − 𝜇 − 𝑀 (cid:48) = 𝑚 − 𝑀 and relabelled 𝑀 (cid:48) → 𝑀 . Defining thecoefficients 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 where 𝛼 = 𝜆 + 𝑁 , this can be written more succinctly as 𝑌 𝑚ℓ ( (cid:154) a + A ) = ∞ ∑︁ 𝛼 = 𝛼 ∑︁ 𝐽 = ℓ + 𝛼 ∑︁ 𝐽 = max ( ,ℓ − 𝛼 ) 𝐽 ∑︁ 𝑀 = − 𝐽 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 (cid:16) 𝑎𝐴 (cid:17) 𝛼 𝑌 𝑀𝐽 ( ˆ a ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ A ) , (B7)where the summation limits will be elaborated upon in the next section. Truncating (B7) at 𝛼 = 𝐾 incurs an error of O( 𝜖 𝐾 + ) , thus, for 𝜖 < This is equivalent to that introduced in Garcia & Slepian (2020) for the three-point correlation function except with one position vector set to zero, affordingsignificant simplification. MNRAS , 1–27 (2021) eyond Yamamoto Correlators B2 Coefficient Properties
From (B6), the shift coefficients 𝜑 are given by 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 = ℓ ∑︁ 𝜆 = 𝜆 ∑︁ 𝜇 = − 𝜆 𝐴 𝜆𝜇ℓ𝑚 ∞ ∑︁ 𝑁 = 𝛿 K 𝛼 ( 𝜆 + 𝑁 ) (cid:98) 𝑁 / (cid:99) ∑︁ 𝑘 = 𝐶 ℓ / 𝑁 𝑘 (− ) 𝑁 − 𝑘 ∑︁ 𝐿 (cid:48) = ( 𝑁 − 𝑘 )↓ ˜ 𝑐 𝐿 (cid:48) 𝑁 − 𝑘 𝜋 𝐿 (cid:48) + × G 𝜇 ( 𝑀 − 𝜇 ) (− 𝑀 ) 𝜆𝐿 (cid:48) 𝐽 G ( 𝑚 − 𝜇 ) ( 𝜇 − 𝑀 ) ( 𝑀 − 𝑚 )( ℓ − 𝜆 ) 𝐿 (cid:48) 𝐽 (− ) 𝑚 + 𝑀 − 𝜇 , (B8)where the Kronecker delta in the first line enforces that 𝜆 + 𝑁 = 𝛼 . Including the various coefficient definitions gives the explicit form 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 = √︁ 𝜋 ( ℓ + )( 𝐽 + )( 𝐽 + ) min ( 𝛼,ℓ ) ∑︁ 𝜆 = 𝜆 ∑︁ 𝜇 = − 𝜆 (cid:18) ℓ + 𝑚𝜆 + 𝜇 (cid:19) / (cid:18) ℓ − 𝑚𝜆 − 𝜇 (cid:19) / (cid:98) 𝛼 − 𝜆 / (cid:99) ∑︁ 𝑘 = ∑︁ 𝐿 (cid:48) = ( 𝛼 − 𝜆 − 𝑘 )↓ (B9) × (− ) 𝑘 + 𝛼 + 𝜆 + 𝑚 + 𝑀 + 𝜇 ( 𝐿 (cid:48) + )√ 𝜋 Γ (cid:16) ℓ + 𝛼 − 𝜆 − 𝑘 (cid:17) Γ ( ℓ / ) Γ (( 𝛼 − 𝜆 − 𝑘 − 𝐿 (cid:48) )/ + ) Γ (( 𝛼 − 𝜆 − 𝑘 + 𝐿 (cid:48) + )/ ) 𝑘 ! × (cid:18) 𝜆 𝐿 (cid:48) 𝐽 𝜇 𝑀 − 𝜇 − 𝑀 (cid:19) (cid:18) ℓ − 𝜆 𝐿 (cid:48) 𝐽 𝑚 − 𝜇 𝜇 − 𝑀 𝑀 − 𝑚 (cid:19) (cid:18) 𝜆 𝐿 (cid:48) 𝐽 (cid:19) (cid:18) ℓ − 𝜆 𝐿 (cid:48) 𝐽 (cid:19) . We consider two special cases. For 𝛼 =
0, we require 𝜆 = 𝑁 =
0, and thus, by the summation limits, 𝜇 = 𝑘 = 𝐿 (cid:48) =
0. This gives 𝜑 ,ℓ𝑚𝐽 𝐽 𝑀 = ( 𝜋 ) / G 𝑀 (− 𝑀 ) 𝐽 G 𝑚 (− 𝑀 ) 𝑀ℓ 𝐽 (− ) 𝑚 + 𝑀 , (B10)using 𝐴 ℓ𝑚 = √ 𝜋 , ˜ 𝑐 = 𝐶 ℓ / = 𝑗 symmetries, such that G 𝑀 𝑀 𝑀 𝐿 𝐿 𝐿 is zero unless | 𝐿 − 𝐿 | <𝐿 < 𝐿 + 𝐿 , implying 𝐽 = 𝑀 = 𝐽 = ℓ . In full, we obtain 𝜑 ,ℓ𝑚𝐽 𝐽 𝑀 = 𝛿 K 𝐽 𝛿 K 𝐽 ℓ 𝛿 K 𝑀 × √ 𝜋. (B11)This gives 𝑌 𝑚ℓ ( (cid:154) a + A ) (cid:12)(cid:12)(cid:12) 𝛼 = = 𝑌 𝑚ℓ ( ˆ A ) , as expected.Secondly, we consider ℓ = 𝑚 =
0. The summation limits enforce 𝜆 = 𝜇 =
0, and, since 𝐶 𝑁 𝑘 = 𝛿 K 𝑁 (A14), 𝑁 = 𝑘 = 𝐿 (cid:48) = 𝛼 = 𝜆 + 𝑁 = 𝜑 𝛼, 𝐽 𝐽 𝑀 = ( 𝜋 ) / G 𝑀 (− 𝑀 ) 𝐽 G (− 𝑀 ) 𝑀 𝐽 (− ) 𝑀 . (B12)In this case, the triangle conditions imply 𝐽 = 𝐽 =
0, thus 𝑀 =
0, and 𝜑 𝛼, 𝐽 𝐽 𝑀 = 𝛿 K 𝐽 𝛿 K 𝐽 𝛿 K 𝑀 𝛿 K 𝛼 × √ 𝜋, (B13)giving 𝑌 ( (cid:154) a + A ) = ( 𝜋 ) − / , as expected.It is useful to consider properties relating to parity. Firstly, consider the parity transformation a → − a , A → − A : 𝑌 𝑚ℓ ( (cid:154) a + A ) → (− ) ℓ 𝑌 𝑚ℓ ( (cid:154) a + A ) (B14) 𝑌 𝑀𝐽 ( ˆ a ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ A ) → (− ) 𝐽 + 𝐽 𝑌 𝑀𝐽 ( ˆ a ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ A ) . By orthogonality of the spherical harmonics, this implies 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 = (− ) ℓ + 𝐽 + 𝐽 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 , and hence that the coefficients are vanishing if ℓ + 𝐽 + 𝐽 is odd. Furthermore, the Gaunt coefficients G 𝑀 𝑀 𝑀 𝐿 𝐿 𝐿 contain 3- 𝑗 symbols with all 𝑀 𝑖 fixed to zero; these vanish if 𝐿 + 𝐿 + 𝐿 is odd, implying that 𝜆 + 𝐿 (cid:48) + 𝐽 and ℓ − 𝜆 + 𝐿 (cid:48) + 𝐽 must be even (B8). Finally, the summation on 𝐿 (cid:48) indicates that 𝐿 (cid:48) and 𝑁 must have thesame sign, hence ℓ + 𝑁 + 𝐽 is even, thus so is 𝛼 + 𝐽 . Our conclusion is that 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 must be vanishing unless both 𝛼 + 𝐽 and ℓ + 𝛼 + 𝐽 areeven. This result proves useful when forming the even parity expansion in §3.3.Finally, we discuss the summation limits in (B7). From the 𝐿 (cid:48) summation, 0 ≤ 𝐿 (cid:48) ≤ 𝑁 (since 0 ≤ 𝑘 ≤ 𝑁 ), and from the Gauntintegrals, | 𝜆 − 𝐿 (cid:48) | ≤ 𝐽 ≤ 𝜆 + 𝐿 (cid:48) , | ℓ − 𝜆 − 𝐿 (cid:48) | ≤ 𝐽 ≤ ℓ − 𝜆 + 𝐿 (cid:48) . Combining the two implies that 0 ≤ 𝐽 ≤ 𝛼 , max ( , ℓ − 𝛼 ) ≤ 𝐽 ≤ ℓ + 𝛼 , | 𝑀 | ≤ 𝐽 , noting that the moduli enforce 𝐽 , 𝐽 to be positive. For the 𝛼 = 𝐽 = 𝛼 + 𝐽 is even)and 𝐽 = ℓ ±
1, or, for 𝛼 = 𝐽 ∈ { , } , 𝐽 ∈ { ℓ, ℓ ± } . APPENDIX C: GENERALIZED SPHERICAL HARMONIC SHIFT THEOREM
Below, we prove a useful generalization of the spherical harmonic shift theorem discussed in Appendix B. To begin, consider the expression (cid:18) | a || A + a | (cid:19) 𝛼 𝑌 𝑀𝐽 ( ˆ a ) 𝑌 𝑚 − 𝑀𝐽 ( (cid:154) a + A ) = (cid:2) + 𝜖 ˆ a · ˆ A + 𝜖 (cid:3) 𝛼 / 𝑌 𝑀𝐽 ( ˆ a ) 𝑌 𝑚 − 𝑀𝐽 ( (cid:154) a + A ) , (C1)for some ˆ a , ˆ A with 𝜖 = 𝑎 / 𝐴 (cid:28)
1. The spherical harmonic in a + A can be expanded identically to Appendix B1; with the prefactor (whichis again a generator for Gegenbauer polynomials), we obtain the same expression except that 𝐶 ℓ / 𝑁 𝑘 is replaced with 𝐶 ( ℓ + 𝛼 )/ 𝑁 𝑘 in (B2) and the
MNRAS000
MNRAS000 , 1–27 (2021) Philcox & Slepian succeeding. In full, we obtain (cid:18) | a || A + a | (cid:19) 𝛼 𝑌 𝑀𝐽 ( ˆ a ) 𝑌 𝑚 − 𝑀𝐽 ( (cid:154) a + A ) = ∑︁ 𝛽𝐽 (cid:48) 𝐽 (cid:48) 𝑀 (cid:48) 𝜔 𝛽,𝐽 ( 𝑚 − 𝑀 ) 𝛼𝐽 (cid:48) 𝐽 (cid:48) 𝑀 (cid:48) 𝜖 𝛼 + 𝛽 𝑌 𝑀𝐽 ( ˆ a ) 𝑌 𝑀 (cid:48) 𝐽 (cid:48) ( ˆ a ) 𝑌 𝑚 − 𝑀 − 𝑀 (cid:48) 𝐽 (cid:48) ( ˆ A ) , (C2)where 𝜔 𝛽,ℓ𝑚𝛼𝐽 𝐽 𝑀 is equal to 𝜑 𝛽,ℓ𝑚𝐽 𝐽 𝑀 (B8), except that 𝐶 ℓ / 𝑁 𝑘 → 𝐶 ( ℓ + 𝛼 )/ 𝑁 𝑘 , implying that 𝜔 𝛽,ℓ𝑚 𝐽 𝐽 𝑀 = 𝜑 𝛽,ℓ𝑚𝐽 𝐽 𝑀 . Next, we can combine the twospherical harmonics in a using the product-to-sum relation (A7) to yield (cid:18) | a || A + a | (cid:19) 𝛼 𝑌 𝑀𝐽 ( ˆ a ) 𝑌 𝑚 − 𝑀𝐽 ( (cid:154) a + A ) = ∑︁ 𝛽𝑆 𝑆 𝑇 (cid:32)∑︁ 𝐽 (cid:48) 𝜔 𝛽,𝐽 ( 𝑚 − 𝑀 ) 𝛼𝐽 (cid:48) 𝑆 ( 𝑇 − 𝑀 ) G 𝑀 ( 𝑇 − 𝑀 ) (− 𝑇 ) 𝐽 𝐽 (cid:48) 𝑆 (− ) 𝑇 (cid:33) 𝜖 𝛼 + 𝛽 𝑌 𝑇𝑆 ( ˆ a ) 𝑌 𝑚 − 𝑇𝑆 ( ˆ A )⇒ (cid:18) | a || A + a | (cid:19) 𝛼 𝑌 𝑀𝐽 ( ˆ a ) 𝑌 𝑚 − 𝑀𝐽 ( (cid:154) a + A ) ≡ ∑︁ 𝛽𝑆 𝑆 𝑇 Ω 𝛽,𝛼𝐽 𝐽 𝑚𝑀𝑆 𝑆 𝑇 𝜖 𝛼 + 𝛽 𝑌 𝑇𝑆 ( ˆ a ) 𝑌 𝑚 − 𝑇𝑆 ( ˆ A ) , (C3)liberally relabeling variables and defining the new coefficient set Ω . The expression is now in the form of a power series combined with twospherical harmonics, just as in (B7).By analogous arguments to before, the 𝜔 𝛽,𝐽 ( 𝑚 − 𝑀 ) 𝛼𝐽 (cid:48) 𝑆 ( 𝑇 − 𝑀 ) coefficient requires 0 ≤ 𝐽 (cid:48) ≤ 𝛽 and max ( , 𝐽 − 𝛽 ) ≤ 𝑆 ≤ 𝐽 + 𝛽 . Coupled withthe triangle conditions on the additional Gaunt integral, we obtain max ( , 𝐽 − 𝛽 ) ≤ 𝑆 ≤ 𝐽 + 𝛽 , max ( , 𝐽 − 𝛽 ) ≤ 𝑆 ≤ 𝐽 + 𝛽 , and | 𝑇 | < 𝑆 .Furthermore, we require 𝛽 + 𝐽 (cid:48) and 𝐽 + 𝛽 + 𝑆 to be even, and, from the Gaunt factor, the same applies for 𝛽 + 𝐽 + 𝑆 . The case 𝛽 = 𝜔 ,𝐽 ( 𝑚 − 𝑀 ) 𝛼𝐽 (cid:48) 𝑆 ( 𝑇 − 𝑀 ) = √ 𝜋𝛿 K 𝐽 (cid:48) 𝛿 K 𝑆 𝐽 𝛿 K 𝑇 𝑀 , thus Ω ,𝛼𝐽 𝐽 𝑚𝑀𝑆 𝑆 𝑇 = √ 𝜋𝛿 K 𝑆 𝐽 𝛿 K 𝑇 𝑀 G 𝑀 (− 𝑀 ) 𝐽 𝑆 (− ) 𝑀 . The Gaunt integralrequires 𝑆 = 𝐽 , thus we obtain Ω ,𝛼𝐽 𝐽 𝑚𝑀𝑆 𝑆 𝑇 = 𝛿 K 𝑆 𝐽 𝛿 K 𝑆 𝐽 𝛿 K 𝑇 𝑀 . (C4) APPENDIX D: ALTERNATIVE DERIVATION FOR THE PARITY EVEN EXPANSION
We briefly present an alternative derivation of the parity-even expansion for spherical harmonics, which, whilst less conceptual, results in aneasier-to-implement formalism. To begin, consider the function 𝑧 𝛼,𝑚ℓ ≡ (cid:104) 𝜖 𝛼 𝑌 𝑚ℓ ( ˆ r ) + (− ) ℓ 𝜖 𝛼 𝑌 𝑚ℓ ( ˆ r ) (cid:105) . (D1)Writing r = R − 𝚫 , r = R + 𝚫 , each term may be expanded using a part of the generalized shift theorem of Appendix C: 𝑧 𝛼,𝑚ℓ ≡ (cid:20)(cid:18) Δ | R − 𝚫 | (cid:19) 𝛼 𝑌 𝑚ℓ ( (cid:155) R − 𝚫 ) + (− ) ℓ (cid:18) Δ | R + 𝚫 | (cid:19) 𝛼 𝑌 𝑚ℓ ( (cid:154) R + 𝚫 ) (cid:21) (D2) = ∑︁ 𝛽𝐽 𝐽 𝑀 𝜔 𝛽,ℓ𝑚𝛼𝐽 𝐽 𝑀 𝜃 𝛼 + 𝛽 𝑌 𝑀𝐽 ( ˆ 𝚫 ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ R ) (cid:20) + (− ) ℓ + 𝐽 (cid:21) . using (C2) and recalling 𝜃 ≡ Δ / 𝑅 . As before, 𝜔 𝛽,ℓ𝑚𝛼𝐽 𝐽 𝑀 is non-zero unless 𝐽 + 𝛽 is even, thus the sum extends only over even 𝛽 (assumingeven ℓ ). Extracting the 𝛽 = 𝑧 𝛼,𝑚ℓ ≡ 𝜃 𝛼 𝑌 𝑚ℓ ( ˆ R ) + ∑︁ even 𝛽> ∑︁ 𝐽 𝐽 𝑀 𝜔 𝛽,ℓ𝑚𝛼𝐽 𝐽 𝑀 𝑌 𝑀𝐽 ( ˆ 𝚫 ) × 𝜃 𝛼 + 𝛽 𝑌 𝑚 − 𝑀𝐽 ( ˆ R ) . (D3)Assuming small 𝜃 , this may be inverted order-by-order to find 𝜃 𝛼 𝑌 𝑚ℓ in terms of 𝑧 𝛼,𝑚ℓ (which is in the form required for the convolutionalpower spectrum estimators). This gives 𝜃 𝛼 𝑌 𝑚ℓ ( ˆ R ) = 𝑧 𝛼,𝑚ℓ − ∑︁ even 𝛽> ∑︁ 𝐽 𝐽 𝑀 𝜔 𝛽,ℓ𝑚𝛼𝐽 𝐽 𝑀 𝑌 𝑀𝐽 ( ˆ 𝚫 ) × 𝜃 𝛼 + 𝛽 𝑌 𝑚 − 𝑀𝐽 ( ˆ R ) (D4) = 𝑧 𝛼,𝑚ℓ − ∑︁ even 𝛽> ∑︁ 𝐽 𝐽 𝑀 𝜔 𝛽,ℓ𝑚𝛼𝐽 𝐽 𝑀 𝑌 𝑀𝐽 ( ˆ 𝚫 ) × 𝑧 ( 𝛼 + 𝛽 ) , ( 𝑚 − 𝑀 ) 𝐽 + ∑︁ even 𝛽> ∑︁ 𝐽 𝐽 𝑀 ∑︁ even 𝛾> ∑︁ 𝑆 𝑆 𝑇 𝜔 𝛽,ℓ𝑚𝛼𝐽 𝐽 𝑀 𝜔 𝛾,𝐽 ( 𝑚 𝑀 ) ( 𝛼 + 𝛽 ) 𝑆 𝑆 𝑇 𝑌 𝑀𝐽 ( ˆ 𝚫 ) 𝑌 𝑇𝑆 ( ˆ 𝚫 ) × 𝑧 ( 𝛼 + 𝛽 + 𝛾 ) , ( 𝑚 − 𝑀 − 𝑇 ) 𝑆 + ... where the first, second and third terms start at zeroth-, second- and fourth-order in 𝜃 respectively. Setting 𝛼 =
0, we obtain a concise form forthe parity-even expansion correct to fourth order:2 𝑌 𝑚ℓ ( ˆ R ) = (cid:2) 𝑌 𝑚ℓ ( ˆ r ) + 𝑌 𝑚ℓ ( ˆ r ) (cid:3) (D5) − ∑︁ even 𝛼> ∑︁ 𝐽 𝐽 𝑀 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 𝑌 𝑀𝐽 ( ˆ 𝚫 ) (cid:104) 𝜖 𝛼 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) + 𝜖 𝛼 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) (cid:105) + ∑︁ even 𝛼> ∑︁ 𝐽 𝐽 𝑀 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 ∑︁ even 𝛽> ∑︁ 𝑆 𝑆 𝑇 Ω 𝛽,𝐽 𝐽 𝑚𝑀 𝛼𝑆 𝑆 𝑇 𝑌 𝑇𝑆 ( ˆ 𝚫 ) (cid:104) 𝜖 𝛼 + 𝛽 𝑌 𝑚 − 𝑀 − 𝑇𝑆 ( ˆ r ) + 𝜖 𝛼 + 𝛽 𝑌 𝑚 − 𝑀 − 𝑇𝑆 ( ˆ r ) (cid:105) + ... MNRAS , 1–27 (2021) eyond Yamamoto Correlators assuming even ℓ , recalling that 𝜔 ,ℓ𝑚𝛼𝐽 𝐽 𝑀 = 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 and contracting the spherical harmonics via (A7). We have additionally inserted thedefinition of Ω from (C3). As before, this can be written as a simple summation over even 𝛼 :2 𝑌 𝑚ℓ ( ˆ R ) = ∑︁ even 𝛼 ∑︁ 𝐽 𝐽 𝑀 Φ 𝛼,ℓ𝑚𝐽 𝐽 𝑀 𝑌 𝑀𝐽 ( ˆ 𝚫 ) (cid:104) 𝜖 𝛼 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) + 𝜖 𝛼 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) (cid:105) , (D6)with Φ ,ℓ𝑚𝐽 𝐽 𝑀 = √ 𝜋𝛿 K 𝐽 𝛿 K 𝐽 ℓ 𝛿 K 𝑀 (D7) Φ ,ℓ𝑚𝐽 𝐽 𝑀 = − 𝜑 ,ℓ𝑚𝐽 𝐽 𝑀 Φ ,ℓ𝑚𝐽 𝐽 𝑀 = − 𝜑 ,ℓ𝑚𝐽 𝐽 𝑀 + ∑︁ 𝐽 (cid:48) 𝐽 (cid:48) 𝑀 (cid:48) 𝜑 ,ℓ𝑚𝐽 (cid:48) 𝐽 (cid:48) 𝑀 (cid:48) Ω , 𝐽 (cid:48) 𝐽 (cid:48) 𝑚𝑀 (cid:48) 𝐽 𝐽 𝑀 . These are equivalent to the Φ coefficients of §3.3, through derived in a different manner. This is straightforward to extend to higher order, andsomewhat faster to compute than the former procedure, since summations are only over even 𝛼 . APPENDIX E: LEGENDRE POLYNOMIAL SHIFT THEOREME1 Derivation
A similar relation to that of Appendix B may be derived for the Legendre polynomial 𝐿 ℓ ( ˆ a · (cid:154) a + A ) . This can be derived in one of two ways;either by expanding in powers of ˆ a · (cid:154) a + A directly, or via the previously proved spherical harmonic shift theorem. For consistency with theabove, we adopt the latter approach, though we note that the two yield consistent results.To perform such an expansion, we first rewrite the Legendre polynomial in terms of spherical harmonics using the addition theorem(A5) and apply (B7), yielding 𝐿 ℓ ( ˆ a · (cid:154) a + A ) = 𝜋 ℓ + ℓ ∑︁ 𝑚 = − ℓ (− ) 𝑚 𝑌 𝑚ℓ ( ˆ a ) 𝑌 𝑚ℓ ( (cid:154) a + A ) (E1) = 𝜋 ℓ + ℓ ∑︁ 𝑚 = − ℓ ∑︁ 𝛼𝐽 𝐽 𝑀 (− ) 𝑚 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 (cid:16) 𝑎𝐴 (cid:17) 𝛼 𝑌 𝑀𝐽 ( ˆ a ) 𝑌 𝑚ℓ ( ˆ a ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ A ) , additionally using 𝑌 𝑚 ∗ ℓ = (− ) 𝑚 𝑌 − 𝑚ℓ . The two spherical harmonics in ˆ a may be contracted using the product-to-sum relation (A7), giving 𝐿 ℓ ( ˆ a · (cid:154) a + A ) = 𝜋 ℓ + ∑︁ 𝛼𝐽 𝐽 𝑀 (cid:48) ∑︁ 𝑚𝐽 (− ) 𝑚 𝜑 𝛼,ℓ𝑚𝐽 𝐽 ( 𝑚 − 𝑀 (cid:48) ) G (− 𝑚 ) ( 𝑚 − 𝑀 (cid:48) ) 𝑀 (cid:48) ℓ𝐽 𝐽 (cid:16) 𝑎𝐴 (cid:17) 𝛼 𝑌 𝑀 (cid:48) ∗ 𝐽 ( ˆ a ) 𝑌 𝑀 (cid:48) 𝐽 ( ˆ A ) (E2) ≡ 𝜋 ℓ + ∑︁ 𝛼𝐽 𝐽 𝑀 (cid:48) 𝑓 𝛼,ℓ𝐽 𝐽 𝑀 (cid:48) (cid:16) 𝑎𝐴 (cid:17) 𝛼 𝑌 𝑀 (cid:48) ∗ 𝐽 ( ˆ a ) 𝑌 𝑀 (cid:48) 𝐽 ( ˆ A ) , where we have relabelled 𝑀 (cid:48) = 𝑚 − 𝑀 , and defined the new coefficient 𝑓 𝛼,ℓ𝐽 𝐽 𝑀 (cid:48) . Whilst this may seem complicated, one can in fact show that 𝑓 𝛼,ℓ𝐽 𝐽 is independent of 𝑀 (cid:48) and non-zero only for 𝐽 = 𝐽 . Defining reduced coefficients 𝑓 𝛼,ℓ𝐽 via 𝑓 𝛼,ℓ𝐽 𝐽 𝑀 (cid:48) ≡ 𝛿 K 𝐽 𝐽 × 𝑓 𝛼,ℓ𝐽 , (E3)and using (A5), this leads to 𝐿 ℓ ( ˆ a · (cid:154) a + A ) = ∞ ∑︁ 𝛼 = ℓ + 𝛼 ∑︁ 𝐽 = max ( ℓ − 𝛼, ) 𝐽 + ℓ + 𝑓 𝛼,ℓ𝐽 (cid:16) 𝑎𝐴 (cid:17) 𝛼 𝐿 𝐽 ( ˆ a · ˆ A ) ; (E4)a far more tractable expansion. Whilst the ( 𝐽 + )/( ℓ + ) factor could be absorbed into 𝑓 𝛼,ℓ𝐽 , we retain it here for later convenience. MNRAS000
A similar relation to that of Appendix B may be derived for the Legendre polynomial 𝐿 ℓ ( ˆ a · (cid:154) a + A ) . This can be derived in one of two ways;either by expanding in powers of ˆ a · (cid:154) a + A directly, or via the previously proved spherical harmonic shift theorem. For consistency with theabove, we adopt the latter approach, though we note that the two yield consistent results.To perform such an expansion, we first rewrite the Legendre polynomial in terms of spherical harmonics using the addition theorem(A5) and apply (B7), yielding 𝐿 ℓ ( ˆ a · (cid:154) a + A ) = 𝜋 ℓ + ℓ ∑︁ 𝑚 = − ℓ (− ) 𝑚 𝑌 𝑚ℓ ( ˆ a ) 𝑌 𝑚ℓ ( (cid:154) a + A ) (E1) = 𝜋 ℓ + ℓ ∑︁ 𝑚 = − ℓ ∑︁ 𝛼𝐽 𝐽 𝑀 (− ) 𝑚 𝜑 𝛼,ℓ𝑚𝐽 𝐽 𝑀 (cid:16) 𝑎𝐴 (cid:17) 𝛼 𝑌 𝑀𝐽 ( ˆ a ) 𝑌 𝑚ℓ ( ˆ a ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ A ) , additionally using 𝑌 𝑚 ∗ ℓ = (− ) 𝑚 𝑌 − 𝑚ℓ . The two spherical harmonics in ˆ a may be contracted using the product-to-sum relation (A7), giving 𝐿 ℓ ( ˆ a · (cid:154) a + A ) = 𝜋 ℓ + ∑︁ 𝛼𝐽 𝐽 𝑀 (cid:48) ∑︁ 𝑚𝐽 (− ) 𝑚 𝜑 𝛼,ℓ𝑚𝐽 𝐽 ( 𝑚 − 𝑀 (cid:48) ) G (− 𝑚 ) ( 𝑚 − 𝑀 (cid:48) ) 𝑀 (cid:48) ℓ𝐽 𝐽 (cid:16) 𝑎𝐴 (cid:17) 𝛼 𝑌 𝑀 (cid:48) ∗ 𝐽 ( ˆ a ) 𝑌 𝑀 (cid:48) 𝐽 ( ˆ A ) (E2) ≡ 𝜋 ℓ + ∑︁ 𝛼𝐽 𝐽 𝑀 (cid:48) 𝑓 𝛼,ℓ𝐽 𝐽 𝑀 (cid:48) (cid:16) 𝑎𝐴 (cid:17) 𝛼 𝑌 𝑀 (cid:48) ∗ 𝐽 ( ˆ a ) 𝑌 𝑀 (cid:48) 𝐽 ( ˆ A ) , where we have relabelled 𝑀 (cid:48) = 𝑚 − 𝑀 , and defined the new coefficient 𝑓 𝛼,ℓ𝐽 𝐽 𝑀 (cid:48) . Whilst this may seem complicated, one can in fact show that 𝑓 𝛼,ℓ𝐽 𝐽 is independent of 𝑀 (cid:48) and non-zero only for 𝐽 = 𝐽 . Defining reduced coefficients 𝑓 𝛼,ℓ𝐽 via 𝑓 𝛼,ℓ𝐽 𝐽 𝑀 (cid:48) ≡ 𝛿 K 𝐽 𝐽 × 𝑓 𝛼,ℓ𝐽 , (E3)and using (A5), this leads to 𝐿 ℓ ( ˆ a · (cid:154) a + A ) = ∞ ∑︁ 𝛼 = ℓ + 𝛼 ∑︁ 𝐽 = max ( ℓ − 𝛼, ) 𝐽 + ℓ + 𝑓 𝛼,ℓ𝐽 (cid:16) 𝑎𝐴 (cid:17) 𝛼 𝐿 𝐽 ( ˆ a · ˆ A ) ; (E4)a far more tractable expansion. Whilst the ( 𝐽 + )/( ℓ + ) factor could be absorbed into 𝑓 𝛼,ℓ𝐽 , we retain it here for later convenience. MNRAS000 , 1–27 (2021) Philcox & Slepian
To see this, we first write out 𝑓 𝛼,ℓ𝐽 𝐽 𝑀 (cid:48) explicitly: 𝑓 𝛼,ℓ𝐽 𝐽 𝑀 (cid:48) = ∑︁ 𝑚𝐽 𝜆𝜇𝑁 𝑘𝐿 (cid:48) 𝐴 𝜆𝜇ℓ𝑚 𝛿 K 𝛼 ( 𝜆 + 𝑁 ) 𝐶 ℓ / 𝑁 𝑘 (− ) 𝑁 − 𝑘 ˜ 𝑐 𝐿 (cid:48) 𝑁 − 𝑘 𝜋 𝐿 (cid:48) + × G 𝜇 ( 𝑚 − 𝑀 (cid:48) − 𝜇 ) ( 𝑀 (cid:48) − 𝑚 ) 𝜆𝐿 (cid:48) 𝐽 G ( 𝑚 − 𝜇 ) ( 𝜇 − 𝑚 + 𝑀 (cid:48) ) (− 𝑀 (cid:48) )( ℓ − 𝜆 ) 𝐿 (cid:48) 𝐽 G (− 𝑚 ) ( 𝑚 − 𝑀 (cid:48) ) 𝑀 (cid:48) ℓ𝐽 𝐽 (− ) 𝑚 − 𝑀 (cid:48) − 𝜇 = ∑︁ 𝜆𝑁 𝑘𝐿 (cid:48) ( ℓ + ) / √︁ ( 𝐽 (cid:48) + )( 𝐽 + ) (cid:18) ℓ 𝜆 (cid:19) / (− ) ℓ 𝛿 K 𝛼 ( 𝜆 + 𝑁 ) 𝐶 ℓ / 𝑁 𝑘 (− ) 𝑁 − 𝑘 ˜ 𝑐 𝐿 (cid:48) 𝑁 − 𝑘 × ∑︁ 𝐽 ( 𝐽 + ) (cid:18) 𝜆 𝐿 (cid:48) 𝐽 (cid:19) (cid:18) ℓ − 𝜆 𝐿 (cid:48) 𝐽 (cid:19) (cid:18) ℓ 𝐽 𝐽 (cid:19) × (cid:34)∑︁ 𝑚𝜇 (− ) 𝑀 (cid:48) − 𝜇 (cid:18) 𝜆 𝐿 (cid:48) 𝐽 𝜇 𝑚 − 𝑀 (cid:48) − 𝜇 𝑀 (cid:48) − 𝑚 (cid:19) (cid:18) ℓ − 𝜆 𝐿 (cid:48) 𝐽 𝑚 − 𝜇 𝜇 − 𝑚 + 𝑀 (cid:48) − 𝑀 (cid:48) (cid:19) (cid:18) ℓ 𝐽 𝐽 − 𝑚 𝑚 − 𝑀 (cid:48) 𝑀 (cid:48) (cid:19) (cid:18) 𝜆 ℓ ℓ − 𝜆𝜇 − 𝑚 𝑚 − 𝜇 (cid:19)(cid:35) , where, in the second line, we have written the Gaunt integrals in terms of Wigner 3- 𝑗 symbols (A8), and noted that 𝐴 𝜆𝜇ℓ𝑚 = (− ) ℓ − 𝑚 √︄ 𝜋 ( ℓ + ) ( 𝜆 + )( ℓ − 𝜆 + ) (cid:18) ℓ 𝜆 (cid:19) / (cid:18) 𝜆 ℓ ℓ − 𝜆𝜇 − 𝑚 𝑚 − 𝜇 (cid:19) . (E6)To simplify, we first employ a summation identity for four 3- 𝑗 symbols, given in Lohmann (2008, Eq. A.24): ∑︁ 𝑚 𝑚 𝑚 𝑚 𝑚 (− ) 𝑗 + 𝑗 + 𝑗 + 𝑚 + 𝑚 + 𝑚 (cid:18) 𝑗 𝑗 𝑗 𝑚 − 𝑚 𝑚 (cid:19) (cid:18) 𝑗 𝑗 𝑗 𝑚 − 𝑚 𝑚 (cid:19) (cid:18) 𝑗 𝑗 𝑗 𝑚 − 𝑚 𝑚 (cid:19) (cid:32) 𝑗 𝑗 𝑗 (cid:48) 𝑚 𝑚 𝑚 (cid:48) (cid:33) = 𝛿 K 𝑗 𝑗 (cid:48) 𝛿 K 𝑚 𝑚 (cid:48) 𝑗 + (cid:26) 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 (cid:27) , (E7)where the quantity in curly parentheses is a Wigner 6- 𝑗 symbol (see NIST DLMF, §34.4). Notably, the RHS of the above equation is independentof 𝑚 and enforces 𝑗 = 𝑗 (cid:48) . Applying this to the final line of (E5), denoted by [ ... ] , with { 𝑗 , 𝑗 , 𝑗 , 𝑗 (cid:48) , 𝑗 , 𝑗 , 𝑗 } = { ℓ, 𝐽 , 𝐽 , 𝐽, 𝐿 (cid:48) , ℓ − 𝜆, 𝜆 } with some massaging via the 3- 𝑗 symmetries gives (cid:2) ... (cid:3) = (− ) 𝐽 + 𝜆 + ℓ 𝛿 K 𝐽 𝐽 𝐽 + (cid:26) ℓ 𝐽 𝐽 𝐿 (cid:48) ℓ − 𝜆 𝜆 (cid:27) , (E8)enforcing 𝐽 = 𝐽 and with no dependence on 𝑀 (cid:48) , as previously mentioned. We may further perform the 𝐽 summation analytically, using the6- 𝑗 summation relation: ∑︁ 𝑗 (− ) 𝑗 + 𝑗 − 𝑗 + 𝑗 + 𝑗 + 𝑗 − 𝑚 − 𝑚 ( 𝑗 + ) (cid:26) 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 (cid:27) (cid:18) 𝑗 𝑗 𝑗 𝑚 𝑚 − 𝑚 (cid:19) (cid:18) 𝑗 𝑗 𝑗 𝑚 𝑚 𝑚 (cid:19) (E9)(Lohmann 2008, Eq. A.26b), which, following application of 6- 𝑗 symmetries, can be used to show ∑︁ 𝐽 ( 𝐽 + )(− ) 𝐽 (cid:18) ℓ 𝐿 (cid:48) 𝐽 (cid:19) (cid:18) ℓ − 𝜆 𝐿 (cid:48) 𝐽 (cid:19) (cid:18) ℓ 𝐽 𝐽 (cid:19) (cid:26) ℓ 𝐽 𝐽𝐿 (cid:48) ℓ − 𝜆 𝜆 (cid:27) (E10) = (− ) 𝐿 (cid:48) + 𝐽 (cid:18) ℓ − 𝜆 𝐿 (cid:48) 𝐽 (cid:19) (cid:18) ℓ 𝜆 ℓ − 𝜆 (cid:19) (cid:18) 𝐿 (cid:48) 𝐽 ℓ − 𝜆 (cid:19) . Inserting into (E5) gives the final form for the reduced coefficients 𝑓 𝛼,ℓ𝐽 : 𝑓 𝛼,ℓ𝐽 = (− ) ℓ ( ℓ + ) / ( 𝛼,ℓ ) ∑︁ 𝜆 = (cid:18) ℓ 𝜆 (cid:19) / (cid:18) ℓ 𝜆 ℓ − 𝜆 (cid:19) (cid:98)( 𝛼 − 𝜆 )/ (cid:99) ∑︁ 𝑘 = 𝐶 ℓ / ( 𝛼 − 𝜆 ) 𝑘 (− ) 𝛼 − 𝜆 − 𝑘 × ∑︁ 𝐿 (cid:48) = ( 𝛼 − 𝜆 − 𝑘 )↓ ˜ 𝑐 𝐿 (cid:48) 𝛼 − 𝜆 − 𝑘 (cid:18) ℓ − 𝜆 𝐿 (cid:48) 𝐽 (cid:19) . (E11)As before, this satisfies 𝑓 ,ℓ𝐽 = 𝛿 K ℓ𝐽 , and 𝑓 𝛼, 𝐽 = 𝛿 K 𝛼 𝛿 K 𝐽 , which, when inserted into (E4) gives 𝐿 ℓ ( ˆ a · (cid:154) a + A ) = 𝐿 ℓ ( ˆ a · ˆ A ) at leading orderand 𝐿 ( ˆ a · (cid:154) a + A ) =
1. Furthermore, we have the parity-rule 𝛼 + ℓ + 𝐽 = even for non-zero 𝑓 𝛼,ℓ𝐽 coefficients. We also note a curious property;the sum of 𝑓 𝛼,ℓ𝐽 over all 𝐽 is equal to zero unless 𝛼 =
0, as demonstrated in (4.3) & (4.11). To prove this, we set A = 𝑘 a in (E4) and recall 𝐿 ℓ ( ) =
1, giving 1 = ∑︁ 𝛼 ∑︁ 𝐽 𝑓 𝛼,ℓ𝐽 (cid:18) 𝑘 (cid:19) 𝛼 = + ∑︁ 𝛼> 𝑘 − 𝛼 ∑︁ 𝐽 𝑓 𝛼,ℓ𝐽 (E12)For this to be true for arbitrary 𝑘 >
1, we require (cid:205) 𝐽 𝑓 𝛼,ℓ𝐽 = 𝛼 > MNRAS , 1–27 (2021) eyond Yamamoto Correlators E2 Parity-Even Form
For the 2PCF estimators considered in §3.3, we use the Legendre polynomial shift theorem (E4) in the following symmetrized form: 𝐿 ℓ ( ˆ 𝚫 · (cid:155) r + r ) = ∞ ∑︁ 𝛼 = ℓ + 𝛼 ∑︁ 𝐽 = max ( ℓ − 𝛼, ) 𝐽 + ℓ + 𝑓 𝛼,ℓ𝐽 (cid:104) 𝜖 𝛼 𝐿 𝐽 ( ˆ 𝚫 · ˆ r ) + (− ) ℓ + 𝐽 𝜖 𝛼 𝐿 𝐽 ( ˆ 𝚫 · ˆ r ) (cid:105) , (E13)for 𝜖 𝑖 ≡ Δ /( 𝑟 𝑖 ) . Just as for the spherical harmonic expansion (§3.3), these may be recast in a manifestly parity-even form, i.e. one onlyinvolving even 𝐽 and 𝛼 (assuming even ℓ ). The derivation of this is analogous to the above, and leads to 𝐿 ℓ ( ˆ 𝚫 · (cid:155) r + r ) = ∞ ∑︁ even 𝛼 = ℓ + 𝛼 ∑︁ 𝐽 = max ( ℓ − 𝛼, ) 𝐽 + ℓ + 𝐹 𝛼,ℓ𝐽 (cid:16) 𝑎𝐴 (cid:17) 𝛼 𝐿 𝐽 ( ˆ a · ˆ A ) , (E14)for 2 𝜖 𝑖 (cid:28)
1, where the parity-even coefficients are defined recursively via 𝐹 𝛼,ℓ𝐽 = 𝑓 𝛼,ℓ𝐽 − ∑︁ odd 𝛽<𝛼 𝛼 − 𝛽 ∑︁ 𝑆 𝑓 𝛽,ℓ𝑆 ℎ ( 𝛼 − 𝛽 ) ,𝛽𝑆𝐽 + ∑︁ odd 𝛽<𝛼 𝛼 − 𝛽 ∑︁ 𝑆 𝑓 𝛽,ℓ𝑆 ∑︁ even 𝛾> ,𝛾<𝛼 − 𝛽 ∑︁ 𝑇 ℎ 𝛾,𝛽𝑆𝑇 ℎ ( 𝛼 − 𝛽 − 𝛾 ) , ( 𝛽 + 𝛾 ) 𝑇𝐽 + ... (E15)using the method of §3.3, or equivalently 𝐹 ,ℓ𝐽 = 𝛿 K ℓ𝐽 (E16) 𝐹 ,ℓ𝐽 = − 𝑓 ,ℓ𝐽 𝐹 ,ℓ𝐽 = − 𝑓 ,ℓ𝐽 + ∑︁ 𝑆 𝑓 ,ℓ𝑆 ℎ , 𝑆𝐽 et cetera , using that of Appendix D. Here, the ℎ 𝛽,𝛼ℓ𝐽 functions are identical to 𝑓 𝛼,ℓ𝐽 (E11), but with ℓ / → ( ℓ + 𝛼 )/ 𝐶 ℓ / 𝑁 𝑘 coefficient.
APPENDIX F: BISECTOR SERIES EXPANSION
To derive the series expansion of §5.1, we start from the bisector vector definition d = 𝑡 r + ( − 𝑡 ) r where 𝑡 = 𝑟 /( 𝑟 + 𝑟 ) , and performthe solid-harmonic expansion (A10): 𝑌 𝑚ℓ ( ˆ d ) = ℓ ∑︁ 𝜆 = 𝜆 ∑︁ 𝜇 = − 𝜆 𝐴 𝜆𝜇ℓ𝑚 (cid:16) 𝑡𝑟 𝑑 (cid:17) ℓ 𝑌 𝜇𝜆 ( ˆ r ) 𝑌 𝑚 − 𝜇ℓ − 𝜆 ( ˆ r ) , (F1)using 𝑡𝑟 = ( − 𝑡 ) 𝑟 , and defining 𝐴 𝜆𝜇ℓ𝑚 as in (B1). Whilst this is not immediately separable in r , r due to the 𝑑 = | d | denominator, wefollow Castorina & White (2018) and note that 𝑑 = 𝑟 𝑟 ( 𝑟 + 𝑟 ) | ˆ r + ˆ r | = 𝑟 𝑡 ( + cos 𝜙 ) , (F2)where cos 𝜙 = ˆ r · ˆ r . Combining with the above: 𝑌 𝑚ℓ ( ˆ d ) = [ ( + cos 𝜙 )] − ℓ / ℓ ∑︁ 𝜆 = 𝜆 ∑︁ 𝜇 = − 𝜆 𝐴 𝜆𝜇ℓ𝑚 𝑌 𝜇𝜆 ( ˆ r ) 𝑌 𝑚 − 𝜇ℓ − 𝜆 ( ˆ r ) . (F3)The former work proceeded to expand the prefactor to second order in cos 𝜙 around cos 𝜙 =
1; here, we note that [ ( + cos 𝜙 )] − ℓ / = − ℓ (cid:20) + cos 𝜙 − (cid:21) − ℓ / = − ℓ ∞ ∑︁ 𝑘 = (cid:18) − ℓ / 𝑘 (cid:19) − 𝑘 ( cos 𝜙 − ) 𝑘 , (F4)via Newton’s generalized binomial theorem around cos 𝜙 = × (cid:0) 𝑎𝑏 (cid:1) = ( 𝑎 ) 𝑏 / 𝑏 !for falling-factorial Pochhammer symbol ( 𝑎 ) 𝑏 , defined in NIST DLMF, §5.2). Using the standard binomial theorem, this gives [ ( + cos 𝜙 )] − ℓ / = − ℓ ∞ ∑︁ 𝑘 = (cid:18) − ℓ / 𝑘 (cid:19) − 𝑘 𝑘 ∑︁ 𝑛 = (cid:18) 𝑘𝑛 (cid:19) (− ) 𝑘 − 𝑛 cos 𝑛 𝜙. (F5)Finally, we may convert cos 𝑛 𝜙 into a sum over Legendre polynomials via (A2) and thus spherical harmonics using (A5). This yields [ ( + cos 𝜙 )] − ℓ / = − ℓ ∞ ∑︁ 𝑘 = (cid:18) − ℓ / 𝑘 (cid:19) − 𝑘 𝑘 ∑︁ 𝑛 = (cid:18) 𝑘𝑛 (cid:19) (− ) 𝑘 − 𝑛 ∑︁ 𝐿 = 𝑛 ↓ ˜ 𝑐 𝐿𝑛 𝜋 𝐿 + 𝐿 ∑︁ 𝑀 = − 𝐿 𝑌 𝑀𝐿 ( ˆ r ) 𝑌 𝑀 ∗ 𝐿 ( ˆ r ) . (F6) MNRAS000
1; here, we note that [ ( + cos 𝜙 )] − ℓ / = − ℓ (cid:20) + cos 𝜙 − (cid:21) − ℓ / = − ℓ ∞ ∑︁ 𝑘 = (cid:18) − ℓ / 𝑘 (cid:19) − 𝑘 ( cos 𝜙 − ) 𝑘 , (F4)via Newton’s generalized binomial theorem around cos 𝜙 = × (cid:0) 𝑎𝑏 (cid:1) = ( 𝑎 ) 𝑏 / 𝑏 !for falling-factorial Pochhammer symbol ( 𝑎 ) 𝑏 , defined in NIST DLMF, §5.2). Using the standard binomial theorem, this gives [ ( + cos 𝜙 )] − ℓ / = − ℓ ∞ ∑︁ 𝑘 = (cid:18) − ℓ / 𝑘 (cid:19) − 𝑘 𝑘 ∑︁ 𝑛 = (cid:18) 𝑘𝑛 (cid:19) (− ) 𝑘 − 𝑛 cos 𝑛 𝜙. (F5)Finally, we may convert cos 𝑛 𝜙 into a sum over Legendre polynomials via (A2) and thus spherical harmonics using (A5). This yields [ ( + cos 𝜙 )] − ℓ / = − ℓ ∞ ∑︁ 𝑘 = (cid:18) − ℓ / 𝑘 (cid:19) − 𝑘 𝑘 ∑︁ 𝑛 = (cid:18) 𝑘𝑛 (cid:19) (− ) 𝑘 − 𝑛 ∑︁ 𝐿 = 𝑛 ↓ ˜ 𝑐 𝐿𝑛 𝜋 𝐿 + 𝐿 ∑︁ 𝑀 = − 𝐿 𝑌 𝑀𝐿 ( ˆ r ) 𝑌 𝑀 ∗ 𝐿 ( ˆ r ) . (F6) MNRAS000 , 1–27 (2021) Philcox & Slepian
Combining the spherical harmonics with those in (F1) via (A7) gives the final form for the expansion: 𝑌 𝑚ℓ ( ˆ d ) = − ℓ ℓ ∑︁ 𝜆 = 𝜆 ∑︁ 𝜇 = − 𝜆 𝐴 𝜆𝜇ℓ𝑚 ∞ ∑︁ 𝑘 = (cid:18) − ℓ / 𝑘 (cid:19) − 𝑘 𝑘 ∑︁ 𝑛 = (cid:18) 𝑘𝑛 (cid:19) (− ) 𝑘 − 𝑛 ∑︁ 𝐿 = 𝑛 ↓ ˜ 𝑐 𝐿𝑛 𝜋 𝐿 + 𝐿 ∑︁ 𝑀 = − 𝐿 (F7) 𝜆 + 𝐿 ∑︁ 𝐽 = | 𝜆 − 𝐿 | (− ) 𝑀 − 𝑚 G 𝜇𝑀 (− 𝜇 − 𝑀 ) 𝜆𝐿𝐽 𝐽 = ℓ − 𝜆 + 𝐿 ∑︁ 𝐽 = | ℓ − 𝜆 − 𝐿 | G ( 𝑚 − 𝜇 ) (− 𝑀 ) ( 𝑀 + 𝜇 − 𝑚 )( ℓ − 𝜆 ) 𝐿𝐽 𝑌 𝜇 + 𝑀𝐽 ( ˆ r ) 𝑌 𝑚 − 𝜇 − 𝑀𝐽 ( ˆ r )⇒ 𝑌 𝑚ℓ ( ˆ d ) ≡ ∞ ∑︁ 𝑘 = ℓ + 𝑘 ∑︁ 𝐽 = ℓ + 𝑘 ∑︁ 𝐽 = 𝐽 ∑︁ 𝑀 = − 𝐽 𝐵 𝑘,ℓ𝑚𝐽 𝐽 𝑀 𝑌 𝑀𝐽 ( ˆ r ) 𝑌 𝑚 − 𝑀𝐽 ( ˆ r ) , defining the coefficients 𝐵 𝑘,ℓ𝑚𝐽 𝐽 𝑀 = − ℓ ℓ ∑︁ 𝜆 = 𝜆 ∑︁ 𝜇 = − 𝜆 𝐴 𝜆𝜇ℓ𝑚 (cid:18) − ℓ / 𝑘 (cid:19) − 𝑘 𝑘 ∑︁ 𝑛 = (cid:18) 𝑘𝑛 (cid:19) (− ) 𝑘 − 𝑛 ∑︁ 𝐿 = 𝑛 ↓ ˜ 𝑐 𝐿𝑛 𝜋 𝐿 + (− ) 𝑀 − 𝜇 − 𝑚 G 𝜇 ( 𝑀 − 𝜇 ) (− 𝑀 ) 𝜆𝐿𝐽 G ( 𝑚 − 𝜇 ) ( 𝜇 − 𝑀 ) ( 𝑀 − 𝑚 )( ℓ − 𝜆 ) 𝐿𝐽 . (F8)The summation limits on 𝐽 , 𝐽 come from the consideration of the Gaunt integral triangle conditions, and the constraints 𝐿 ≤ 𝑛 ≤ 𝑘 , as inAppendix B2. Furthermore, the Gaunt integrals require even 𝜆 + 𝐿 + 𝐽 and ℓ − 𝜆 + 𝐿 + 𝐽 , thus ℓ + 𝐽 + 𝐽 is even, but we do not require 𝐽 and 𝐽 themselves to be even.For the special case of 𝑘 =
0, we obtain 𝐵 ,ℓ𝑚𝐽 𝐽 𝑀 = − ℓ ℓ ∑︁ 𝜆 = 𝜆 ∑︁ 𝜇 = − 𝜆 𝐴 𝜆𝜇ℓ𝑚 ( 𝜋 )(− ) 𝑀 − 𝜇 − 𝑚 G 𝜇 ( 𝑀 − 𝜇 ) (− 𝑀 ) 𝜆 𝐽 G ( 𝑚 − 𝜇 ) ( 𝜇 − 𝑀 ) ( 𝑀 − 𝑚 )( ℓ − 𝜆 ) 𝐽 (F9) = 𝛿 K 𝐽 ( ℓ − 𝐽 ) × − ℓ 𝐴 𝐽 𝑀ℓ𝑚 . using (cid:0) − ℓ / (cid:1) =
1, ˜ 𝑐 = 𝐽 , 𝐽 ≤ ℓ and | 𝑀 | ≤ 𝐽 . Notably, 𝑌 𝑚ℓ ( ˆ d ) still involves asummation over 𝐽 , 𝐽 , 𝑀 at zeroth-order in 𝑘 . Considering ℓ = 𝑚 =
0, the coefficient becomes 𝐵 𝑘, 𝐽 𝐽 𝑀 = 𝐴 − 𝑘 𝑘 ∑︁ 𝑛 = (cid:18) 𝑘𝑛 (cid:19) (− ) 𝑘 − 𝑛 ∑︁ 𝐿 = 𝑛 ↓ ˜ 𝑐 𝐿𝑛 𝜋 𝐿 + (− ) 𝑀 G 𝑀 (− 𝑀 ) 𝐿𝐽 G (− 𝑀 ) 𝑀 𝐿𝐽 (F10) = 𝛿 K 𝐽 𝛿 K 𝐽 𝛿 K 𝑀 × √ 𝜋 using 𝐴 = √ 𝜋 and (cid:0) 𝑘 (cid:1) = 𝛿 K 𝑘 . As expected, this recovers 𝑌 ( ˆ d ) = /√ 𝜋 . REFERENCES
Abramowitz M., Stegun I. A., 1964, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover, New York CityAlam S., et al., 2017, MNRAS, 470, 2617Arfken G., Weber H., Harris F., 2013, Mathematical Methods for Physicists: A Comprehensive Guide. Elsevier ScienceBeutler F., et al., 2012, MNRAS, 423, 3430Beutler F., et al., 2017, MNRAS, 466, 2242Beutler F., Castorina E., Zhang P., 2019, J. Cosmology Astropart. Phys., 2019, 040Bianchi D., Gil-Marín H., Ruggeri R., Percival W. J., 2015, MNRAS, 453, L11Bonvin C., Durrer R., 2011, Phys. Rev. D, 84, 063505Castorina E., White M., 2018, MNRAS, 476, 4403Chudaykin A., Ivanov M. M., 2019, J. Cosmology Astropart. Phys., 2019, 034DESI Collaboration et al., 2016, arXiv e-prints, p. arXiv:1611.00036Datta K. K., Choudhury T. R., Bharadwaj S., 2007, MNRAS, 378, 119Doré O., et al., 2014, arXiv e-prints, p. arXiv:1412.4872Eisenstein D. J., et al., 2011, AJ, 142, 72Feldman H. A., Kaiser N., Peacock J. A., 1994, ApJ, 426, 23Garcia K., Slepian Z., 2020, arXiv e-prints, p. arXiv:2011.03503Gil-Marín H., Percival W. J., Verde L., Brownstein J. R., Chuang C.-H., Kitaura F.-S., Rodríguez-Torres S. A., Olmstead M. D., 2017, MNRAS, 465, 1757Gradshteyn I. S., Ryzhik I. M., 1994, Table of integrals, series and products. Elsevier ScienceHamilton A. J. S., 1992, ApJ, 385, L5Hamilton A. J. S., 1998, Linear Redshift Distortions: a Review. p. 185, doi:10.1007/978-94-011-4960-0_17Hamilton A. J. S., Culhane M., 1996, MNRAS, 278, 73Hand N., Li Y., Slepian Z., Seljak U., 2017, J. Cosmology Astropart. Phys., 2017, 002Hand N., Feng Y., Beutler F., Li Y., Modi C., Seljak U., Slepian Z., 2018, AJ, 156, 160Kaiser N., 1987, MNRAS, 227, 1Kitaura F.-S., et al., 2016, MNRAS, 456, 4156 MNRAS , 1–27 (2021) eyond Yamamoto Correlators Landy S. D., Szalay A. S., 1993, ApJ, 412, 64Laureijs R., et al., 2011, arXiv e-prints, p. arXiv:1110.3193Lesgourgues J., Pastor S., 2006, Phys. Rep., 429, 307Lohmann B., 2008, Angle and spin resolved Auger emission: Theory and applications to atoms and molecules. Vol. 46, Springer Science & Business MediaMatsubara T., 2000, ApJ, 535, 1NIST DLMF, NIST Digital Library of Mathematical Functions. http://dlmf.nist.gov/
Pápai P., Szapudi I., 2008, MNRAS, 389, 292Philcox O. H. E., 2020, arXiv e-prints, p. arXiv:2012.09389Philcox O. H. E., 2021, MNRAS, 501, 4004Philcox O. H. E., Eisenstein D. J., 2020, MNRAS, 492, 1214Philcox O. H. E., Massara E., Spergel D. N., 2020, Phys. Rev. D, 102, 043516Raccanelli A., Bertacca D., Doré O., Maartens R., 2014, J. Cosmology Astropart. Phys., 2014, 022Reimberg P., Bernardeau F., Pitrou C., 2016, J. Cosmology Astropart. Phys., 2016, 048Rodríguez-Torres S. A., et al., 2016, MNRAS, 460, 1173Samushia L., Percival W. J., Raccanelli A., 2012, MNRAS, 420, 2102Samushia L., Branchini E., Percival W. J., 2015, MNRAS, 452, 3704Samushia L., Slepian Z., Villaescusa-Navarro F., 2021, arXiv e-prints, p. arXiv:2102.01696Scoccimarro R., 2015, Phys. Rev. D, 92, 083532Shaw J. R., Lewis A., 2008, Phys. Rev. D, 78, 103512Sinha M., Garrison L. H., 2020, MNRAS, 491, 3022Slepian Z., Eisenstein D. J., 2015, arXiv e-prints, p. arXiv:1510.04809Slepian Z., Eisenstein D. J., 2016, MNRAS, 455, L31Slepian Z., et al., 2017, MNRAS, 469, 1738Sugiyama N. S., Saito S., Beutler F., Seo H.-J., 2019, MNRAS, 484, 364Szalay A. S., Matsubara T., Landy S. D., 1998, ApJ, 498, L1Szapudi I., 2004, ApJ, 614, 51Tough R. J. A., Stone A. J., 1977, Journal of Physics A: Mathematical and General, 10, 1261Weinberg D. H., Mortonson M. J., Eisenstein D. J., Hirata C., Riess A. G., Rozo E., 2013, Phys. Rep., 530, 87Yamamoto K., Nakamichi M., Kamino A., Bassett B. A., Nishioka H., 2006, PASJ, 58, 93Yoo J., Seljak U., 2015, MNRAS, 447, 1789Zaroubi S., Hoffman Y., 1996, ApJ, 462, 25eBOSS Collaboration et al., 2020, arXiv e-prints, p. arXiv:2007.08991This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000