Using galaxy pairs to investigate the three-point correlation function in the squeezed limit
MMNRAS , 1–16 (2017) Preprint 13 September 2018 Compiled using MNRAS L A TEX style file v3.0
Using galaxy pairs to investigate the three-pointcorrelation function in the squeezed limit
Sihan Yuan, (cid:63) Daniel J. Eisenstein, and Lehman H. Garrison Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02138, USA
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We investigate the three-point correlation function (3PCF) in the squeezed limit byconsidering galaxy pairs as discrete objects and cross-correlating them with the galaxyfield. We develop an efficient algorithm using Fast Fourier Transforms to compute suchcross-correlations and their associated pair-galaxy bias b p , g and the squeezed 3PCFcoefficient Q eff . We implement our method using N-body cosmological simulations anda fiducial Halo Occupation Distribution (HOD) and present the results in both thereal space and redshift space. In real space, we observe a peak in b p , g and Q eff at pairseparation of ∼ b p , g and Q eff signals that track the large-scale filamentary structure. In redshift space,both the 2 Mpc peak and the anisotropy are significantly smeared out along the line-of-sight due to Finger-of-God effect. In both the real space and redshift space, thesqueezed 3PCF shows a factor of 2 variation, contradicting the hierarchical ansatzbut offering rich information on the galaxy-halo connection. Thus, we explore thepossibility of using the squeezed 3PCF to constrain the HOD. When we compare twosimple HOD models that are closely matched in their projected two-point correlationfunction (2PCF), we do not yet see a strong variation in the 3PCF that is clearlydisentangled from variations in the projected 2PCF. Nevertheless, we propose thatmore complicated HOD models, e.g. those incorporating assembly bias, can breakdegeneracies in the 2PCF and show a distinguishable squeezed 3PCF signal. Key words: cosmology: large-scale structure of Universe – cosmology: dark matter– galaxies: haloes – methods: analytical
Understanding galaxy clustering is one of the primary goalsof cosmology, and the advent of a new generation of galaxysurveys has ushered in the new era of high precision galaxyclustering studies. A key objective of such studies is to un-derstand the relationship between galaxies and the under-lying dark matter halos since such relationship is critical inhelping us infer the linear dark matter power spectrum fromthe observed galaxy power spectrum. The Halo OccupationDistribution (HOD) is the most popular framework for de-scribing how the number of galaxies per halo depends thehalo mass (Seljak 2000; Peacock & Smith 2000; Scoccimarroet al. 2001; Berlind & Weinberg 2002; Cooray & Sheth 2002;Kravtsov et al. 2004; Zheng et al. 2005). For the simple HODwe consider in this paper, galaxies populating a halo are di-vided into central and satellite galaxies, hereafter referred to (cid:63)
E-mail: [email protected] simply as centrals and satellites (Berlind et al. 2003; Zhenget al. 2005).There has been a wealth of studies that utilize the two-point correlation functions (2PCF) or higher-order statisticsto constrain the HOD (e.g. Zehavi et al. 2005b; Kulkarniet al. 2007; Blake et al. 2008; Zheng et al. 2009; Whiteet al. 2011; Parejko et al. 2013; More et al. 2015; Guoet al. 2016; Rodr´ıguez-Torres et al. 2016). Another popularway to constrain the HOD is by probing the halo mass ofgalaxies via weak gravitational lensing, commonly know asgalaxy-galaxy lensing, (Mandelbaum et al. 2006; Yoo et al.2006; Li et al. 2009; Leauthaud et al. 2012; Tinker et al.2012; Mandelbaum et al. 2013; Miyatake et al. 2015; Moreet al. 2015). Other methods for constraining the HOD in-clude a pair-counting based technique called the Counts-in-Cylinders (CiC) (Reid & Spergel 2009), and direct countingof the number of galaxies in clusters as a function of halomass (Ho et al. 2009; Hoshino et al. 2015).For this paper, we are interested in constraining the © a r X i v : . [ a s t r o - ph . C O ] S e p S. Yuan, D. J. Eisenstein and L. H. Garrison high mass regime of the HOD ( M ∼ − M (cid:12) ), where eachdark matter halo hosts multiple galaxies close together. Thissuggests that if we consider close pairs of galaxies as objectsin their own right, they would be biased tracers of high masshalos (refer to Figure 4).The first step is to formulate the proper statistic on closegalaxy pairs to calculate their bias. One might choose to usethe pair-pair auto-correlation function. However, to preservethe pair separation dependence of our statistic, we need toauto-correlate much smaller subsamples of pairs binned bypair separation. Auto-correlating smaller samples leads toa noisier statistic. Thus, instead of auto-correlating closepairs, we investigate the cross-correlation of close pairs withgalaxy singlets to calculate the pair-galaxy cross-correlationfunction. The large-scale pair-galaxy cross-correlation thenreveals the linear bias (Fry 1996; Tegmark & Peebles 1998;Scherrer & Weinberg 1998) of the pair field relative to thegalaxy field, i.e. the pair-galaxy bias, which we denote as b p , g , given by b p , g = ξ p , g / ξ g , g , (1)where ξ p , g and ξ g , g are the pair-galaxy cross-correlationand the galaxy-galaxy cross-correlation, respectively. Whenthe pair separation is much less than the distance to thegalaxy singlet, the pair-galaxy cross-correlation is equiva-lent to the three-point correlation function (3PCF) in thesqueezed limit. For conciseness, we refer to the 3PCF in thesqueezed limit as the squeezed 3PCF in the rest of the paper.Generally speaking, the 3PCF, ζ ( r , r , r ) , describesthe probability distribution of finding three galaxies forminga triangle with three sides given by vectors r , r , and r (SeeBernardeau et al. (2002) for a comprehensive review). In thehierarchy of correlation functions, the 3PCF is the secondlowest order correlation function after the 2PCF, and is thelowest order correlation function to provide shape informa-tion (Peebles 1980). The squeezed limit of the 3PCF refersto the limit where one side of the triangle is much shorterthan the other two sides, i.e. r (cid:28) r ≈ r .The spatial distribution of matter in the Universe canbe fully described by the 2PCF if the matter density field isfully Gaussian. However, we expect non-linear gravitationalevolution to produce non-Gaussian signatures in the galaxydistribution that we measure today (Bernardeau et al. 2002).The 3PCF is thus needed to complement the informationprovided by the 2PCF when describing late-time galaxy dis-tribution. A series of recent papers have focused on mea-suring and analysing the 3PCF in modern cosmological sur-veys, such as Gazta˜naga et al. (2009); Mar´ın (2011); Gil-Mar´ın et al. (2015); Slepian et al. (2017a,b). The additionalinformation provided by the 3PCF is often used to breakmodel degeneracies describing cosmology and galaxy bias(Gaztanaga & Frieman 1994; Jing & B¨orner 2004; Kulkarniet al. 2007; Zheng & Weinberg 2007; McBride et al. 2011).Additionally, the 3PCF and its Fourier space counterpart,the bispectrum, are also commonly used to probe primordialnon-gaussianities as a test of single-field slow-roll inflation-ary models (e.g. Dalal et al. 2008; Fergusson & Shellard 2009;Sefusatti et al. 2010; Baldauf et al. 2011; Tasinato et al. 2014;Hashimoto et al. 2016). Tellarini et al. (2016) offers a nicereview of this topic.Measuring the 3PCF is computationally intense as thenumber of triplets scales as N and thus becomes prohibitive for large samples of galaxies. The 3PCF in the squeezedlimit is advantageous because its computation can be re-duced to two separate pair finding problems of complexity O ( N ) . Thus, even though we phrased our statistic as a cross-correlation between the pair field and the galaxy field, thesame statistic can be interpreted as a squeezed 3PCF. Weshow the equivalence of the two approaches in Section 2.2.In the following sections of the paper, we present a fastand efficient way of calculating the pair-galaxy bias usingFast Fourier Transforms (FFTs) in Section 2.1, and how ittranslates to the squeezed 3PCF in Section 2.2. We discussimplementing our methodology for Luminous Red Galax-ies (LRGs) using N-body cosmological simulations in Sec-tion 3, and present the results in Section 4. In Section 5, wediscuss the possibility of constraining HOD with the pair-galaxy bias and the squeezed 3PCF. Finally, we draw a fewconcluding remarks in Section 6. In this section, we develop a fast and efficient method ofcomputing b p , g using FFT and then formulate the exactrelationship between b p , g and the squeezed 3PCF startingfrom the triplet counting definition of the 3PCF. b p , g Consider a galaxy overdensity field δ g ( x ) , a galaxy pair over-density field δ p ( x ) , and their Fourier space counterparts ˜ δ g ( k ) and ˜ δ p ( k ) . Assuming white noise denoted by σ , fora single wave vector k , we can write out a χ objectivefunction χ = | ˜ δ p ( k ) − b p , g ˜ δ g ( k )| σ . (2)Then by minimizing χ against b p , g , we find a best fit for b p , g and its uncertainties at wave vector k b p , g ( k ) = Re ( ˜ δ g ∗ ( k ) ˜ δ p ( k )) ˜ δ g ∗ ( k ) ˜ δ g ( k ) , σ − b ( k ) = ˜ δ g ∗ ( k ) ˜ δ g ( k ) σ . (3)Then we average over a range of k modes, weighted by in-verse variance multiplied with some weighting function ˜ W ( k ) .We get b p , g = (cid:205) k b p , g ( k ) σ − b ( k ) ˜ W ( k ) (cid:205) k σ − b ( k ) ˜ W ( k ) = (cid:205) k Re ( ˜ δ g ∗ ( k ) ˜ δ p ( k )) ˜ W ( k ) (cid:205) k ˜ δ g ∗ ( k ) ˜ δ g ( k ) ˜ W ( k ) = ∫ d k ˜ W ( k ) Re ( ˜ δ g ∗ ( k ) ˜ δ p ( k )) ∫ d k ˜ W ( k ) P g ( k ) , (4)where P g ( k ) = ˜ δ g ∗ ( k ) ˜ δ g ( k ) is the galaxy power spectrum. Adetailed discussion on the choice of weighting function ˜ W ( k ) is provided in Section 3.6.Equation 4 shows that b p , g is linear in δ p , which meansthat we can compute the integrals for a single pair at atime and then sum over all pairs to obtain b p , g . For the i -thgalaxy pair, we can define its position x p , i as the middle of MNRAS000
E-mail: [email protected] simply as centrals and satellites (Berlind et al. 2003; Zhenget al. 2005).There has been a wealth of studies that utilize the two-point correlation functions (2PCF) or higher-order statisticsto constrain the HOD (e.g. Zehavi et al. 2005b; Kulkarniet al. 2007; Blake et al. 2008; Zheng et al. 2009; Whiteet al. 2011; Parejko et al. 2013; More et al. 2015; Guoet al. 2016; Rodr´ıguez-Torres et al. 2016). Another popularway to constrain the HOD is by probing the halo mass ofgalaxies via weak gravitational lensing, commonly know asgalaxy-galaxy lensing, (Mandelbaum et al. 2006; Yoo et al.2006; Li et al. 2009; Leauthaud et al. 2012; Tinker et al.2012; Mandelbaum et al. 2013; Miyatake et al. 2015; Moreet al. 2015). Other methods for constraining the HOD in-clude a pair-counting based technique called the Counts-in-Cylinders (CiC) (Reid & Spergel 2009), and direct countingof the number of galaxies in clusters as a function of halomass (Ho et al. 2009; Hoshino et al. 2015).For this paper, we are interested in constraining the © a r X i v : . [ a s t r o - ph . C O ] S e p S. Yuan, D. J. Eisenstein and L. H. Garrison high mass regime of the HOD ( M ∼ − M (cid:12) ), where eachdark matter halo hosts multiple galaxies close together. Thissuggests that if we consider close pairs of galaxies as objectsin their own right, they would be biased tracers of high masshalos (refer to Figure 4).The first step is to formulate the proper statistic on closegalaxy pairs to calculate their bias. One might choose to usethe pair-pair auto-correlation function. However, to preservethe pair separation dependence of our statistic, we need toauto-correlate much smaller subsamples of pairs binned bypair separation. Auto-correlating smaller samples leads toa noisier statistic. Thus, instead of auto-correlating closepairs, we investigate the cross-correlation of close pairs withgalaxy singlets to calculate the pair-galaxy cross-correlationfunction. The large-scale pair-galaxy cross-correlation thenreveals the linear bias (Fry 1996; Tegmark & Peebles 1998;Scherrer & Weinberg 1998) of the pair field relative to thegalaxy field, i.e. the pair-galaxy bias, which we denote as b p , g , given by b p , g = ξ p , g / ξ g , g , (1)where ξ p , g and ξ g , g are the pair-galaxy cross-correlationand the galaxy-galaxy cross-correlation, respectively. Whenthe pair separation is much less than the distance to thegalaxy singlet, the pair-galaxy cross-correlation is equiva-lent to the three-point correlation function (3PCF) in thesqueezed limit. For conciseness, we refer to the 3PCF in thesqueezed limit as the squeezed 3PCF in the rest of the paper.Generally speaking, the 3PCF, ζ ( r , r , r ) , describesthe probability distribution of finding three galaxies forminga triangle with three sides given by vectors r , r , and r (SeeBernardeau et al. (2002) for a comprehensive review). In thehierarchy of correlation functions, the 3PCF is the secondlowest order correlation function after the 2PCF, and is thelowest order correlation function to provide shape informa-tion (Peebles 1980). The squeezed limit of the 3PCF refersto the limit where one side of the triangle is much shorterthan the other two sides, i.e. r (cid:28) r ≈ r .The spatial distribution of matter in the Universe canbe fully described by the 2PCF if the matter density field isfully Gaussian. However, we expect non-linear gravitationalevolution to produce non-Gaussian signatures in the galaxydistribution that we measure today (Bernardeau et al. 2002).The 3PCF is thus needed to complement the informationprovided by the 2PCF when describing late-time galaxy dis-tribution. A series of recent papers have focused on mea-suring and analysing the 3PCF in modern cosmological sur-veys, such as Gazta˜naga et al. (2009); Mar´ın (2011); Gil-Mar´ın et al. (2015); Slepian et al. (2017a,b). The additionalinformation provided by the 3PCF is often used to breakmodel degeneracies describing cosmology and galaxy bias(Gaztanaga & Frieman 1994; Jing & B¨orner 2004; Kulkarniet al. 2007; Zheng & Weinberg 2007; McBride et al. 2011).Additionally, the 3PCF and its Fourier space counterpart,the bispectrum, are also commonly used to probe primordialnon-gaussianities as a test of single-field slow-roll inflation-ary models (e.g. Dalal et al. 2008; Fergusson & Shellard 2009;Sefusatti et al. 2010; Baldauf et al. 2011; Tasinato et al. 2014;Hashimoto et al. 2016). Tellarini et al. (2016) offers a nicereview of this topic.Measuring the 3PCF is computationally intense as thenumber of triplets scales as N and thus becomes prohibitive for large samples of galaxies. The 3PCF in the squeezedlimit is advantageous because its computation can be re-duced to two separate pair finding problems of complexity O ( N ) . Thus, even though we phrased our statistic as a cross-correlation between the pair field and the galaxy field, thesame statistic can be interpreted as a squeezed 3PCF. Weshow the equivalence of the two approaches in Section 2.2.In the following sections of the paper, we present a fastand efficient way of calculating the pair-galaxy bias usingFast Fourier Transforms (FFTs) in Section 2.1, and how ittranslates to the squeezed 3PCF in Section 2.2. We discussimplementing our methodology for Luminous Red Galax-ies (LRGs) using N-body cosmological simulations in Sec-tion 3, and present the results in Section 4. In Section 5, wediscuss the possibility of constraining HOD with the pair-galaxy bias and the squeezed 3PCF. Finally, we draw a fewconcluding remarks in Section 6. In this section, we develop a fast and efficient method ofcomputing b p , g using FFT and then formulate the exactrelationship between b p , g and the squeezed 3PCF startingfrom the triplet counting definition of the 3PCF. b p , g Consider a galaxy overdensity field δ g ( x ) , a galaxy pair over-density field δ p ( x ) , and their Fourier space counterparts ˜ δ g ( k ) and ˜ δ p ( k ) . Assuming white noise denoted by σ , fora single wave vector k , we can write out a χ objectivefunction χ = | ˜ δ p ( k ) − b p , g ˜ δ g ( k )| σ . (2)Then by minimizing χ against b p , g , we find a best fit for b p , g and its uncertainties at wave vector k b p , g ( k ) = Re ( ˜ δ g ∗ ( k ) ˜ δ p ( k )) ˜ δ g ∗ ( k ) ˜ δ g ( k ) , σ − b ( k ) = ˜ δ g ∗ ( k ) ˜ δ g ( k ) σ . (3)Then we average over a range of k modes, weighted by in-verse variance multiplied with some weighting function ˜ W ( k ) .We get b p , g = (cid:205) k b p , g ( k ) σ − b ( k ) ˜ W ( k ) (cid:205) k σ − b ( k ) ˜ W ( k ) = (cid:205) k Re ( ˜ δ g ∗ ( k ) ˜ δ p ( k )) ˜ W ( k ) (cid:205) k ˜ δ g ∗ ( k ) ˜ δ g ( k ) ˜ W ( k ) = ∫ d k ˜ W ( k ) Re ( ˜ δ g ∗ ( k ) ˜ δ p ( k )) ∫ d k ˜ W ( k ) P g ( k ) , (4)where P g ( k ) = ˜ δ g ∗ ( k ) ˜ δ g ( k ) is the galaxy power spectrum. Adetailed discussion on the choice of weighting function ˜ W ( k ) is provided in Section 3.6.Equation 4 shows that b p , g is linear in δ p , which meansthat we can compute the integrals for a single pair at atime and then sum over all pairs to obtain b p , g . For the i -thgalaxy pair, we can define its position x p , i as the middle of MNRAS000 , 1–16 (2017)
PCF in squeezed limit the pair. Then the fractional overdensity function of pairs isgiven by δ p , i ( x ) = δ D ( x − x p , i ) ¯ n p − , (5)where δ D ( x ) is the Dirac delta function. ¯ n p denotes the meannumber density of galaxy pairs. The Fourier Transform isgiven by ˜ δ p , i ( k ) = e − i k · x p , i ¯ n p − δ D ( k )( π ) (6)Then we calculate the value Re ( ˜ δ g ∗ ( k ) ˜ δ p ( k )) in equation 4asRe ( ˜ δ g ∗ ( k ) ˜ δ p ( k )) = Re ( ˜ δ g ∗ ( k ) (cid:213) i ˜ δ p , i ( k )) = (cid:213) i ( ˜ δ g ∗ ( k ) ˜ δ p , i ( k ) + ˜ δ g ( k ) ˜ δ ∗ p , i ( k )) =
12 ¯ n p (cid:213) i ( ˜ δ g ∗ ( k ) e − i k · x p , i + ˜ δ g ( k ) e i k · x p , i )− N p δ D ( k )( π ) ( ˜ δ g ( k ) + ˜ δ g ∗ ( k )) , (7)where N p is the total number of galaxy pairs. When we plugthe second term of equation 7 into equation 4, we get a resultproportional to ( ˜ δ g ( ) + ˜ δ g ∗ ( )) , which gives ˜ δ g ( ) = ˜ δ g ∗ ( ) = ∫ d x δ g ( x ) = . (8)Thus, the second term of equation 7 does not contribute tothe numerator in equation 4.By combining equation 7 and equation 4, we get b p , g (cid:20)∫ d k ˜ W ( k ) P g ( k ) (cid:21) =
12 ¯ n p (cid:213) i (cid:20)∫ d k ˜ W ( k )( ˜ δ g ∗ ( k ) e − i k · x p , i + ˜ δ g ( k ) e i k · x p , i ) (cid:21) =
12 ¯ n p (cid:213) i (cid:20)∫ d k ˜ W (− k ) ˜ δ g (− k ) e − i k · x p , i + ∫ d k ˜ W ( k ) ˜ δ g ( k ) e i k · x p , i (cid:21) = n p (cid:213) i (cid:20)∫ d k ˜ W ( k ) ˜ δ g ( k ) e i k · x p , i (cid:21) = n p (cid:213) i ( W ∗ δ g )( x p , i ) . (9)Note that during this derivation, we invoke that W ( x ) is evenand real and that δ g ( x ) is real. Similarly, we can evaluatethe denominator of equation 4 ∫ d k ˜ W ( k ) P g ( k ) = n g (cid:213) j ( W ∗ δ g )( x g , j ) , (10)where we sum over all the galaxy positions x g , j , and ¯ n g isthe mean galaxy number density.Thus, we get the final expression for pair-galaxy bias b p , g = n p (cid:205) i ( W ∗ δ g )( x p , i ) n g (cid:205) j ( W ∗ δ g )( x g , j ) . (11)This formula gives a fast and efficient way of calculating b p , g in linear time as a function of the number of galaxy pairs.The numerator allows us to compute only one convolutionfor the full pair sample and simply evaluate it at the po-sitions of any subsample of pairs whereas the denominatoralso only requires one convolution. This is much faster thanexplicitly cross-correlating each subsample of pairs with theoverall galaxy field. b p , g to the the squeezed limit of the3PCF To relate b p , g to the squeezed 3PCF, we first evaluate thedenominator of equation 11 n g (cid:213) j ( W ∗ δ g )( x g , j ) = n g (cid:213) j ∫ d x (cid:48) δ g ( x (cid:48) + x g , j ) W ( x (cid:48) ) = n g (cid:213) j ∫ d x (cid:48) (cid:18) n g ( x (cid:48) + x g , j | x g , j ) ¯ n g − (cid:19) W ( x (cid:48) ) , (12)where n g ( x (cid:48) + x g , j | x g , j ) is the number of galaxies at x (cid:48) + x g , j given an arbitrary galaxy at x g , j . By definition of 2PCF ξ ( x (cid:48) ) , we have n g ( x (cid:48) + x g , j | x g , j ) = ¯ n g ( + ξ ( x (cid:48) )) (13)Thus, combining equation 12 and equation 13, we have n g (cid:213) j ( W ∗ δ g )( x g , j ) = ∫ d x (cid:48) ξ ( x (cid:48) ) W ( x (cid:48) ) . (14)Similarly to equation 12, we can evaluate the numeratorof equation 11 as n p (cid:213) i ( W ∗ δ g )( x p , i ) = n p (cid:213) i ∫ d x (cid:48) (cid:18) n g ( x (cid:48) + x p , i | x p , i ) ¯ n g − (cid:19) W ( x (cid:48) ) , (15)where n g ( x (cid:48) + x p , i | x p , i ) is the number of galaxies at x (cid:48) + x p , i given an arbitrary galaxy pair at position x p , i . This numberrelates to the 3PCF by the following n g ( x (cid:48) + x p , i | x p , i ) = ¯ n g ( + ξ ( d p ) + ξ ( x (cid:48) + d p ) + ξ ( x (cid:48) − d p ) + ζ ( d p , x (cid:48) + d p , x (cid:48) − d p ))/( + ξ ( d p )) , (16)where d p is the separation vector of the galaxy pair. Now wecombine equation 15 and equation 16 in the squeezed limit.Then, we have n p (cid:213) i ( W ∗ δ g )( x p , i ) = n p (cid:213) i ∫ d x (cid:48) (cid:18) ξ ( x (cid:48) ) + ξ ( x (cid:48) ) + ζ ( d p , x (cid:48) , x (cid:48) ) + ξ ( d p ) (cid:19) W ( x (cid:48) ) (17)where we have defined x (cid:48) = x (cid:48) + d p / and x (cid:48) = x (cid:48) − d p / .For our purpose, we are only summing over pairs withina specific d p bin. Thus, to good approximation, we can treat d p as approximately constant in the summation, which re-moves all dependence on pair position inside the summation. MNRAS , 1–16 (2017)
S. Yuan, D. J. Eisenstein and L. H. Garrison
Then we have n p (cid:213) i ( W ∗ δ g )( x p , i ) = ∫ d x (cid:48) W ( x (cid:48) )( ξ ( x (cid:48) ) + ξ ( x (cid:48) )) + ξ ( d p ) + ∫ d x (cid:48) W ( x (cid:48) ) ζ ( d p , x (cid:48) , x (cid:48) ) + ξ ( d p ) . (18)Then, using equation 10, equation 11, and equation 18, wecan express the pair-galaxy bias in terms of the 2PCF andthe 3PCF b p , g = A + B + ξ ( d p ) , (19)where A = ∫ d x (cid:48) W ( x (cid:48) )( ξ ( x (cid:48) ) + ξ ( x (cid:48) )) ∫ d x (cid:48) W ( x (cid:48) ) ξ ( x (cid:48) ) , B = ∫ d x (cid:48) W ( x (cid:48) ) ζ ( d p , x (cid:48) , x (cid:48) ) ∫ d x (cid:48) W ( x (cid:48) ) ξ ( x (cid:48) ) . So now we have expressed the pair-galaxy bias as given byequation 11 in terms of general case 2PCFs and 3PCFs.To proceed, we introduce the reduced 3PCF Q eff , de-fined by Q eff ( x , x , x ) = ζ ( x , x , x ) ξ ( x ) ξ ( x ) + ξ ( x ) ξ ( x ) + ξ ( x ) ξ ( x ) , (20)where x , x , x are the three separation vectors for any3PCF. In the original hierarchical model (Peebles & Groth1975; Soneira & Peebles 1976; Groth & Peebles 1977; Peebles1980), the Q eff factor is a constant, which has been disprovenin later surveys and simulations (e.g. Kulkarni et al. 2007;Mar´ın et al. 2008; McBride et al. 2011). As we will see, ourresults in the squeezed limit also show a strongly variable Q eff .In the strictly squeezed limit, | x (cid:48) | (cid:29) | d p | and ξ ( x (cid:48) ) ≈ ξ ( x (cid:48) ) ≈ ξ ( x (cid:48) ) . A typical 2PCF assumes a power-law form ξ ( x (cid:48) ) ∝ x (cid:48) α , where the exponent α ≈ − . This means that, inthe strictly squeezed limit, ξ ( x (cid:48) ) (cid:28) ξ ( d p ) , and thus ζ ( d p , x (cid:48) , x (cid:48) ) ≈ Q eff ( ξ ( x (cid:48) ) ξ ( d p ) + ξ ( x (cid:48) ) )≈ Q eff ξ ( x (cid:48) ) ξ ( d p ) . (21)Then combining equation 19 and equation 21, we expresspair-galaxy bias as b p , g = ( + Q eff ξ ( d p )) + ξ ( d p ) . (22)Then we can express the reduced 3PCF Q eff as a functionof pair-galaxy bias b p , g Q eff = ( + ξ ( d p )) b p , g − ξ ( d p ) . (23) Q eff is not a function of x (cid:48) as a manifestation of the local biaslimit, i.e. the pair separation scale is much smaller than theseparation to the third galaxy so that Q eff ’s dependence on x (cid:48) converges to scale with ξ ( x (cid:48) ) . Also note that when pairseparation is extremely small, i.e. | d p | ∼ Mpc, we have ξ ( d p ) (cid:29) , and Q eff would tend to b p , g / .However, as we will see in Section 3, the pairs we areconsidering have a maximum separation of up to 10 or 30Mpc, where the | x (cid:48) | (cid:29) | d p | limit does not always hold. Specifically, we can no longer assume ξ ( x (cid:48) ) ≈ ξ ( x (cid:48) ) ≈ ξ ( x (cid:48) ) and also the ξ ( x (cid:48) ) ξ ( x (cid:48) ) term contributes non-trivially toequation 21. We use the general form of equation 20 ζ ( d p , x (cid:48) , x (cid:48) ) = Q eff (cid:2) ξ ( x (cid:48) ) ξ ( d p ) + ξ ( x (cid:48) ) ξ ( d p ) + ξ ( x (cid:48) ) ξ ( x (cid:48) ) (cid:3) . (24)Assuming ξ ( r ) ∝ r − and treating d p / x (cid:48) as a small num-ber, we Taylor expand the first two terms, ξ ( x (cid:48) ) ξ ( d p ) and ξ ( x (cid:48) ) ξ ( d p ) , to the fourth order O (cid:104)(cid:0) d p / x (cid:48) (cid:1) (cid:105) . We also pre-serve the zero-th order approximation of the Taylor expan-sion of the last term ξ ( x (cid:48) ) ξ ( x (cid:48) ) ≈ ξ ( x (cid:48) ) ∝ ( d p / x (cid:48) ) . TheTaylor expanded form of the first two terms still carries adependence on the angle between x (cid:48) and d p . As a roughestimate, we average over the angle between x (cid:48) and d p as-suming isotropy, we get ξ ( x (cid:48) ) + ξ ( x (cid:48) ) ≈ ξ ( x (cid:48) ) + ξ ( x (cid:48) ) ξ ( d p ) , (25)and ζ ( d p , x (cid:48) , x (cid:48) ) ≈ Q eff (cid:20) ξ ( x (cid:48) ) ξ ( d p ) + ξ ( x (cid:48) ) (cid:21) . (26)Note that since we are assuming isotropy and a simplifiedform for ξ ( r ) , the coefficient in front of ξ ( x (cid:48) ) is only anestimate of the correction. However, a correction to thecorrection term comes in as a higher order correction to ζ ( d p , x (cid:48) , x (cid:48) ) , so we ignore these high order effects for thisderivation. Then, in the not strictly squeezed limit, the pair-galaxy bias is given by b p , g = ( + Q eff ξ ( d p )) + ξ ( d p ) + δ b . (27)where δ b = ξ ( d p ) + Q eff + ξ ( d p ) ∫ d x (cid:48) W ( x (cid:48) ) ξ ( x (cid:48) ) ∫ d x (cid:48) W ( x (cid:48) ) ξ ( x (cid:48) ) (28)Given a weighting function and a 2PCF, the fractionbetween the two integrals yields a constant. For ease ofnotation, we express this term as the 2PCF evaluated atsome characteristic separation d w , which is dependent onthe weighting function and the shape of the 2PCF. Specifi-cally, we define ξ ( d w ) ≡ ∫ d x (cid:48) W ( x (cid:48) ) ξ ( x (cid:48) ) ∫ d x (cid:48) W ( x (cid:48) ) ξ ( x (cid:48) ) . (29)To calculate the characteristic separation, we employ thepower-law form of the galaxy-galaxy 2PCF, ξ ( r ) ∝ r − p . Thenequation 27 gives the following formula to calculate d w d w = (cid:34) ∫ d x (cid:48) W ( x (cid:48) )| x (cid:48) | − p ∫ d x (cid:48) W ( x (cid:48) )| x (cid:48) | − p (cid:35) / p . (30)Having introduced ξ ( d w ) , the pair-galaxy bias can then beconveniently expressed with b p , g = + ξ ( d w ) ξ ( d p ) + Q eff ( ξ ( d p ) + ξ ( d w )) + ξ ( d p ) , (31)which translates to a modified estimator of the reduced3PCF Q eff = ( + ξ ( d p )) b p , g − − ξ ( d w ) ξ ( d p ) ξ ( d p ) + ξ ( d w ) . (32) MNRAS000
Then we have n p (cid:213) i ( W ∗ δ g )( x p , i ) = ∫ d x (cid:48) W ( x (cid:48) )( ξ ( x (cid:48) ) + ξ ( x (cid:48) )) + ξ ( d p ) + ∫ d x (cid:48) W ( x (cid:48) ) ζ ( d p , x (cid:48) , x (cid:48) ) + ξ ( d p ) . (18)Then, using equation 10, equation 11, and equation 18, wecan express the pair-galaxy bias in terms of the 2PCF andthe 3PCF b p , g = A + B + ξ ( d p ) , (19)where A = ∫ d x (cid:48) W ( x (cid:48) )( ξ ( x (cid:48) ) + ξ ( x (cid:48) )) ∫ d x (cid:48) W ( x (cid:48) ) ξ ( x (cid:48) ) , B = ∫ d x (cid:48) W ( x (cid:48) ) ζ ( d p , x (cid:48) , x (cid:48) ) ∫ d x (cid:48) W ( x (cid:48) ) ξ ( x (cid:48) ) . So now we have expressed the pair-galaxy bias as given byequation 11 in terms of general case 2PCFs and 3PCFs.To proceed, we introduce the reduced 3PCF Q eff , de-fined by Q eff ( x , x , x ) = ζ ( x , x , x ) ξ ( x ) ξ ( x ) + ξ ( x ) ξ ( x ) + ξ ( x ) ξ ( x ) , (20)where x , x , x are the three separation vectors for any3PCF. In the original hierarchical model (Peebles & Groth1975; Soneira & Peebles 1976; Groth & Peebles 1977; Peebles1980), the Q eff factor is a constant, which has been disprovenin later surveys and simulations (e.g. Kulkarni et al. 2007;Mar´ın et al. 2008; McBride et al. 2011). As we will see, ourresults in the squeezed limit also show a strongly variable Q eff .In the strictly squeezed limit, | x (cid:48) | (cid:29) | d p | and ξ ( x (cid:48) ) ≈ ξ ( x (cid:48) ) ≈ ξ ( x (cid:48) ) . A typical 2PCF assumes a power-law form ξ ( x (cid:48) ) ∝ x (cid:48) α , where the exponent α ≈ − . This means that, inthe strictly squeezed limit, ξ ( x (cid:48) ) (cid:28) ξ ( d p ) , and thus ζ ( d p , x (cid:48) , x (cid:48) ) ≈ Q eff ( ξ ( x (cid:48) ) ξ ( d p ) + ξ ( x (cid:48) ) )≈ Q eff ξ ( x (cid:48) ) ξ ( d p ) . (21)Then combining equation 19 and equation 21, we expresspair-galaxy bias as b p , g = ( + Q eff ξ ( d p )) + ξ ( d p ) . (22)Then we can express the reduced 3PCF Q eff as a functionof pair-galaxy bias b p , g Q eff = ( + ξ ( d p )) b p , g − ξ ( d p ) . (23) Q eff is not a function of x (cid:48) as a manifestation of the local biaslimit, i.e. the pair separation scale is much smaller than theseparation to the third galaxy so that Q eff ’s dependence on x (cid:48) converges to scale with ξ ( x (cid:48) ) . Also note that when pairseparation is extremely small, i.e. | d p | ∼ Mpc, we have ξ ( d p ) (cid:29) , and Q eff would tend to b p , g / .However, as we will see in Section 3, the pairs we areconsidering have a maximum separation of up to 10 or 30Mpc, where the | x (cid:48) | (cid:29) | d p | limit does not always hold. Specifically, we can no longer assume ξ ( x (cid:48) ) ≈ ξ ( x (cid:48) ) ≈ ξ ( x (cid:48) ) and also the ξ ( x (cid:48) ) ξ ( x (cid:48) ) term contributes non-trivially toequation 21. We use the general form of equation 20 ζ ( d p , x (cid:48) , x (cid:48) ) = Q eff (cid:2) ξ ( x (cid:48) ) ξ ( d p ) + ξ ( x (cid:48) ) ξ ( d p ) + ξ ( x (cid:48) ) ξ ( x (cid:48) ) (cid:3) . (24)Assuming ξ ( r ) ∝ r − and treating d p / x (cid:48) as a small num-ber, we Taylor expand the first two terms, ξ ( x (cid:48) ) ξ ( d p ) and ξ ( x (cid:48) ) ξ ( d p ) , to the fourth order O (cid:104)(cid:0) d p / x (cid:48) (cid:1) (cid:105) . We also pre-serve the zero-th order approximation of the Taylor expan-sion of the last term ξ ( x (cid:48) ) ξ ( x (cid:48) ) ≈ ξ ( x (cid:48) ) ∝ ( d p / x (cid:48) ) . TheTaylor expanded form of the first two terms still carries adependence on the angle between x (cid:48) and d p . As a roughestimate, we average over the angle between x (cid:48) and d p as-suming isotropy, we get ξ ( x (cid:48) ) + ξ ( x (cid:48) ) ≈ ξ ( x (cid:48) ) + ξ ( x (cid:48) ) ξ ( d p ) , (25)and ζ ( d p , x (cid:48) , x (cid:48) ) ≈ Q eff (cid:20) ξ ( x (cid:48) ) ξ ( d p ) + ξ ( x (cid:48) ) (cid:21) . (26)Note that since we are assuming isotropy and a simplifiedform for ξ ( r ) , the coefficient in front of ξ ( x (cid:48) ) is only anestimate of the correction. However, a correction to thecorrection term comes in as a higher order correction to ζ ( d p , x (cid:48) , x (cid:48) ) , so we ignore these high order effects for thisderivation. Then, in the not strictly squeezed limit, the pair-galaxy bias is given by b p , g = ( + Q eff ξ ( d p )) + ξ ( d p ) + δ b . (27)where δ b = ξ ( d p ) + Q eff + ξ ( d p ) ∫ d x (cid:48) W ( x (cid:48) ) ξ ( x (cid:48) ) ∫ d x (cid:48) W ( x (cid:48) ) ξ ( x (cid:48) ) (28)Given a weighting function and a 2PCF, the fractionbetween the two integrals yields a constant. For ease ofnotation, we express this term as the 2PCF evaluated atsome characteristic separation d w , which is dependent onthe weighting function and the shape of the 2PCF. Specifi-cally, we define ξ ( d w ) ≡ ∫ d x (cid:48) W ( x (cid:48) ) ξ ( x (cid:48) ) ∫ d x (cid:48) W ( x (cid:48) ) ξ ( x (cid:48) ) . (29)To calculate the characteristic separation, we employ thepower-law form of the galaxy-galaxy 2PCF, ξ ( r ) ∝ r − p . Thenequation 27 gives the following formula to calculate d w d w = (cid:34) ∫ d x (cid:48) W ( x (cid:48) )| x (cid:48) | − p ∫ d x (cid:48) W ( x (cid:48) )| x (cid:48) | − p (cid:35) / p . (30)Having introduced ξ ( d w ) , the pair-galaxy bias can then beconveniently expressed with b p , g = + ξ ( d w ) ξ ( d p ) + Q eff ( ξ ( d p ) + ξ ( d w )) + ξ ( d p ) , (31)which translates to a modified estimator of the reduced3PCF Q eff = ( + ξ ( d p )) b p , g − − ξ ( d w ) ξ ( d p ) ξ ( d p ) + ξ ( d w ) . (32) MNRAS000 , 1–16 (2017)
PCF in squeezed limit Thus the reduced 3PCF estimator only differs from thestrictly squeezed limit by an extra factor of ξ ( d w ) ξ ( d p ) on thenumerator and a factor of ξ ( d w ) on the denominator. Forthe rest of this paper, since we mostly do not operate in thestrictly squeezed limit, we use equation 32 to calculate Q eff . In this section, we present the implementation of our b p , g and Q eff statistic for Luminous Red Galaxies (LRGs) using aseries of N-body simulation boxes. LRGs are highly luminousearly-type massive galaxies consisting of mainly old starsand quenched of star formation. These luminous galaxies arepopular targets for large-scale structure studies and surveys.We study the statistics of LRGs specifically in this paper inpreparation for applying our methods to Baryon OscillationSpectroscopic Survey (BOSS; Eisenstein et al. 2011; Dawsonet al. 2013) data in follow up papers.We also present some basic characteristics of the simula-tion results in this section, namely the 2PCF, pair separationdistribution, and halo mass distribution of the LRGs andtheir pairs. It should be pointed out that we are only inter-ested in showcasing our methodology in this paper, and wedo not attempt any detailed matching with observed LRGstatistics. For the purpose of this analysis, we utilize the cosmologicalsimulations (Garrison et al., in preparation) generated bythe fast and high precision
Abacus
N-body code (Garrisonet al. 2016, Ferrer et al., in preparation; Metchnik & Pinto, inpreparation). Specifically we use a series of 16 cyclic boxeswith Planck 2016 cosmology (Planck Collaboration et al.2016) at redshift z = . , where each box is of size 1100 h − Mpc, and contains 1440 dark matter particles of mass × h − M (cid:12) . Dark matter halos are found and charac-terized using the friends-of-friends (FoF; Davis et al. 1985)halo finder. Different halo finders, specifically ROCKSTAR(Behroozi et al. 2013), behave similarly for halos of at least afew hundred particles (Garrison et al. 2016). For the rest ofthis paper, we use halo catalogues generated using the FoFhalo finder with linking length . b , where b is the averageinterparticle spacing. We have also repeated the key analysiswith mock catalogues based on halos from ROCKSTAR andrecovered consistent results.Although FoF finder is fast, one main concern with FoFfinder is the effect of FoF overbridging (Lukic 2008; Knebeet al. 2011), which misidentifies a small fraction of physicallyseparate halos as a single halo through long chain of links.This effect occurs at all scales and eventually leads to asmall over-estimate of the number of galaxy pairs inhabitingthe same halo. For our purpose, the correction is small andwould not impact the qualitative results of this paper. An important effect to consider is Redshift Space Distortion(RSD, Hamilton 1998). Essentially, the peculiar velocity ofthe galaxies in clusters along the line-of-sight (LOS) leads to a deviation in distance measurements when using Hubble’slaw. This is commonly known as the “finger-of-god” (FoG)effect. There is also the more subtle effect known as theKaiser effect (Kaiser 1987), where the coherent infall velocityof galaxies towards central masses introduces a flattening inthe observed structure on the large scale. It is importantto include RSD in our calculations to give results that arecomparable to observables. To include RSD, we shift the x (cid:107) coordinate of each galaxy by x (cid:48)(cid:107) = x (cid:107) + v (cid:107) H ( z ) , (33)where v (cid:107) is the LOS velocity of the galaxy. For simplicity,we set the LOS velocity of the galaxy equal to the LOS ve-locity of the host dark matter particle, ignoring the ∼ relative velocity bias due to the shear motion between thegalaxies and their underlying halos (Guo et al. 2015a). Thenwe can derive the pair-galaxy bias and the squeezed 3PCFin redshift space using the distorted positions of the galax-ies. Several previous papers have also studied the behav-iors of redshift space 3PCF in the context of galaxy surveysand simulations (e.g. Gazta˜naga & Scoccimarro 2005; Mar´ınet al. 2008; Guo et al. 2015b).For the rest of this paper, we calculate the squeezed3PCF in both real space and redshift space. To avoid con-fusion, we explicit state in the figure caption whether thefigures are made in real space or redshift space. The first step is to generate mock LRG catalogue from thehalo catalogue using a HOD prescription. The HOD we useis a slightly modified version of the the simple HOD fitted tothe LRG 2PCF in the Sloan Digital Sky Survey (SDSS; Yorket al. 2000; Zehavi et al. 2002; Abazajian et al. 2003; Zehaviet al. 2004; Eisenstein et al. 2005a; Zehavi et al. 2005a) andsummarized in Zheng et al. (2009); Kwan et al. (2015). Theform of our fiducial HOD (equation 34; Kravtsov et al. 2004;Zheng et al. 2005, 2007) is simple in the sense that it onlyhas five parameters, only depends on the halo mass, and ig-nores effects of halo assembly and halo environment. A briefdiscussion on more sophisticated halo occupation models isdeferred to Section 5.2, but for the purpose of providingan simple application of our methodology, this simple HODthat solely depends on halo mass is sufficient (e.g. Kauff-mann et al. 2004; Mo et al. 2004; Blanton et al. 2006; Tinkeret al. 2008).Specifically, this HOD prescription gives the number ofcentrals and satellites as (cid:104) n cent (cid:105) =
12 erfc (cid:20) ln ( M cut / M )√ σ (cid:21) , (cid:104) n sat (cid:105) = (cid:18) M − κ M cut M (cid:19) α , (34)where the parameters are chosen as, M cut ≈ . M (cid:12) , M ≈ . M (cid:12) , σ = . , α = , and κ = (Zheng et al. 2009;Kwan et al. 2015). The halo mass threshold for satellites isgiven by M min = κ M cut . This HOD is plotted in Figure 1.From here on, we refer to this HOD as the Z09 design. Alsonote that we truncate the HOD at M = × M (cid:12) to dis-regard halos with fewer than 100 particles because smallerhalos are unreliable and dependent on which halo finder we MNRAS , 1–16 (2017)
S. Yuan, D. J. Eisenstein and L. H. Garrison
Figure 1.
The HOD plotted as number of LRGs per halo for agiven halo mass M (Zheng et al. 2009; Kwan et al. 2015). Theblue dashed line is the number of centrals whereas the red solidline is the number of satellites. The total number of LRGs perhalo is given by the black solid line. use. The central is assigned to the centre of mass of the halowith probability (cid:104) n cen (cid:105) . The satellites are assigned to eachparticle within the halo with probability p = (cid:104) n sat (cid:105)/ n particle ,where n particle is the total number of particles within the halo.We assign satellite galaxies at particle positions instead ofby fitting a smooth profile (e.g. Navarro-Frenk-White profile;Navarro et al. 1996, 1997) since we want to allow satellitesto trace halo substructure and avoid assumptions about theisotropy or equilibrium state of the halos.Looping through all the dark matter particles and as-signing galaxies with a small probability introduces signif-icant shot noise, which leads to uncertainties in our b p , g and Q eff signals not greater than (refer to Figure 8 andFigure 9). For most of our analysis, the b p , g and Q eff sig-nals are strong enough that this shot noise is not important.However, in some cases we require much higher signal-to-noise, such as in Section 5, where we investigate percentlevel perturbations to the b p , g and Q eff signal due to vari-ations in HOD parameters. To suppress shot noise, we seedthe Python pseudo random number generator with the samenumber at the start of generating the mock galaxy cata-logue for each simulation. Thus, when we loop through allthe dark matter particles to populate them with galaxies,each dark matter particle is called with the same pseudorandom number on each run. This way, not only do we getreproducible results when we re-run the same HOD, but alsosmall changes in the HOD make small changes in the galaxydistribution, thereby sharply reducing the noise in the nu-merical derivatives of the correlation functions with respectto HOD parameters.To further suppress the shot noise, we populate eachof the 16 simulation boxes with the same HOD 16 differenttimes, each time with a different seed. Then we take theaverage Q eff of the realizations. Note that the reseeding isonly done for Figure 10 and Figure 13, due to computationallimit.Once we have compiled the mock LRG catalogue fromour HOD, the LRG pairs are generated with a pair search on Figure 2.
The real-space pair separation distribution up to d p , max = Mpc for one of the 16 simulation boxes. The red his-togram represents the 1-halo pairs, whereas the blue histogramrepresents the 2-halo pairs, and the black histogram is the sumof both. a kd-tree of LRG locations. For our simulation boxes, thisHOD typically generates around . × LRGs per box. Formaximum pair separation d p , max = Mpc in real space, wegenerate around . × pairs per box, of which around . × are identified as 1-halo pairs, i.e. the two LRGsforming the pair occupy the same halo. We clarify that allpairs within maximum pair separation d p , max are counted inour pair catalogues, e.g. we include all 3 pairs generated bya close triplet. Also for the purpose of this paper, we do notconsider effects of fiber collisions or survey boundaries.Figure 2 shows the distribution of pair separation in all16 simulation boxes in real space. We see that at d p < Mpc,1-halo pairs dominate, while at d p > Mpc, 2-halo pairsdominate.In Figure 3, we plot the halo mass corresponding to theLRG pair as a function of pair separation d p in real space.For 1-halo pairs, the halo mass M is simply the mass of thehost halo, and for 2-halo pairs, M is the sum of the massof the two halos. The width of the shaded region shown onthe plot is twice the sample standard deviation at that d p .We see a distinct peak at around d p ≈ Mpc, suggestingthat close LRG pairs at approximately 2 Mpc apart samplethe most massive halos, whereas at separation less or greaterthan d p ∼ Mpc, the halo mass of the pair decreases.Figure 4 shows the halo mass probability density distri-bution for galaxies (in red) and for pairs (in blue), with theaverage halo mass denoted by the dashed curves. The pairsample shown on this figure has a maximum pair separa-tion d p , max = Mpc in real space. We see that the galaxiessample halos of approximately . M (cid:12) whereas the pairssample more massive halos at around . M (cid:12) . This plotmotivates the point in Section 1 that close pairs are biasedtracers of more massive halos and that statistics on the pairfield can be used to constrain the high mass regime of theHOD. MNRAS000
The real-space pair separation distribution up to d p , max = Mpc for one of the 16 simulation boxes. The red his-togram represents the 1-halo pairs, whereas the blue histogramrepresents the 2-halo pairs, and the black histogram is the sumof both. a kd-tree of LRG locations. For our simulation boxes, thisHOD typically generates around . × LRGs per box. Formaximum pair separation d p , max = Mpc in real space, wegenerate around . × pairs per box, of which around . × are identified as 1-halo pairs, i.e. the two LRGsforming the pair occupy the same halo. We clarify that allpairs within maximum pair separation d p , max are counted inour pair catalogues, e.g. we include all 3 pairs generated bya close triplet. Also for the purpose of this paper, we do notconsider effects of fiber collisions or survey boundaries.Figure 2 shows the distribution of pair separation in all16 simulation boxes in real space. We see that at d p < Mpc,1-halo pairs dominate, while at d p > Mpc, 2-halo pairsdominate.In Figure 3, we plot the halo mass corresponding to theLRG pair as a function of pair separation d p in real space.For 1-halo pairs, the halo mass M is simply the mass of thehost halo, and for 2-halo pairs, M is the sum of the massof the two halos. The width of the shaded region shown onthe plot is twice the sample standard deviation at that d p .We see a distinct peak at around d p ≈ Mpc, suggestingthat close LRG pairs at approximately 2 Mpc apart samplethe most massive halos, whereas at separation less or greaterthan d p ∼ Mpc, the halo mass of the pair decreases.Figure 4 shows the halo mass probability density distri-bution for galaxies (in red) and for pairs (in blue), with theaverage halo mass denoted by the dashed curves. The pairsample shown on this figure has a maximum pair separa-tion d p , max = Mpc in real space. We see that the galaxiessample halos of approximately . M (cid:12) whereas the pairssample more massive halos at around . M (cid:12) . This plotmotivates the point in Section 1 that close pairs are biasedtracers of more massive halos and that statistics on the pairfield can be used to constrain the high mass regime of theHOD. MNRAS000 , 1–16 (2017)
PCF in squeezed limit Figure 3.
The halo mass M of the pairs as a function of real-space pair separation d p in real space. For 1-halo pairs, M issimply the mass of the host halo, while for 2-halo pairs, M is thesum of the mass of the two halos. The values plotted are averagedacross all 16 simulations. The linewidth shows the sample stan-dard deviation at that d p . We see that the halo mass stronglypeaks at around 2 Mpc, where we have mostly 1-halo pairs butalso a significant fraction of 2-halo pairs. Figure 4.
The halo mass M probability density distribution forgalaxies (in red) and for pairs (in blue). The dashed curves show-case the average halo mass of the galaxies and the pairs. Thepair sample used for this plot has a maximum pair separation d p , max = Mpc in real space. We see that the pairs at this sep-aration trace more massive halos than galaxies alone.
By the pair-counting definition of 2PCF, we can calculatethe 2PCF as ξ ( d p ) = N p ( d p ) ¯ n g N g ∆ V shell − , (35)where N p ( d p ) is the number of pairs with separation d p and ∆ V shell = π d p ∆ d p is the shell volume. ¯ n g denotes the LRGnumber density, and N g denotes the total number of LRGsin the sample. Figure 5.
The real space 2PCF as a function of pair separa-tion generated from the simulation box is in red solid line. The2PCFs from SDSS CMASS sample in two different redshift binsare plotted in dashed blue line and dashed green line respectively.
If we consider RSD, the LOS separation of the pair d (cid:107) becomes largely meaningless due to velocity dispersionalong the LOS. Thus, we calculate the LOS projected 2PCF, w ( d ⊥ ) , as a function of the perpendicular (to LOS) compo-nent of pair separation d ⊥ , w ( d ⊥ ) = ∫ d (cid:107) , max − d (cid:107) , max ξ ( d (cid:107) , d ⊥ ) d ( d (cid:107) ) = d (cid:107) , max (cid:18) σ p ( d ⊥ ) ¯ n g N g ∆ V ⊥ − (cid:19) , (36)where d (cid:107) , max = (cid:113) d p , max − d ⊥ is the maximum pair sepa-ration along the LOS given a maximum pair separation d p , max . σ p ( d ⊥ ) is the integrated column pair density, and ∆ V ⊥ = π ∆ d ⊥ × d (cid:107) , max . Unlike the real space 2PCF givenby equation 35, the LOS projected 2PCF can be inferredfrom observed LRG distributions without redshift informa-tion. For this reason, we use projected 2PCF instead of thefull 2PCF in Section 5.1 as well.Figure 5 shows w ( d ⊥ ) calculated from our LRG paircatalogue. We also plot in dashed blue and green curves theprojected 2PCFs from the CMASS luminous red galaxiessample of Sloan Digital Sky Survey III (SDSS-III; Eisen-stein et al. 2011; Dawson et al. 2013) within the redshiftrange . < z < . and . < z < . , respectively. These2PCFs are directly adapted from Guo et al. (2013). Our2PCF is roughly consistent with the observed 2PCFs fromCMASS red galaxies. The inconsistency in shape can be at-tributed to a combination of differences in cosmology, prob-lems with halo finders, and oversimplifications in our HODmodel itself. However, for the purpose of this methodologypaper, the exact form of the HOD is not essential, thus adetailed fit to the observed 2PCF is not important for ouranalysis. Applications of our methodology to more sophisti-cated halo occupation models are deferred to future papers. MNRAS , 1–16 (2017)
S. Yuan, D. J. Eisenstein and L. H. Garrison
To be able to utilize FFTs to efficiently compute the powerspectra of the LRG and LRG pair fields, we require a massassignment function that assigns LRGs (and pairs) onto aregular 3D grid in such a way that generates a smoothLRG (and pair) density field. For the purpose of this pa-per, we assign LRGs and pairs onto an evenly spaced grid spanning the full simulation box using the TriangularShaped Cloud (TSC) method (Jeong 2010). For a simula-tion box size L box = h − Mpc with 512 , the cell size is ∆ x ≈ . h − Mpc. Thus, a TSC of base length ∆ x spans27 cells on the grid, with 3 cells along each dimension. Thefact that each TSC can span at most 3 cells along each di-rection becomes important in the next section as we designour custom weighting functions. In Section 2, we present our method of efficiently comput-ing the cross-correlation between galaxy pairs and galaxies,from which we infer the squeezed 3PCF. This routine re-quires a robust weighting function for the pair-galaxy cross-correlation that satisfies the following three criteria.(a) The weighting function should prevent self-counting,where one of the two galaxies within the pair is consideredas the third galaxy when performing the cross-correlation.This means that the weighting function should have a hollowspherical cavity in the centre. The radius of the cavity shouldbe r cut ≈ d p , max / + ∆ x , where d p , max is the maximum pairseparation and ∆ x is the cell size. The extra ∆ x comes fromavoiding the TSC of the pair, which can extend a maximumof 2 cells outwards from either galaxy along the axis of thepair.(b) The weighting function should minimize the FoG ef-fect when we apply our methodology to redshift space data,where the correlation along the LOS becomes largely mean-ingless. To disregard correlation between pairs and galaxiesalong the LOS, we “drill” a cylindrical hole through the cen-tre of the weighting function parallel to the LOS. The radiusof the cylindrical hole should be r cut ≈ d p , max / + ∆ x for thesame reason as explained in (a).(c) The weighting function should sample scales approx-imately in the range k = . ∼ . h Mpc − . The scalerange is limited on the large scale by the simulation boxsize, where we cannot meaningfully interpret scales largerthan the box size, which translates to a minimum frequencyof k min = π / L box ≈ . h Mpc − . The scale range is limitedon the small scale by self-counting, so we do not want thethird galaxy to be too close to the pair that their TSCs over-lap. We also avoid small scale contribution since we want tobe computing bias in the linear regime.In addition to these requirements, we also want theweighting function to be mostly smooth and analytic for easeof performing FFTs and to avoid high k modes in the weight-ing function. Specifically, our weighting function is given inreal space by, W ( r ⊥ , r (cid:107) ) = N + exp [− s ( r ⊥ − r s )] exp (cid:20) − (cid:18) r ⊥ σ + r (cid:107) σ (cid:19)(cid:21) , r ⊥ ≥ r cut , r ⊥ < r cut (37) Figure 6.
The weighting functions for d p , max = Mpc (redsolid line) and d p , max = Mpc (blue dashed line) as defined inequation 37. We only show a 2D cutaway of the 3D function at r (cid:107) = . The dotted black lines mark r cut = Mpc and r cut = Mpc. where N is a normalization constant. r ⊥ and r (cid:107) are distancebetween the third galaxy and the centre of the pair, pro-jected perpendicular and parallel to the LOS respectively.The function is parameterized by σ , r s , s and r cut , whichare determined by the maximum pair separtion d p , max thatwe choose. Although the function is divided into two seg-ment and thus not fully analytical, such engineering is nec-essary to remove contamination from pair self counting. Wechoose the parameters such that the jump in W ( r ⊥ , r (cid:107) ) at r cut is small (refer to Figure 6).For the rest of the analysis in this paper, we use d p , max = Mpc when working in real space and d p , max = Mpcwhen working in redshift space. We choose d p , max = Mpcin redshift space because the typical velocity dispersionalong the LOS translates to ∼ Mpc displacement, so weextend the pair search in both direction by Mpc to cap-ture most of the pairs. The parameters of the weighting func-tion for both cases are summarized in Table 1.Figure 6 shows a 2D cutaway of the weighting functionsat r (cid:107) = in the positive half of r ⊥ . The solid red curve showsthe weighting function for d p , max = Mpc and the dashedblue curve shows the weighting function for d p , max = Mpc.The dotted black lines denote the location of r cut to showthat the jump at r cut is small.According to equation 11, the weighting function is mul-tiplied with the pair density field in Fourier space as a sam-pling weight for different k modes. The Fourier Transformsof the weighting functions are plotted in Figure 7, where thetop panel shows the weighting function for d p , max = Mpcwhereas the middle panel shows the weighting function for d p , max = Mpc. The bottom panel shows cutaways of thetwo cases at k (cid:107) = . The colorbar gives ˜ W k ⊥ / k , which isproportional to the number of modes sampled by the weight-ing function in Fourier space. We see that both weight-ing functions are indeed heavily sampling scales between . ∼ . h Mpc − . The first and second order ringing seenat higher frequency are due to the steep drop-off in the MNRAS000
The weighting functions for d p , max = Mpc (redsolid line) and d p , max = Mpc (blue dashed line) as defined inequation 37. We only show a 2D cutaway of the 3D function at r (cid:107) = . The dotted black lines mark r cut = Mpc and r cut = Mpc. where N is a normalization constant. r ⊥ and r (cid:107) are distancebetween the third galaxy and the centre of the pair, pro-jected perpendicular and parallel to the LOS respectively.The function is parameterized by σ , r s , s and r cut , whichare determined by the maximum pair separtion d p , max thatwe choose. Although the function is divided into two seg-ment and thus not fully analytical, such engineering is nec-essary to remove contamination from pair self counting. Wechoose the parameters such that the jump in W ( r ⊥ , r (cid:107) ) at r cut is small (refer to Figure 6).For the rest of the analysis in this paper, we use d p , max = Mpc when working in real space and d p , max = Mpcwhen working in redshift space. We choose d p , max = Mpcin redshift space because the typical velocity dispersionalong the LOS translates to ∼ Mpc displacement, so weextend the pair search in both direction by Mpc to cap-ture most of the pairs. The parameters of the weighting func-tion for both cases are summarized in Table 1.Figure 6 shows a 2D cutaway of the weighting functionsat r (cid:107) = in the positive half of r ⊥ . The solid red curve showsthe weighting function for d p , max = Mpc and the dashedblue curve shows the weighting function for d p , max = Mpc.The dotted black lines denote the location of r cut to showthat the jump at r cut is small.According to equation 11, the weighting function is mul-tiplied with the pair density field in Fourier space as a sam-pling weight for different k modes. The Fourier Transformsof the weighting functions are plotted in Figure 7, where thetop panel shows the weighting function for d p , max = Mpcwhereas the middle panel shows the weighting function for d p , max = Mpc. The bottom panel shows cutaways of thetwo cases at k (cid:107) = . The colorbar gives ˜ W k ⊥ / k , which isproportional to the number of modes sampled by the weight-ing function in Fourier space. We see that both weight-ing functions are indeed heavily sampling scales between . ∼ . h Mpc − . The first and second order ringing seenat higher frequency are due to the steep drop-off in the MNRAS000 , 1–16 (2017)
PCF in squeezed limit d p , max (Mpc) σ (Mpc) r s (Mpc) s (Mpc − ) r cut (Mpc) d w (Mpc) ξ ( d w )
10 30 20 0.8 10 30.9 0.36230 50 35 0.5 20 52.6 0.104
Table 1.
Parameters of the weighting functions for d p , max = Mpc and d p , max = Mpc. Again the cutoff is given by r cut ≈ d p , max / + ∆ x . Figure 7.
The Fourier Transforms of the weighting function for d p , max = Mpc (top panel) and d p , max = Mpc (middle panel).The bottom panel shows cutaways of ˜ W ( k ) k ⊥ k at k (cid:107) = for the d p , max = Mpc case (red solid curve) and the d p , max = Mpccase (blue dashed curve). real space weighting function towards small r ⊥ . However,the ringings have much lower integrated weight in k -spaceand are not expected to affect our analysis significantly.Now that we have defined the weighting functions, wecan calculate d w (equation 30) and ξ ( d w ) (equation 29),which we need to compute the squeezed 3PCF coefficient Q eff from the pair-galaxy bias b p , g . To do so, we employ the2PCF given by Figure 2 and Figure 3 in Eisenstein et al.(2005b), which is constructed from the linear power spec-tra calculated using CMBFAST (Seljak & Zaldarriaga 1996;Zaldarriaga et al. 1998; Zaldarriaga & Seljak 2000). For anyof our weighting functions, we fit the 2PCF with a power lawin the scale range sampled by the weighting function. Thenwe use the best fit exponent p to calculate d w and ξ ( d w ) .Specifically, the d w and ξ ( d w ) for our d p , max = Mpc and30 Mpc are given in the last two columns of Table 1. Itshould be pointed out that the analysis in Eisenstein et al.(2005b) uses a slightly different cosmology, where h = . , Ω m h = . and Ω b h = . . However, the correction issmall and does not meaningfully change our calculations.Now that we have defined our mass assignment func-tion and our weighting functions, we apply the TSC to theLRGs and LRG pairs in all 16 simulation boxes to generatesmooth density fields of the LRGs and pairs on our grid. Theresulting density fields are then convolved with the appropri-ate weighting function to compute the pair-galaxy bias b p , g (equation 11) and the reduced 3PCF Q eff (equation 32). In this section, we present the results of implementing ourefficient method of calculating b p , g and Q eff using the N-body simulation and the Z09 HOD. The results shown areaveraged over the 16 cosmological boxes to reduce noise.Section 4.1 presents the results in real space, and Section 4.2presents the results in redshift space. We average the b p , g calculated from all 16 simulation boxesat z = . and plot the averaged real space b p , g in the toppanel of Figure 8. We have chosen maximum pair separation d p , max = Mpc and a Gauss-Sigmoid weighting functionwith the following parameters: σ = Mpc, r s = Mpc, s = . Mpc − and r cut = Mpc (Refer to Section 3.6 fordetails). The associated uncertainty map is given in middlepanel. The bottom panel shows the reduced 3PCF Q eff inreal space, as calculated from b p , g using equation 32.The top panel of Figure 8 shows that the pair-galaxybias is close to 1.3 at pair separation less than ∼ Mpc, butdramatically increases by a factor of 2 and peaks at ∼ Mpc,before decreasing again towards greater pair separation. Thebottom panel shows that the squeezed 3PCF coefficient Q eff recovers the same behavior. We can analytically estimate thevalue of b p , g by computing the halo bias as a function of halomass and then take the ratio of the halo bias averaged overall pairs and the halo bias averaged over all galaxies (usingthe halo mass function of the pairs and galaxies shown inFigure 4). Our estimate yields b p , g ≈ . , which is consistentwith the values we see in the top panel of Figure 8.Referring to Figure 2 and Figure 3, we see that pairswith separation < Mpc are mostly 1-halo pairs with lowerhalo mass. b p , g is approximately 1.3 at d p < Mpc. The factthat the pair-galaxy bias b p , g is greater than 1 shows thatclose galaxy pairs are more biased than individual galaxies.Pairs with separation ∼ Mpc consist of ∼ larger 1-halopairs and ∼ b p , g and Q eff at large pair separation d p > Mpc. The anisotropy
MNRAS , 1–16 (2017) S. Yuan, D. J. Eisenstein and L. H. Garrison
Figure 8. b p , g and Q eff results in real space. The top panelplots the pair-galaxy bias distribution as a function of real-spacepair separation. The middle panel plots the corresponding un-certainty of the pair-galaxy bias distribution. The uncertainty iscalculated as the standard deviation of the pair-galaxy bias cal-culated for all 16 simulation boxes. The bottom panel plots thecorresponding reduced 3PCF Q eff . The weighting function is theGaussian-sigmoid function for d p , max = Mpc as described inSection. 3.6. is due to the fact that the weighting function samples favor-ably volume in the transverse direction and avoids volumesalong the LOS direction. When a pair is oriented transverseto the LOS, it establishes the transverse direction as thepreferred direction for the large-scale filament. Thus, whena pair is perpendicular to the LOS, the weighting function overlaps the most with the large-scale filament, giving riseto much stronger clustering. Conversely, when a pair is ori-ented along the direction of LOS, the filament is preferen-tially aligned with the LOS, which does not overlap signifi-cantly with the weighting function, thus giving rise to weakerclustering. The large b p , g value at d p > Mpc suggests thatthe filament, not just the halo, is also strongly biased againstthe galaxy field.The middle panel of Figure 8 shows that most of the b p , g values have percent level uncertainties, except the onestowards the top left corner of the plot, where the small sam-ple size contributes to the high noise.The bottom panel of Figure 8 shows that the reduced3PCF recovers the same behavior as seen in the pair-galaxybias b p , g . This is consistent with our conclusion that in thestrictly squeezed limit, the 3PCF would basically reduce to Q eff = b p , g / . At pair separation d p > Mpc, we continue tosee the anisotropy in bias due to the filamentary structurealigning with the preferred direction of the weighting func-tion. However, we also see suppression of a few percent in Q eff due to the ξ ( d w ) term. The associated uncertainty mapof Q eff is omitted, but shares the same properties as theuncertainty in b p , g , but with approximately half the ampli-tude. Similar to the last section, we can incorporate RSD andcalculate the b p , g and Q eff in redshift space, as shown inFigure 9. The top panel shows the the pair-galaxy bias b p , g as a function of pair separation. The middle panel shows thecorresponding uncertainty on b p , g . The bottom panel showsthe reduced 3PCF Q eff in redshift space. For these plots, theweighting function is chosen to be the d p , max = Mpc. Fig-ure 10 shows a zoom-in version of the Q eff signal in redshiftspace with finer bins on pair separation to be comparableto the real space Q eff signal plotted in the bottom panelof Figure 8. To compensate for the loss in signal-to-noiseby choosing finer bins, we use 16 different seeds on top of16 different simulation boxes, as discussed in Section 3.3.One should disregard the high values along the outer circu-lar edge of plots as these values are calculated from smallsample of pairs and are thus dominated by noise.The top and bottom panel of Figure 9 show a vertical“ridge” at d ⊥ ≈ Mpc in b p , g and Q eff that corresponds tothe ∼ Mpc peak in the real space b p , g and Q eff signal (Fig-ure 8). This ridge tracks the halo mass peak at d p ∼ Mpcin Figure 3, but smeared out along the LOS by FoG. Thestrong clustering along the ridge all the way to d (cid:107) ∼ Mpcsuggests that the galaxy pair-wise peculiar velocities reach ∼ km/s in the most massive halos. The zoom-in inFigure 10 shows that the ridge is more clustered at largerLOS separation, which is due to the fact that larger LOSpair separation corresponds to larger peculiar velocities ofthe galaxies, which in turn corresponds to more massive ha-los. We also point out that the infall velocity of mergingdark matter halos may also contribute to this ridge. How-ever, we do not expect the infall velocity to frequently reach ∼ km/s and speculate that the strongest Q eff signal to-wards the top of the ridge should be completely dominatedby FoG.We also see weaker anisotropy at large pair separation MNRAS000
Figure 8. b p , g and Q eff results in real space. The top panelplots the pair-galaxy bias distribution as a function of real-spacepair separation. The middle panel plots the corresponding un-certainty of the pair-galaxy bias distribution. The uncertainty iscalculated as the standard deviation of the pair-galaxy bias cal-culated for all 16 simulation boxes. The bottom panel plots thecorresponding reduced 3PCF Q eff . The weighting function is theGaussian-sigmoid function for d p , max = Mpc as described inSection. 3.6. is due to the fact that the weighting function samples favor-ably volume in the transverse direction and avoids volumesalong the LOS direction. When a pair is oriented transverseto the LOS, it establishes the transverse direction as thepreferred direction for the large-scale filament. Thus, whena pair is perpendicular to the LOS, the weighting function overlaps the most with the large-scale filament, giving riseto much stronger clustering. Conversely, when a pair is ori-ented along the direction of LOS, the filament is preferen-tially aligned with the LOS, which does not overlap signifi-cantly with the weighting function, thus giving rise to weakerclustering. The large b p , g value at d p > Mpc suggests thatthe filament, not just the halo, is also strongly biased againstthe galaxy field.The middle panel of Figure 8 shows that most of the b p , g values have percent level uncertainties, except the onestowards the top left corner of the plot, where the small sam-ple size contributes to the high noise.The bottom panel of Figure 8 shows that the reduced3PCF recovers the same behavior as seen in the pair-galaxybias b p , g . This is consistent with our conclusion that in thestrictly squeezed limit, the 3PCF would basically reduce to Q eff = b p , g / . At pair separation d p > Mpc, we continue tosee the anisotropy in bias due to the filamentary structurealigning with the preferred direction of the weighting func-tion. However, we also see suppression of a few percent in Q eff due to the ξ ( d w ) term. The associated uncertainty mapof Q eff is omitted, but shares the same properties as theuncertainty in b p , g , but with approximately half the ampli-tude. Similar to the last section, we can incorporate RSD andcalculate the b p , g and Q eff in redshift space, as shown inFigure 9. The top panel shows the the pair-galaxy bias b p , g as a function of pair separation. The middle panel shows thecorresponding uncertainty on b p , g . The bottom panel showsthe reduced 3PCF Q eff in redshift space. For these plots, theweighting function is chosen to be the d p , max = Mpc. Fig-ure 10 shows a zoom-in version of the Q eff signal in redshiftspace with finer bins on pair separation to be comparableto the real space Q eff signal plotted in the bottom panelof Figure 8. To compensate for the loss in signal-to-noiseby choosing finer bins, we use 16 different seeds on top of16 different simulation boxes, as discussed in Section 3.3.One should disregard the high values along the outer circu-lar edge of plots as these values are calculated from smallsample of pairs and are thus dominated by noise.The top and bottom panel of Figure 9 show a vertical“ridge” at d ⊥ ≈ Mpc in b p , g and Q eff that corresponds tothe ∼ Mpc peak in the real space b p , g and Q eff signal (Fig-ure 8). This ridge tracks the halo mass peak at d p ∼ Mpcin Figure 3, but smeared out along the LOS by FoG. Thestrong clustering along the ridge all the way to d (cid:107) ∼ Mpcsuggests that the galaxy pair-wise peculiar velocities reach ∼ km/s in the most massive halos. The zoom-in inFigure 10 shows that the ridge is more clustered at largerLOS separation, which is due to the fact that larger LOSpair separation corresponds to larger peculiar velocities ofthe galaxies, which in turn corresponds to more massive ha-los. We also point out that the infall velocity of mergingdark matter halos may also contribute to this ridge. How-ever, we do not expect the infall velocity to frequently reach ∼ km/s and speculate that the strongest Q eff signal to-wards the top of the ridge should be completely dominatedby FoG.We also see weaker anisotropy at large pair separation MNRAS000 , 1–16 (2017)
PCF in squeezed limit Figure 9. b p , g and Q eff results in redshift space. The top panelplots the pair-galaxy bias distribution as a function of redshift-space pair separation. The middle panel plots the correspondinguncertainty of the pair-galaxy bias distribution. The uncertaintyis calculated as the standard deviation of the pair-galaxy biascalculated for all 16 simulation boxes. The bottom panel plotsthe corresponding reduced 3PCF Q eff . The weighting function isthe Gaussian-sigmoid function for d p , max = Mpc as describedin Section. 3.6. in redshift space. This is due to the RSD smearing out thefilamentary structures, thus also randomizing the preferreddirection given an LRG pair. The rise in bias value towardslarger pair separation perpendicular to the LOS can be ex-plained by the fact that more distant pairs along the samefilament are better at picking out the preferred direction
Figure 10.
A zoom-in version of the Q eff signal in redshift space.Specifically, we use finer bins on pair separation and zoom in onthe d ⊥ axis so to be comparable to the real space Q eff signal plot-ted in the bottom panel of Figure 8. To compensate for the lossin signal-to-noise due to using finer bins, we re-use each of the 16different simulation boxes 16 times using different random seedsfor the HOD. The weighting function is the Gaussian-sigmoidfunction for d p , max = Mpc as described in Section. 3.6. of the filament, which gives rise to higher bias given theanisotropy of the weighting function.Both the real space and redshift space b p , g signal showa factor of 2 variation. This suggests that the hierarchicalansatz on the 3PCF does not hold in the squeezed limit.However, as we have discussed in this section, the variationsin the squeezed 3PCF signal provides rich information onthe galaxy-halo connection on the small scale and on thefilamentary structure on the large scale. In the followingsection, we explore using the squeezed 3PCF as a diagnosticson the galaxy-halo connection. MNRAS , 1–16 (2017) S. Yuan, D. J. Eisenstein and L. H. Garrison
We next discuss the relationship between the squeezed 3PCFand the HOD. Specifically, we explore the possibility of usingthe squeezed 3PCF to break degeneracies in the 2PCF andconstrain the high mass end of the HOD.
To investigate how the squeezed 3PCF breaks degeneraciesin the 2PCF, we first need to understand how the HOD con-trols the 2PCF and devise HOD recipes that would changethe 3PCF while maintaining the same 2PCF. First we at-tempt a linear emulator to model the 2PCF as a function ofHOD parameters. As discussed in Section 3.4, we study theprojected 2PCF w ( d ⊥ ) as given by equation 36. To the firstorder, we can expand the projected 2PCF as w i ( p ) ≈ w i ( p ) + (cid:213) j = w i , p j ( p j − p , j ) , (38)where p = [ p , p , p , p ] = (cid:2) log M cut , log M , σ, α (cid:3) are theHOD parameters that we use to control the 2PCF. Notethat we omitted parameter κ because our tests show that κ perturbs the 2PCF insignificantly for the variation ampli-tude that we are interested in. p denotes the Z09 design ofthe HOD, as quoted in Section 3.3. w i ( p ) denotes the valueof the projected 2PCF integrated over the i -th bin in d ⊥ .By choosing an uniform set of bins on the projected 2PCF,we form a vector [ w i ( p )] that holds the shape informationof the projected 2PCF. (cid:104) w i , p j (cid:105) records the matrix of thederivatives of the 2PCF integrated over the i -th bin, w i ( p ) ,against the j -th HOD parameter, p j . We can determine eachderivative w i , p j by independently varying each parameter ofthe HOD and tracking the variation incurred onto w i ( p ) .Once we have (cid:104) w i , p j (cid:105) , we solve for the set of HOD parame-ters p that minimizes | [ w i ( p )] − [ w i ( p )] | . However, it turnsout that the linear model does not produce a good emula-tor for the projected 2PCF. Our tests show that the linearapproximation breaks down when varying α by more than ∼ .Thus, we develop an emulator that incorporates second-order derivatives. To formulate a second order emulator of [ w i ( p )] , we expand equation 38 to the 2nd order w i ( p ) ≈ w i ( p ) + (cid:213) j = w i , p j ( p j − p , j ) + (cid:213) j = (cid:213) k = j w i , p j p k ( p j − p , j )( p k − p , k ) , (39)where w i corresponds to the integrated value of the projected2PCF over the i -th bin. w i , p j denotes the first derivativeswhile w i , p j p k denotes the second derivatives. Since for each i , we have 4 w i , p j and 10 w i , p j p k unknowns, we use a total of14 different and independent evaluations of [ w i ( p )] to solvefor all 14 derivatives.Once we have solved for the first and second deriva-tives on the projected 2PCF, we have determined the fullparametrized form of our nd order emulator for [ w i ( p )] .Then we vary the HOD parameters by first varying α (wechoose α because it controls the number of satellites) by + α log M cut log M σ Table 2.
The HOD parameters corresponding to δα = ± . thatbest reproduce the Z09 design 2PCF. and − and then perturbing the remaining HOD parame-ters to find the [ w i ( p )] that best emulates [ w i ( p )] . Specifi-cally, we minimize χ with an error on δ w i that scales withinverse sample volume. Considering that our full sample livein a sphere of diameter d p , max in pair separation space, theerror is given by σ i = (cid:113) d p , max − d ⊥ , i / d ⊥ , i for an even binningin d ⊥ , i .For both δα = . and δα = − . , we find local min-ima of χ from an arbitrarily selected set of initial locationsin the HOD parameter space. Then we visually select the setof HOD parameters that best emulates the Z09 2PCF. Ta-ble 2 shows the resulting parameters of these best fit HODs.The minimization is done using the Powell algorithm (Pow-ell 1964), Nelder-Mead algorithm and the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm (Shanno 1970; Shanno& Kettler 1970), which yield consistent best fits.We plot these two best fit HODs together with the Z09design in Figure 11. The number of centrals and satellitesare plotted in dotted red curve. The HOD correspondingto α = . is plotted in blue and the HOD correspondingto α = . is plotted in magenta. The solid curves show n cent and the dashed curves show n sat . The bottom panelshows the absolute difference between the α = . HODand the α = . HOD, where again the solid line shows thedifference in n cent and the dashed line shows the differencein n sat . We see that by increasing α , we are increasing thenumber of satellites at the high-mass end of the HOD asexpected. To compensate for such changes so as to maintainthe same 2PCF, we see a decrease in the number of centralsat halo mass M < M (cid:12) . Thus, while we are increasingthe number of LRG pairs in the largest halos, our algorithmis removing central-central pairs in the more common butsmaller halos to compensate.Figure 12 top panel shows the projected 2PCF corre-sponding to the Z09 design. The bottom panel shows the rel-ative change in the projected 2PCFs after perturbing HODparameters. The dashed curves give the relative change inprojected 2PCF after only changing α , with blue correspond-ing to α = . and magenta corresponding to α = . . Thesolid curves show the relative change in projected 2PCF af-ter also perturbing the other HOD parameters to suppressthe deviations in the projected 2PCF from that of the Z09design. We see that for a change in α , the best-fit 2PCFsshows an offset of ∼ at most, with the largest fractionaldeviations occurring between d ⊥ ≈ Mpc. The fractionaldeviation is below beyond 4 Mpc. This change is con-siderably smaller than the ∼ change that occurs in theprojected 2PCF before we adjust the other HOD parame-ters. MNRAS000
The HOD parameters corresponding to δα = ± . thatbest reproduce the Z09 design 2PCF. and − and then perturbing the remaining HOD parame-ters to find the [ w i ( p )] that best emulates [ w i ( p )] . Specifi-cally, we minimize χ with an error on δ w i that scales withinverse sample volume. Considering that our full sample livein a sphere of diameter d p , max in pair separation space, theerror is given by σ i = (cid:113) d p , max − d ⊥ , i / d ⊥ , i for an even binningin d ⊥ , i .For both δα = . and δα = − . , we find local min-ima of χ from an arbitrarily selected set of initial locationsin the HOD parameter space. Then we visually select the setof HOD parameters that best emulates the Z09 2PCF. Ta-ble 2 shows the resulting parameters of these best fit HODs.The minimization is done using the Powell algorithm (Pow-ell 1964), Nelder-Mead algorithm and the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm (Shanno 1970; Shanno& Kettler 1970), which yield consistent best fits.We plot these two best fit HODs together with the Z09design in Figure 11. The number of centrals and satellitesare plotted in dotted red curve. The HOD correspondingto α = . is plotted in blue and the HOD correspondingto α = . is plotted in magenta. The solid curves show n cent and the dashed curves show n sat . The bottom panelshows the absolute difference between the α = . HODand the α = . HOD, where again the solid line shows thedifference in n cent and the dashed line shows the differencein n sat . We see that by increasing α , we are increasing thenumber of satellites at the high-mass end of the HOD asexpected. To compensate for such changes so as to maintainthe same 2PCF, we see a decrease in the number of centralsat halo mass M < M (cid:12) . Thus, while we are increasingthe number of LRG pairs in the largest halos, our algorithmis removing central-central pairs in the more common butsmaller halos to compensate.Figure 12 top panel shows the projected 2PCF corre-sponding to the Z09 design. The bottom panel shows the rel-ative change in the projected 2PCFs after perturbing HODparameters. The dashed curves give the relative change inprojected 2PCF after only changing α , with blue correspond-ing to α = . and magenta corresponding to α = . . Thesolid curves show the relative change in projected 2PCF af-ter also perturbing the other HOD parameters to suppressthe deviations in the projected 2PCF from that of the Z09design. We see that for a change in α , the best-fit 2PCFsshows an offset of ∼ at most, with the largest fractionaldeviations occurring between d ⊥ ≈ Mpc. The fractionaldeviation is below beyond 4 Mpc. This change is con-siderably smaller than the ∼ change that occurs in theprojected 2PCF before we adjust the other HOD parame-ters. MNRAS000 , 1–16 (2017)
PCF in squeezed limit Figure 11.
The HODs for α = . (in blue) and α = . (inmagenta) that best fit the Z09 design HOD (in red). For the Z09design, both n cent and n sat are plotted in dotted line. For the othertwo HODs, n cent is plotted in solid curves whereas n sat is plottedin dashed curves. The bottom panel shows the difference betweenthe α = . HOD and the α = . HOD, where again the solidline shows the difference in n cent and the dashed line shows thedifference in n sat . Now we compute the change in the reduced 3PCF Q eff cor-responding to the two best-fit HODs that we found in thelast section.We plot the relative change in Q eff in Figure 13, where δ Q eff is defined as δ Q eff = [ Q eff ( α = . ) − Q eff ( α = . )] Q eff ( α = ) . (40)The δ Q eff signal is averaged over 16 simulation boxes, eachwith 16 different seeds to suppress shot noise. We see a rel-ative change in Q eff not greater than at d ⊥ ∼ d (cid:107) . Thissuggests that Q eff is most sensitive to δα along the 2 Mpcridge (refer to the bottom panel of Figure 9) that corre-sponds to the most massive halos. This makes sense since α controls the number of galaxies in massive halos.In detail, the morphology of δ Q eff shown in Figure 13recalls the variations in the projected 2PCF (Figure 12): anexcess of at 2 Mpc scale, dominating at large d (cid:107) wherethe most massive halo pairs occur, followed by < differ- Figure 12.
The top panel shows the projected 2PCFs corre-sponding to the Z09 design. The bottom panel show the changein projected 2PCFs corresponding to perturbed HODs, given interms of ∆ [ d ⊥ w ( d ⊥ )] = d ⊥ w perturbed − d ⊥ w Z09 . The dashed curvescorrespond to only changing HOD parameter α , with α = . plotted in blue and α = . plotted in magenta. The solid curvescorrespond to the best-fit HODs after perturbing the other HODparameters to best reproduce the original 2PCF. The best fitssuppress the difference between the α = . case and the α = . case from < down to < . ences at larger d ⊥ . We therefore wanted to check how the Q eff variations in Figure 13 compared to the variations in the fullanisotropic 2PCF. This is shown in Figure 14, where we havecomputed ξ ( d ⊥ , d (cid:107) ) for both the α = . and α = . HODsand then shown the fractional difference. One sees that whilethe broadest patterns coincide, the details vary significantly.Specifically, the maximum variation in ξ is ∼ , twice aslarge as that of Q eff . We see a strong signal in δ Q eff at thesmallest d ⊥ but not in δξ . The ridge in δ Q eff extends allthe way down to the smallest d (cid:107) whereas the δξ signal onlyextends down to about d (cid:107) ∼ Mpc. There is also a consis-tently negative region in δξ at d ⊥ ∼ δ Q eff signal. The angular gradient at largepair-separation is also different between δ Q eff and δξ . Thus,we conclude that the reduced 3PCF Q eff ( d ⊥ , d (cid:107) ) is not a sim-ple rescaling of the anisotropic 2PCF ξ ( d ⊥ , d (cid:107) ) . However, asthe variations in the reduced 3PCF Q eff ( < ) are similarin size and scale to those in the projected 2PCF, it is unclearwhether the squeezed 3PCF is breaking a degeneracy withinthis particular space of HOD models. MNRAS , 1–16 (2017) S. Yuan, D. J. Eisenstein and L. H. Garrison
Figure 13.
The relative change in reduced 3PCF δ Q eff = [ Q eff ( α = . ) − Q eff ( α = . )]/ Q eff ( α = ) plotted as a func-tion of pair separation parallel to LOS d (cid:107) and pair separationperpendicular to LOS d ⊥ . The Q eff is given in redshift space. Figure 14.
The change in the redshift-space 2PCF as a functionof pair separation. Specifically, we define δξ ( d p ) = ( ξ ( α = . ) − ξ ( α = . ))/ ξ ( α = ) . While this test was inconclusive within this modelspace, we expect that there could be reasonable extensions ofthe HOD where the variation to Q eff is disentangled from thevariation to the projected 2PCF. For these scenarios, Q eff would be a new diagnostic to distinguish among differentHOD models. However, it is also possible that the squeezed3PCF is well predicted by the 2PCF for certain classes ofHOD models. Such a predicative relationship would be in-teresting in itself, as any discrepancy between the observedsqueezed 3PCF and the 2PCF could then be used to falsifysuch classes of HOD models.In future works, we intend to apply our methodologyto more sophisticated halo occupation models that incorpo-rate dependence on more parameters than the one we usedfor this paper, such as models that incorporate halo assem-bly bias (Wechsler et al. 2002; Gao et al. 2005; Wechsler et al. 2006; Croton et al. 2007; Zentner et al. 2014; Lehmannet al. 2017). In particular, Hearin et al. (2016) parameter-ize the assembly bias in an HOD-like model (“DecoratedHOD”). Recently, there has also been a series of studies ex-ploring the dependence of the HOD on halo environment(e.g. Mandelbaum et al. 2006; Gil-Mar´ın et al. 2011; Croftet al. 2012; Tonnesen & Cen 2015) and halo geometry (vanDaalen et al. 2012). It is plausible that the pair-galaxy biasand the squeezed 3PCF can offer extra constraints for thesesophisticated halo occupation models. To summarize, this paper presents the methodology of cross-correlating the galaxy pair field with the galaxy field as asqueezed 3PCF. We develop a fast and efficient method us-ing FFTs to compute the pair-galaxy bias b p , g and its as-sociated reduced squeezed 3PCF Q eff . We implement ourmethod using a series of N-body cosmology simulations andshow the b p , g and Q eff as functions of pair separation inreal space and in redshift space. In real space, a significantpeak in b p , g and Q eff at pair separation ∼ ∼ b p , g and Q eff signals that tracksthe large-scale filamentary structure.In redshift space, the FoG effect smears out the 2 Mpcpeak into a ridge along the LOS. The anisotropy signal isalso smeared out along the LOS, but persists at large pairseparation. Both in real space and redshift, we see a factorof 2 variation in the squeezed 3PCF, contradicting the hi-erarchical ansatz. However, the variations in the squeezed3PCF offers rich information on the galaxy-halo connectionand the large-scale filamentary structure.Thus, we investigate using the reduced squeezed 3PCF Q eff to constrain the high mass end of the HOD while break-ing degeneracies in the 2PCF. We develop a second order em-ulator to model the projected 2PCF around the Z09 designHOD. Using the emulator, we find two simple HOD prescrip-tions that constrain the variations in the projected 2PCF to < when handed a perturbation to HOD parameter α . The resulting perturbation to the squeezed 3PCF is nogreater than , which cannot be immediately disentangledfrom the variations in the 2PCF. Nevertheless, we proposethat more sophisticated HOD models or different variationsto the HOD could further disentangle the squeezed 3PCFfrom the 2PCF.In future papers, we intend to construct a more sophis-ticated emulator of the 2PCF as a function of HOD pa-rameters. Another key objective is to extend our methodol-ogy to more sophisticated halo occupation models, such asones that incorporate assembly bias, and explore how thepair-galaxy bias and the squeezed 3PCF can be used to dis-tinguish among these more complex models. Eventually wewant to apply our methodology to observations to detectthese signals. ACKNOWLEDGEMENTS
We would also like to thank Dr. Zachary Slepian and JoshuaSpeagle for useful discussions. DJE and LG are supported
MNRAS000
MNRAS000 , 1–16 (2017)
PCF in squeezed limit by National Science Foundation grant AST-1313285. DJE isadditionally supported by U.S. Department of Energy grantDE-SC0013718 and as a Simons Foundation Investigator. REFERENCES
Abazajian K., et al., 2003, AJ, 126, 2081Baldauf T., Seljak U., Senatore L., 2011, J. Cosmology Astropart.Phys., 4, 006Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013, ApJ, 762, 109Berlind A. A., Weinberg D. H., 2002, ApJ, 575, 587Berlind A. A., et al., 2003, ApJ, 593, 1Bernardeau F., Colombi S., Gazta˜naga E., Scoccimarro R., 2002,Phys. Rep., 367, 1Blake C., Collister A., Lahav O., 2008, MNRAS, 385, 1257Blanton M. R., Eisenstein D., Hogg D. W., Zehavi I., 2006, ApJ,645, 977Cooray A., Sheth R., 2002, Phys. Rep., 372, 1Croft R. A. C., Matteo T. D., Khandai N., Springel V., Jana A.,Gardner J. P., 2012, MNRAS, 425, 2766Croton D. J., Gao L., White S. D. M., 2007, MNRAS, 374, 1303Dalal N., Dor´e O., Huterer D., Shirokov A., 2008, Phys. Rev. D,77, 123514Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ,292, 371Dawson K. S., et al., 2013, AJ, 145, 10Eisenstein D. J., Blanton M., Zehavi I., Bahcall N., BrinkmannJ., Loveday J., Meiksin A., Schneider D., 2005a, ApJ, 619,178Eisenstein D. J., et al., 2005b, ApJ, 633, 560Eisenstein D. J., et al., 2011, AJ, 142, 72Fergusson J. R., Shellard E. P. S., 2009, Phys. Rev. D, 80, 043510Fry J. N., 1996, ApJ, 461, L65Gao L., Springel V., White S. D. M., 2005, MNRAS, 363, L66Garrison L. H., Eisenstein D. J., Ferrer D., Metchnik M. V., PintoP. A., 2016, MNRAS, 461, 4125Gazta˜naga E., Scoccimarro R., 2005, MNRAS, 361, 824Gazta˜naga E., Cabr´e A., Castander F., Crocce M., Fosalba P.,2009, MNRAS, 399, 801Gaztanaga E., Frieman J. A., 1994, ApJ, 437, L13Gil-Mar´ın H., Jimenez R., Verde L., 2011, MNRAS, 414, 1207Gil-Mar´ın H., Nore˜na J., Verde L., Percival W. J., Wagner C.,Manera M., Schneider D. P., 2015, MNRAS, 451, 539Groth E. J., Peebles P. J. E., 1977, ApJ, 217, 385Guo H., et al., 2013, ApJ, 767, 122Guo H., et al., 2015a, MNRAS, 446, 578Guo H., et al., 2015b, MNRAS, 449, L95Guo H., et al., 2016, ApJ, 831, 3Hamilton A. J. S., 1998, in Hamilton D., ed., Astrophysicsand Space Science Library Vol. 231, The Evolving Universe.p. 185 ( arXiv:astro-ph/9708102 ), doi:10.1007/978-94-011-4960-0 17Hashimoto I., Taruya A., Matsubara T., Namikawa T., YokoyamaS., 2016, Phys. Rev. D, 93, 103537Hearin A. P., Zentner A. R., van den Bosch F. C., Campbell D.,Tollerud E., 2016, MNRAS, 460, 2552Ho S., Lin Y.-T., Spergel D., Hirata C. M., 2009, ApJ, 697, 1358Hoshino H., et al., 2015, MNRAS, 452, 998Jeong D., 2010, PhD thesis, University of Texas at Austin,[email protected] Y. P., B¨orner G., 2004, ApJ, 607, 140Kaiser N., 1987, MNRAS, 227, 1Kauffmann G., White S. D. M., Heckman T. M., M´enard B.,Brinchmann J., Charlot S., Tremonti C., Brinkmann J., 2004,MNRAS, 353, 713Knebe A., et al., 2011, MNRAS, 415, 2293 Kravtsov A. V., Berlind A. A., Wechsler R. H., Klypin A. A.,Gottl¨ober S., Allgood B., Primack J. R., 2004, ApJ, 609, 35Kulkarni G. V., Nichol R. C., Sheth R. K., Seo H.-J., EisensteinD. J., Gray A., 2007, MNRAS, 378, 1196Kwan J., Heitmann K., Habib S., Padmanabhan N., Lawrence E.,Finkel H., Frontiere N., Pope A., 2015, ApJ, 810, 35Leauthaud A., et al., 2012, ApJ, 744, 159Lehmann B. V., Mao Y.-Y., Becker M. R., Skillman S. W., Wech-sler R. H., 2017, ApJ, 834, 37Li R., Mo H. J., Fan Z., Cacciato M., van den Bosch F. C., YangX., More S., 2009, MNRAS, 394, 1016Lukic Z., 2008, PhD thesis, University of Illinois at Urbana-ChampaignMandelbaum R., Seljak U., Kauffmann G., Hirata C. M.,Brinkmann J., 2006, MNRAS, 368, 715Mandelbaum R., Slosar A., Baldauf T., Seljak U., Hirata C. M.,Nakajima R., Reyes R., Smith R. E., 2013, MNRAS, 432, 1544Mar´ın F., 2011, ApJ, 737, 97Mar´ın F. A., Wechsler R. H., Frieman J. A., Nichol R. C., 2008,ApJ, 672, 849McBride C. K., Connolly A. J., Gardner J. P., Scranton R., New-man J. A., Scoccimarro R., Zehavi I., Schneider D. P., 2011,ApJ, 726, 13Miyatake H., et al., 2015, ApJ, 806, 1Mo H. J., Yang X., van den Bosch F. C., Jing Y. P., 2004, MN-RAS, 349, 205More S., Miyatake H., Mandelbaum R., Takada M., Spergel D. N.,Brownstein J. R., Schneider D. P., 2015, ApJ, 806, 2Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493Parejko J. K., et al., 2013, MNRAS, 429, 98Peacock J. A., Smith R. E., 2000, MNRAS, 318, 1144Peebles P. J. E., 1980, The large-scale structure of the universePeebles P. J. E., Groth E. J., 1975, ApJ, 196, 1Planck Collaboration et al., 2016, A&A, 594, A13Powell M. J. D., 1964, The Computer Journal, 7, 155Reid B. A., Spergel D. N., 2009, ApJ, 698, 143Rodr´ıguez-Torres S. A., et al., 2016, MNRAS, 460, 1173Scherrer R. J., Weinberg D. H., 1998, ApJ, 504, 607Scoccimarro R., Sheth R. K., Hui L., Jain B., 2001, ApJ, 546, 20Sefusatti E., Crocce M., Desjacques V., 2010, MNRAS, 406, 1014Seljak U., 2000, MNRAS, 318, 203Seljak U., Zaldarriaga M., 1996, ApJ, 469, 437Shanno D. F., 1970, Mathematics of computation, 24, 647Shanno D. F., Kettler P. C., 1970, Mathematics of Computation,24, 657Slepian Z., et al., 2017a, MNRAS, 468, 1070Slepian Z., et al., 2017b, MNRAS, 469, 1738Soneira R. M., Peebles P. J. E., 1976, in Bulletin of the AmericanAstronomical Society. p. 354Tasinato G., Tellarini M., Ross A. J., Wands D., 2014, J. Cosmol-ogy Astropart. Phys., 3, 032Tegmark M., Peebles P. J. E., 1998, ApJ, 500, L79Tellarini M., Ross A. J., Tasinato G., Wands D., 2016, J. Cosmol-ogy Astropart. Phys., 6, 014Tinker J. L., Conroy C., Norberg P., Patiri S. G., Weinberg D. H.,Warren M. S., 2008, ApJ, 686, 53Tinker J. L., et al., 2012, ApJ, 745, 16Tonnesen S., Cen R., 2015, ApJ, 812, 104Wechsler R. H., Bullock J. S., Primack J. R., Kravtsov A. V.,Dekel A., 2002, ApJ, 568, 52Wechsler R. H., Zentner A. R., Bullock J. S., Kravtsov A. V.,Allgood B., 2006, ApJ, 652, 71White M., et al., 2011, ApJ, 728, 126Yoo J., Tinker J. L., Weinberg D. H., Zheng Z., Katz N., Dav´eR., 2006, ApJ, 652, 26York D. G., et al., 2000, AJ, 120, 1579Zaldarriaga M., Seljak U., 2000, ApJS, 129, 431MNRAS , 1–16 (2017) S. Yuan, D. J. Eisenstein and L. H. Garrison
Zaldarriaga M., Seljak U., Bertschinger E., 1998, ApJ, 494, 491Zehavi I., et al., 2002, ApJ, 571, 172Zehavi I., et al., 2004, ApJ, 608, 16Zehavi I., et al., 2005a, ApJ, 621, 22Zehavi I., et al., 2005b, ApJ, 630, 1Zentner A. R., Hearin A. P., van den Bosch F. C., 2014, MNRAS,443, 3044Zheng Z., Weinberg D. H., 2007, ApJ, 659, 1Zheng Z., et al., 2005, ApJ, 633, 791Zheng Z., Coil A. L., Zehavi I., 2007, ApJ, 667, 760Zheng Z., Zehavi I., Eisenstein D. J., Weinberg D. H., Jing Y. P.,2009, ApJ, 707, 554van Daalen M. P., Angulo R. E., White S. D. M., 2012, MNRAS,424, 2954 MNRAS000