Using galaxy pairs to investigate the three-point correlation function in the squeezed limit

TL;DR

This paper develops an efficient FFT-based method to analyze the squeezed limit of the three-point correlation function using galaxy pairs, revealing insights into galaxy-halo connections and the impact of large-scale structures.

Contribution

It introduces a novel algorithm for computing the squeezed 3PCF with galaxy pairs, applied to simulations, and explores its potential to constrain complex HOD models.

Findings

Peak in pair-galaxy bias at ~2 Mpc separation

Strong anisotropy in signals tracking filamentary structures

Squeezed 3PCF varies significantly, offering additional galaxy-halo information

Abstract

We investigate the three-point correlation function (3PCF) in the squeezed limit by considering galaxy pairs as discrete objects and cross-correlating them with the galaxy field. We develop an efficient algorithm using Fast Fourier Transforms to compute such cross-correlations and their associated pair-galaxy bias bpg and the squeezed 3PCF coefficient Qeff. We implement our method using N-body cosmological simulations and a fiducial Halo Occupation Distribution (HOD) and present the results in both the real space and redshift space. In real space, we observe a peak in bpg and Qeff at pair separation of ~ 2 Mpc, attributed to the fact that galaxy pairs at 2 Mpc separation trace the most massive dark matter halos. We also see strong anisotropy in the bpg and Qeff 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, the squeezed 3PCF shows a factor of 2 variation, contradicting the hierarchical ansatz but offering rich information on the galaxy-halo connection. Thus, we explore the possibility of using the squeezed 3PCF to constrain the HOD. When we compare two simple HOD models that are closely matched in their projected two-point correlation function (2PCF), we do not yet see a strong variation in the 3PCF that is clearly disentangled from variations in the projected 2PCF. Nevertheless, we propose that more complicated HOD models, e.g. those incorporating assembly bias, can break degeneracies in the 2PCF and show a distinguishable squeezed 3PCF signal.

Tables

Table 1. Table 1: Parameters of the weighting functions for d p , max = 10 subscript 𝑑 𝑝 max 10 d_{p,\rm{max}}=10  Mpc and d p , max = 30 subscript 𝑑 𝑝 max 30 d_{p,\rm{max}}=30  Mpc. Again the cutoff is given by r cut ≈ d p , max / 2 + 2 ​ Δ ​ x subscript 𝑟 cut subscript 𝑑 𝑝 max 2 2 Δ 𝑥 r_{\rm{cut}}\approx d_{p,\rm{max}}/2+2\Delta x .

$d_{p,\max }$ (Mpc)	$\sigma$ (Mpc)	$r_s$ (Mpc)	$s_0$ (Mpc-1)	$r_{\mathrm{cut}}$ (Mpc)	$d_w$ (Mpc)	$\xi \left(d_w\right)$
10	30	20	0.8	10	30.9	0.362
30	50	35	0.5	20	52.6	0.104

Table 2. Table 2: The HOD parameters corresponding to δ ​ α = ± 0.05 𝛿 𝛼 plus-or-minus 0.05 \delta\alpha=\pm 0.05 that best reproduce the Z09 design 2PCF.

$\alpha$	${\log }_{10}{M}_{\mathrm{cut}}$	${\log }_{10}{M}_1$	$\sigma$
1+0.05	13.35–0.086	13.8–0.018	0.85–0.16
1–0.05	13.35+0.070	13.8–0.024	0.85+0.093

Equations

$$b_{p,g}={\xi }_{p,g}/{\xi }_{g,g},$$

$${\chi }^2=\frac{|\widetilde{{\delta }_p}(\text{k})-b_{p,g}\widetilde{{\delta }_g}(\text{k}){|}^2}{{\sigma }^2}.$$

$$b_{p,g}(\text{k})=\frac{\text{Re}({\widetilde{{\delta }_g}}^{\ast }(\text{k})\widetilde{{\delta }_p}(\text{k}))}{{\widetilde{{\delta }_g}}^{\ast }(\text{k})\widetilde{{\delta }_g}(\text{k})},{\sigma }_b^{-2}(\text{k})=\frac{{\widetilde{{\delta }_g}}^{\ast }(\text{k})\widetilde{{\delta }_g}(\text{k})}{{\sigma }^2}.$$

$$b_{p,g}$$

$$=\frac{{\sum }_{\text{k}}\text{Re}({\widetilde{{\delta }_g}}^{\ast }(\text{k})\widetilde{{\delta }_p}(\text{k}))\tilde{W}(\text{k})}{{\sum }_{\text{k}}{\widetilde{{\delta }_g}}^{\ast }(\text{k})\widetilde{{\delta }_g}(\text{k})\tilde{W}(\text{k})}$$

$$=\frac{\int d^3\text{k}\tilde{W}(\text{k})\text{Re}({\widetilde{{\delta }_g}}^{\ast }(\text{k})\widetilde{{\delta }_p}(\text{k}))}{\int d^3\text{k}\tilde{W}(\text{k})P_g(\text{k})},$$

$${\delta }_{p,i}(\text{x})=\frac{{\delta }_D(\text{x}-{\text{x}}_{p,i})}{{\bar{n}}_p}-1,$$

$${\tilde{\delta }}_{p,i}(\text{k})=\frac{e^{-i\text{k}\cdot {\text{x}}_{p,i}}}{{\bar{n}}_p}-\frac{{\delta }_D(\text{k})}{(2\pi {)}^3}$$

$$\text{Re}({\widetilde{{\delta }_g}}^{\ast }(\text{k})\widetilde{{\delta }_p}(\text{k}))=\text{Re}({\widetilde{{\delta }_g}}^{\ast }(\text{k})\sum _i{\tilde{\delta }}_{p,i}(\text{k}))$$

$$=\frac{1}{2}\sum _i({\widetilde{{\delta }_g}}^{\ast }(\text{k}){\tilde{\delta }}_{p,i}(\text{k})+\widetilde{{\delta }_g}(\text{k}){\tilde{\delta }}_{p,i}^{\ast }(\text{k}))$$

$$=\frac{1}{2{\bar{n}}_p}\sum _i({\widetilde{{\delta }_g}}^{\ast }(\text{k})e^{-i\text{k}\cdot {\text{x}}_{p,i}}+\widetilde{{\delta }_g}(\text{k})e^{i\text{k}\cdot {\text{x}}_{p,i}})$$

$$-\frac{N_p}{2}\frac{{\delta }_D(\text{k})}{(2\pi {)}^3}(\widetilde{{\delta }_g}(\text{k})+{\widetilde{{\delta }_g}}^{\ast }(\text{k})),$$

$$\widetilde{{\delta }_g}(\text{0})={\widetilde{{\delta }_g}}^{\ast }(\text{0})=\int d^3\text{x}{\delta }_g(\text{x})=0.$$

$$b_{p,g}\left[\int d^3\text{k}\tilde{W}(\text{k})P_g(\text{k})\right]$$

$$=\frac{1}{2{\bar{n}}_p}\sum _i\left[\int d^3\text{k}\tilde{W}(\text{k})({\widetilde{{\delta }_g}}^{\ast }(\text{k})e^{-i\text{k}\cdot {\text{x}}_{p,i}}+\widetilde{{\delta }_g}(\text{k})e^{i\text{k}\cdot {\text{x}}_{p,i}})\right]$$

$$=\frac{1}{2{\bar{n}}_p}\sum _i[\int d^3\text{k}\tilde{W}(-\text{k})\widetilde{{\delta }_g}(-\text{k})e^{-i\text{k}\cdot {\text{x}}_{p,i}}$$

$$+\int d^3\text{k}\tilde{W}(\text{k})\widetilde{{\delta }_g}(\text{k})e^{i\text{k}\cdot {\text{x}}_{p,i}}]$$

$$=\frac{1}{{\bar{n}}_p}\sum _i\left[\int d^3\text{k}\tilde{W}(\text{k})\widetilde{{\delta }_g}(\text{k})e^{i\text{k}\cdot {\text{x}}_{p,i}}\right]$$

$$=\frac{1}{{\bar{n}}_p}\sum _i(W\ast {\delta }_g)({\text{x}}_{p,i}).$$

$$\int d^3\text{k}\tilde{W}(\text{k})P_g(\text{k})=\frac{1}{{\bar{n}}_g}\sum _j(W\ast {\delta }_g)({\text{x}}_{g,j}),$$

$$b_{p,g}=\frac{\frac{1}{{\bar{n}}_p}{\sum }_i(W\ast {\delta }_g)({\text{x}}_{p,i})}{\frac{1}{{\bar{n}}_g}{\sum }_j(W\ast {\delta }_g)({\text{x}}_{g,j})}.$$

$$\frac{1}{{\bar{n}}_g}\sum _j(W\ast {\delta }_g)({\text{x}}_{g,j})$$

$$=\frac{1}{{\bar{n}}_g}\sum _j\int d^3{\text{x}}^{\prime }{\delta }_g({\text{x}}^{\prime }+{\text{x}}_{g,j})W({\text{x}}^{\prime })$$

$$=\frac{1}{{\bar{n}}_g}\sum _j\int d^3{\text{x}}^{\prime }\left(\frac{n_g({\text{x}}^{\prime }+{\text{x}}_{g,j}|{\text{x}}_{g,j})}{{\bar{n}}_g}-1\right)W({\text{x}}^{\prime }),$$

$$n_g({\text{x}}^{\prime }+{\text{x}}_{g,j}|{\text{x}}_{g,j})={\bar{n}}_g(1+\xi ({\text{x}}^{\prime }))$$

$$\frac{1}{{\bar{n}}_g}\sum _j(W\ast {\delta }_g)({\text{x}}_{g,j})=\int d^3{\text{x}}^{\prime }\xi ({\text{x}}^{\prime })W({\text{x}}^{\prime }).$$

$$\frac{1}{{\bar{n}}_p}\sum _i(W\ast {\delta }_g)({\text{x}}_{p,i})$$

$$=\frac{1}{{\bar{n}}_p}\sum _i\int d^3{\text{x}}^{\prime }\left(\frac{n_g({\text{x}}^{\prime }+{\text{x}}_{p,i}|{\text{x}}_{p,i})}{{\bar{n}}_g}-1\right)W({\text{x}}^{\prime }),$$

$$n_g({\text{x}}^{\prime }+{\text{x}}_{p,i}|{\text{x}}_{p,i})={\bar{n}}_g(1+\xi ({\text{d}}_p)+\xi ({\text{x}}^{\prime }+\frac{{\text{d}}_p}{2})$$

$$+\xi ({\text{x}}^{\prime }-\frac{{\text{d}}_p}{2})+\zeta ({\text{d}}_p,{\text{x}}^{\prime }+\frac{{\text{d}}_p}{2},{\text{x}}^{\prime }-\frac{{\text{d}}_p}{2}))/(1+\xi ({\text{d}}_p)),$$

$$\frac{1}{{\bar{n}}_p}\sum _i(W\ast {\delta }_g)({\text{x}}_{p,i})$$

$$=\frac{1}{{\bar{n}}_p}\sum _i\int d^3{\text{x}}^{\prime }\left(\frac{\xi ({\text{x}}_1^{\prime })+\xi ({\text{x}}_2^{\prime })+\zeta ({\text{d}}_p,{\text{x}}_1^{\prime },{\text{x}}_2^{\prime })}{1+\xi ({\text{d}}_p)}\right)W({\text{x}}^{\prime })$$

$$\frac{1}{{\bar{n}}_p}\sum _i(W\ast {\delta }_g)({\text{x}}_{p,i})$$

$$=\frac{\int d^3{\text{x}}^{\prime }W({\text{x}}^{\prime })(\xi ({\text{x}}_1^{\prime })+\xi ({\text{x}}_2^{\prime }))}{1+\xi ({\text{d}}_p)}+\frac{\int d^3{\text{x}}^{\prime }W({\text{x}}^{\prime })\zeta ({\text{d}}_p,{\text{x}}_1^{\prime },{\text{x}}_2^{\prime })}{1+\xi ({\text{d}}_p)}.$$

$$b_{p,g}=\frac{A+B}{1+\xi ({\text{d}}_p)},$$

$$A$$

$$B$$

$$Q_{\mathrm{eff}}({\text{x}}_1,{\text{x}}_2,{\text{x}}_3)=\frac{\zeta ({\text{x}}_1,{\text{x}}_2,{\text{x}}_3)}{\xi ({\text{x}}_1)\xi ({\text{x}}_2)+\xi ({\text{x}}_1)\xi ({\text{x}}_3)+\xi ({\text{x}}_2)\xi ({\text{x}}_3)},$$

$$\zeta ({\text{d}}_p,{\text{x}}_1^{\prime },{\text{x}}_2^{\prime })$$

Full text

††pubyear: 2017††pagerange: Using galaxy pairs to investigate the three-point correlation function in the squeezed limit–References

1 Introduction

Understanding galaxy clustering is one of the primary goals of cosmology, and the advent of a new generation of galaxy surveys has ushered in the new era of high precision galaxy clustering studies. A key objective of such studies is to understand the relationship between galaxies and the underlying dark matter halos since such relationship is critical in helping us infer the linear dark matter power spectrum from the observed galaxy power spectrum. The Halo Occupation Distribution (HOD) is the most popular framework for describing how the number of galaxies per halo depends the halo mass (Seljak, 2000; Peacock & Smith, 2000; Scoccimarro et al., 2001; Berlind & Weinberg, 2002; Cooray & Sheth, 2002; Kravtsov et al., 2004; Zheng et al., 2005). For the simple HOD we consider in this paper, galaxies populating a halo are divided into central and satellite galaxies, hereafter referred to simply as centrals and satellites (Berlind et al., 2003; Zheng et al., 2005).

There has been a wealth of studies that utilize the two-point correlation functions (2PCF) or higher-order statistics to constrain the HOD (e.g. Zehavi et al., 2005b; Kulkarni et al., 2007; Blake et al., 2008; Zheng et al., 2009; White et al., 2011; Parejko et al., 2013; More et al., 2015; Guo et al., 2016; Rodríguez-Torres et al., 2016). Another popular way to constrain the HOD is by probing the halo mass of galaxies via weak gravitational lensing, commonly know as galaxy-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; More et al., 2015). Other methods for constraining the HOD include a pair-counting based technique called the Counts-in-Cylinders (CiC) (Reid & Spergel, 2009), and direct counting of the number of galaxies in clusters as a function of halo mass (Ho et al., 2009; Hoshino et al., 2015).

For this paper, we are interested in constraining the high mass regime of the HOD ($M\sim 10^{14-15}M_{\odot}$), where each dark matter halo hosts multiple galaxies close together. This suggests that if we consider close pairs of galaxies as objects in their own right, they would be biased tracers of high mass halos (refer to Figure 4).

The first step is to formulate the proper statistic on close galaxy pairs to calculate their bias. One might choose to use the pair-pair auto-correlation function. However, to preserve the pair separation dependence of our statistic, we need to auto-correlate much smaller subsamples of pairs binned by pair separation. Auto-correlating smaller samples leads to a noisier statistic. Thus, instead of auto-correlating close pairs, we investigate the cross-correlation of close pairs with galaxy singlets to calculate the pair-galaxy cross-correlation function. The large-scale pair-galaxy cross-correlation then reveals the linear bias (Fry, 1996; Tegmark & Peebles, 1998; Scherrer & Weinberg, 1998) of the pair field relative to the galaxy field, i.e. the pair-galaxy bias, which we denote as $b_{p,g}$, given by

[TABLE]

where $\xi_{p,g}$ and $\xi_{g,g}$ are the pair-galaxy cross-correlation and the galaxy-galaxy cross-correlation, respectively. When the pair separation is much less than the distance to the galaxy singlet, the pair-galaxy cross-correlation is equivalent to the three-point correlation function (3PCF) in the squeezed limit. For conciseness, we refer to the 3PCF in the squeezed limit as the squeezed 3PCF in the rest of the paper.

Generally speaking, the 3PCF, $\zeta(\textbf{r}_{1},\textbf{r}_{2},\textbf{r}_{3})$, describes the probability distribution of finding three galaxies forming a triangle with three sides given by vectors $\textbf{r}_{1}$, $\textbf{r}_{2}$, and $\textbf{r}_{3}$ (See Bernardeau et al. (2002) for a comprehensive review). In the hierarchy of correlation functions, the 3PCF is the second lowest order correlation function after the 2PCF, and is the lowest order correlation function to provide shape information (Peebles, 1980). The squeezed limit of the 3PCF refers to the limit where one side of the triangle is much shorter than the other two sides, i.e. $r_{1}\ll r_{2}\approx r_{3}$.

The spatial distribution of matter in the Universe can be fully described by the 2PCF if the matter density field is fully Gaussian. However, we expect non-linear gravitational evolution to produce non-Gaussian signatures in the galaxy distribution that we measure today (Bernardeau et al., 2002). The 3PCF is thus needed to complement the information provided by the 2PCF when describing late-time galaxy distribution. A series of recent papers have focused on measuring and analysing the 3PCF in modern cosmological surveys, such as Gaztañaga et al. (2009); Marín (2011); Gil-Marín et al. (2015); Slepian et al. (2017a, b). The additional information provided by the 3PCF is often used to break model degeneracies describing cosmology and galaxy bias (Gaztanaga & Frieman, 1994; Jing & Börner, 2004; Kulkarni et 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 primordial non-gaussianities as a test of single-field slow-roll inflationary 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 nice review of this topic.

Measuring the 3PCF is computationally intense as the number of triplets scales as $N^{3}$ and thus becomes prohibitive for large samples of galaxies. The 3PCF in the squeezed limit is advantageous because its computation can be reduced to two separate pair finding problems of complexity $O(N^{2})$. Thus, even though we phrased our statistic as a cross-correlation between the pair field and the galaxy field, the same statistic can be interpreted as a squeezed 3PCF. We show the equivalence of the two approaches in Section 2.2.

In the following sections of the paper, we present a fast and efficient way of calculating the pair-galaxy bias using Fast Fourier Transforms (FFTs) in Section 2.1, and how it translates to the squeezed 3PCF in Section 2.2. We discuss implementing our methodology for Luminous Red Galaxies (LRGs) using N-body cosmological simulations in Section 3, and present the results in Section 4. In Section 5, we discuss the possibility of constraining HOD with the pair-galaxy bias and the squeezed 3PCF. Finally, we draw a few concluding remarks in Section 6.

2 Methodology

In this section, we develop a fast and efficient method of computing $b_{p,g}$ using FFT and then formulate the exact relationship between $b_{p,g}$ and the squeezed 3PCF starting from the triplet counting definition of the 3PCF.

2.1 Computing pair-galaxy bias $b_{p,g}$

Consider a galaxy overdensity field $\delta_{g}(\textbf{x})$, a galaxy pair overdensity field $\delta_{p}(\textbf{x})$, and their Fourier space counterparts $\tilde{\delta_{g}}(\textbf{k})$ and $\tilde{\delta_{p}}(\textbf{k})$. Assuming white noise denoted by $\sigma$, for a single wave vector k, we can write out a $\chi^{2}$ objective function

[TABLE]

Then by minimizing $\chi^{2}$ against $b_{p,g}$, we find a best fit for $b_{p,g}$ and its uncertainties at wave vector k

[TABLE]

Then we average over a range of $k$ modes, weighted by inverse variance multiplied with some weighting function $\tilde{W}(k)$. We get

[TABLE]

where $P_{g}(\textbf{k})=\tilde{\delta_{g}}^{*}(\textbf{k})\tilde{\delta_{g}}(\textbf{k})$ is the galaxy power spectrum. A detailed discussion on the choice of weighting function $\tilde{W}(k)$ is provided in Section 3.6.

Equation 4 shows that $b_{p,g}$ is linear in $\delta_{p}$, which means that we can compute the integrals for a single pair at a time and then sum over all pairs to obtain $b_{p,g}$. For the $i$-th galaxy pair, we can define its position $\textbf{x}_{p,i}$ as the middle of the pair. Then the fractional overdensity function of pairs is given by

[TABLE]

where $\delta_{D}(\textbf{x})$ is the Dirac delta function. $\bar{n}_{p}$ denotes the mean number density of galaxy pairs. The Fourier Transform is given by

[TABLE]

Then we calculate the value $\textrm{Re}(\tilde{\delta_{g}}^{*}(\textbf{k})\tilde{\delta_{p}}(\textbf{k}))$ in equation 4 as

[TABLE]

where $N_{p}$ is the total number of galaxy pairs. When we plug the second term of equation 7 into equation 4, we get a result proportional to $(\tilde{\delta_{g}}(\textbf{0})+\tilde{\delta_{g}}^{*}(\textbf{0}))$, which gives

[TABLE]

Thus, the second term of equation 7 does not contribute to the numerator in equation 4.

By combining equation 7 and equation 4, we get

[TABLE]

Note that during this derivation, we invoke that $W(\textbf{x})$ is even and real and that $\delta_{g}(\textbf{x})$ is real. Similarly, we can evaluate the denominator of equation 4

[TABLE]

where we sum over all the galaxy positions $\textbf{x}_{g,j}$, and $\bar{n}_{g}$ is the mean galaxy number density.

Thus, we get the final expression for pair-galaxy bias

[TABLE]

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 convolution for the full pair sample and simply evaluate it at the positions of any subsample of pairs whereas the denominator also only requires one convolution. This is much faster than explicitly cross-correlating each subsample of pairs with the overall galaxy field.

2.2 Relating $b_{p,g}$ to the the squeezed limit of the 3PCF

To relate $b_{p,g}$ to the squeezed 3PCF, we first evaluate the denominator of equation 11

[TABLE]

where $n_{g}(\textbf{x}^{\prime}+\textbf{x}_{g,j}|\textbf{x}_{g,j})$ is the number of galaxies at $\textbf{x}^{\prime}+\textbf{x}_{g,j}$ given an arbitrary galaxy at $\textbf{x}_{g,j}$. By definition of 2PCF $\xi(\textbf{x}^{\prime})$, we have

[TABLE]

Thus, combining equation 12 and equation 13, we have

[TABLE]

Similarly to equation 12, we can evaluate the numerator of equation 11 as

[TABLE]

where $n_{g}(\textbf{x}^{\prime}+\textbf{x}_{p,i}|\textbf{x}_{p,i})$ is the number of galaxies at $\textbf{x}^{\prime}+\textbf{x}_{p,i}$ given an arbitrary galaxy pair at position $\textbf{x}_{p,i}$. This number relates to the 3PCF by the following

[TABLE]

where $\textbf{d}_{p}$ is the separation vector of the galaxy pair. Now we combine equation 15 and equation 16 in the squeezed limit. Then, we have

[TABLE]

where we have defined $\textbf{x}^{\prime}_{1}=\textbf{x}^{\prime}+\textbf{d}_{p}/2$ and $\textbf{x}^{\prime}_{2}=\textbf{x}^{\prime}-\textbf{d}_{p}/2$.

For our purpose, we are only summing over pairs within a specific $\textbf{d}_{p}$ bin. Thus, to good approximation, we can treat $\textbf{d}_{p}$ as approximately constant in the summation, which removes all dependence on pair position inside the summation. Then we have

[TABLE]

Then, using equation 10, equation 11, and equation 18, we can express the pair-galaxy bias in terms of the 2PCF and the 3PCF

[TABLE]

where

[TABLE]

So now we have expressed the pair-galaxy bias as given by equation 11 in terms of general case 2PCFs and 3PCFs.

To proceed, we introduce the reduced 3PCF $Q_{\rm{eff}}$, defined by

[TABLE]

where $\textbf{x}_{1},\textbf{x}_{2},\textbf{x}_{3}$ are the three separation vectors for any 3PCF. In the original hierarchical model (Peebles & Groth, 1975; Soneira & Peebles, 1976; Groth & Peebles, 1977; Peebles, 1980), the $Q_{\rm{eff}}$ factor is a constant, which has been disproven in later surveys and simulations (e.g. Kulkarni et al., 2007; Marín et al., 2008; McBride et al., 2011). As we will see, our results in the squeezed limit also show a strongly variable $Q_{\rm{eff}}$.

In the strictly squeezed limit, $|\textbf{x}^{\prime}|\gg|\textbf{d}_{p}|$ and $\xi(\textbf{x}^{\prime}_{1})\approx\xi(\textbf{x}^{\prime}_{2})\approx\xi(\textbf{x}^{\prime})$. A typical 2PCF assumes a power-law form $\xi(\textbf{x}^{\prime})\propto x^{\prime\alpha}$, where the exponent $\alpha\approx-2$. This means that, in the strictly squeezed limit, $\xi(\textbf{x}^{\prime})\ll\xi(\textbf{d}_{p})$, and thus

[TABLE]

Then combining equation 19 and equation 21, we express pair-galaxy bias as

[TABLE]

Then we can express the reduced 3PCF $Q_{\rm eff}$ as a function of pair-galaxy bias $b_{p,g}$

[TABLE]

$Q_{\rm{eff}}$ is not a function of $\textbf{x}^{\prime}$ as a manifestation of the local bias limit, i.e. the pair separation scale is much smaller than the separation to the third galaxy so that $Q_{\rm{eff}}$’s dependence on $\textbf{x}^{\prime}$ converges to scale with $\xi(\textbf{x}^{\prime})$. Also note that when pair separation is extremely small, i.e. $|\textbf{d}_{p}|\sim 1$ Mpc, we have $\xi(\textbf{d}_{p})\gg 1$, and $Q_{\rm eff}$ would tend to $b_{p,g}/2$.

However, as we will see in Section 3, the pairs we are considering have a maximum separation of up to 10 or 30 Mpc, where the $|\textbf{x}^{\prime}|\gg|\textbf{d}_{p}|$ limit does not always hold. Specifically, we can no longer assume $\xi(\textbf{x}^{\prime}_{1})\approx\xi(\textbf{x}^{\prime}_{2})\approx\xi(\textbf{x}^{\prime})$ and also the $\xi(\textbf{x}^{\prime}_{1})\xi(\textbf{x}^{\prime}_{2})$ term contributes non-trivially to equation 21. We use the general form of equation 20

[TABLE]

Assuming $\xi(r)\propto r^{-2}$ and treating $d_{p}/x^{\prime}$ as a small number, we Taylor expand the first two terms, $\xi(\textbf{x}^{\prime}_{1})\xi(\textbf{d}_{p})$ and $\xi(\textbf{x}^{\prime}_{2})\xi(\textbf{d}_{p})$, to the fourth order $O\left[\left(d_{p}/x^{\prime}\right)^{4}\right]$. We also preserve the zero-th order approximation of the Taylor expansion of the last term $\xi(\textbf{x}^{\prime}_{1})\xi(\textbf{x}^{\prime}_{2})\approx\xi(\textbf{x}^{\prime})^{2}\propto(d_{p}/x^{\prime})^{4}$. The Taylor expanded form of the first two terms still carries a dependence on the angle between $\textbf{x}^{\prime}$ and $\textbf{d}_{p}$. As a rough estimate, we average over the angle between $\textbf{x}^{\prime}$ and $\textbf{d}_{p}$ assuming isotropy, we get

[TABLE]

and

[TABLE]

Note that since we are assuming isotropy and a simplified form for $\xi(r)$, the coefficient in front of $\xi(\textbf{x}^{\prime})^{2}$ is only an estimate of the correction. However, a correction to the correction term comes in as a higher order correction to $\zeta(\textbf{d}_{p},\textbf{x}^{\prime}_{1},\textbf{x}^{\prime}_{2})$, so we ignore these high order effects for this derivation. Then, in the not strictly squeezed limit, the pair-galaxy bias is given by

[TABLE]

where

[TABLE]

Given a weighting function and a 2PCF, the fraction between the two integrals yields a constant. For ease of notation, we express this term as the 2PCF evaluated at some characteristic separation $d_{w}$, which is dependent on the weighting function and the shape of the 2PCF. Specifically, we define

[TABLE]

To calculate the characteristic separation, we employ the power-law form of the galaxy-galaxy 2PCF, $\xi(r)\propto r^{-p}$. Then equation 27 gives the following formula to calculate $d_{w}$

[TABLE]

Having introduced $\xi(d_{w})$, the pair-galaxy bias can then be conveniently expressed with

[TABLE]

which translates to a modified estimator of the reduced 3PCF

[TABLE]

Thus the reduced 3PCF estimator only differs from the strictly squeezed limit by an extra factor of $\frac{\xi(d_{w})}{2\xi(\textbf{d}_{p})}$ on the numerator and a factor of $\frac{3}{2}\xi(d_{w})$ on the denominator. For the rest of this paper, since we mostly do not operate in the strictly squeezed limit, we use equation 32 to calculate $Q_{\rm{eff}}$.

3 Implementation

In this section, we present the implementation of our $b_{p,g}$ and $Q_{\rm{eff}}$ statistic for Luminous Red Galaxies (LRGs) using a series of N-body simulation boxes. LRGs are highly luminous early-type massive galaxies consisting of mainly old stars and quenched of star formation. These luminous galaxies are popular targets for large-scale structure studies and surveys. We study the statistics of LRGs specifically in this paper in preparation for applying our methods to Baryon Oscillation Spectroscopic Survey (BOSS; Eisenstein et al. 2011; Dawson et al. 2013) data in follow up papers.

We also present some basic characteristics of the simulation results in this section, namely the 2PCF, pair separation distribution, and halo mass distribution of the LRGs and their pairs. It should be pointed out that we are only interested in showcasing our methodology in this paper, and we do not attempt any detailed matching with observed LRG statistics.

3.1 N-body cosmology simulation

For the purpose of this analysis, we utilize the cosmological simulations (Garrison et al., in preparation) generated by the fast and high precision Abacus N-body code (Garrison et al. 2016, Ferrer et al., in preparation; Metchnik $\&$ Pinto, in preparation). Specifically we use a series of 16 cyclic boxes with Planck 2016 cosmology (Planck Collaboration et al., 2016) at redshift $z=0.5$, where each box is of size 1100 $h^{-1}$ Mpc, and contains 14403 dark matter particles of mass $4\times 10^{10}$ $h^{-1}\ M_{\odot}$. Dark matter halos are found and characterized 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 a few hundred particles (Garrison et al., 2016). For the rest of this paper, we use halo catalogues generated using the FoF halo finder with linking length $0.15b$, where $b$ is the average interparticle spacing. We have also repeated the key analysis with mock catalogues based on halos from ROCKSTAR and recovered consistent results.

Although FoF finder is fast, one main concern with FoF finder is the effect of FoF overbridging (Lukic, 2008; Knebe et al., 2011), which misidentifies a small fraction of physically separate halos as a single halo through long chain of links. This effect occurs at all scales and eventually leads to a small over-estimate of the number of galaxy pairs inhabiting the same halo. For our purpose, the correction is small and would not impact the qualitative results of this paper.

3.2 Redshift space distortion

An important effect to consider is Redshift Space Distortion (RSD, Hamilton 1998). Essentially, the peculiar velocity of the galaxies in clusters along the line-of-sight (LOS) leads to a deviation in distance measurements when using Hubble’s law. This is commonly known as the “finger-of-god" (FoG) effect. There is also the more subtle effect known as the Kaiser effect (Kaiser, 1987), where the coherent infall velocity of galaxies towards central masses introduces a flattening in the observed structure on the large scale. It is important to include RSD in our calculations to give results that are comparable to observables. To include RSD, we shift the $x_{\parallel}$ coordinate of each galaxy by

[TABLE]

where $v_{\parallel}$ is the LOS velocity of the galaxy. For simplicity, we set the LOS velocity of the galaxy equal to the LOS velocity of the host dark matter particle, ignoring the $\sim 10\%$ relative velocity bias due to the shear motion between the galaxies and their underlying halos (Guo et al., 2015a). Then we can derive the pair-galaxy bias and the squeezed 3PCF in redshift space using the distorted positions of the galaxies. Several previous papers have also studied the behaviors of redshift space 3PCF in the context of galaxy surveys and simulations (e.g. Gaztañaga & Scoccimarro, 2005; Marín et al., 2008; Guo et al., 2015b).

For the rest of this paper, we calculate the squeezed 3PCF in both real space and redshift space. To avoid confusion, we explicit state in the figure caption whether the figures are made in real space or redshift space.

3.3 Mock LRG and LRG pair catalogue

The first step is to generate mock LRG catalogue from the halo catalogue using a HOD prescription. The HOD we use is a slightly modified version of the the simple HOD fitted to the LRG 2PCF in the Sloan Digital Sky Survey (SDSS; York et al. 2000; Zehavi et al. 2002; Abazajian et al. 2003; Zehavi et al. 2004; Eisenstein et al. 2005a; Zehavi et al. 2005a) and summarized in Zheng et al. (2009); Kwan et al. (2015). The form of our fiducial HOD (equation 34; Kravtsov et al. 2004; Zheng et al. 2005, 2007) is simple in the sense that it only has five parameters, only depends on the halo mass, and ignores effects of halo assembly and halo environment. A brief discussion on more sophisticated halo occupation models is deferred to Section 5.2, but for the purpose of providing an simple application of our methodology, this simple HOD that solely depends on halo mass is sufficient (e.g. Kauffmann et al., 2004; Mo et al., 2004; Blanton et al., 2006; Tinker et al., 2008).

Specifically, this HOD prescription gives the number of centrals and satellites as

[TABLE]

where the parameters are chosen as, $M_{\rm{cut}}\approx 10^{13.35}M_{\odot}$, $M_{1}\approx 10^{13.8}M_{\odot}$, $\sigma=0.85$, $\alpha=1$, and $\kappa=1$ (Zheng et al., 2009; Kwan et al., 2015). The halo mass threshold for satellites is given by $M_{\rm{min}}=\kappa M_{\rm cut}$. This HOD is plotted in Figure 1. From here on, we refer to this HOD as the Z09 design. Also note that we truncate the HOD at $M=4\times 10^{12}M_{\odot}$ to disregard halos with fewer than 100 particles because smaller halos are unreliable and dependent on which halo finder we use. The central is assigned to the centre of mass of the halo with probability $\langle n_{\rm cen}\rangle$. The satellites are assigned to each particle within the halo with probability $p=\langle n_{\rm sat}\rangle/n_{\rm particle}$, where $n_{\rm particle}$ is the total number of particles within the halo. We assign satellite galaxies at particle positions instead of by fitting a smooth profile (e.g. Navarro-Frenk-White profile; Navarro et al. 1996, 1997) since we want to allow satellites to trace halo substructure and avoid assumptions about the isotropy or equilibrium state of the halos.

Looping through all the dark matter particles and assigning galaxies with a small probability introduces significant shot noise, which leads to uncertainties in our $b_{p,g}$ and $Q_{\rm{eff}}$ signals not greater than $10\%$ (refer to Figure 8 and Figure 9). For most of our analysis, the $b_{p,g}$ and $Q_{\rm{eff}}$ signals 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 percent level perturbations to the $b_{p,g}$ and $Q_{\rm{eff}}$ signal due to variations in HOD parameters. To suppress shot noise, we seed the Python pseudo random number generator with the same number at the start of generating the mock galaxy catalogue for each simulation. Thus, when we loop through all the dark matter particles to populate them with galaxies, each dark matter particle is called with the same pseudo random number on each run. This way, not only do we get reproducible results when we re-run the same HOD, but also small changes in the HOD make small changes in the galaxy distribution, thereby sharply reducing the noise in the numerical derivatives of the correlation functions with respect to HOD parameters.

To further suppress the shot noise, we populate each of the 16 simulation boxes with the same HOD 16 different times, each time with a different seed. Then we take the average $Q_{\rm{eff}}$ of the $16^{2}$ realizations. Note that the reseeding is only done for Figure 10 and Figure 13, due to computational limit.

Once we have compiled the mock LRG catalogue from our HOD, the LRG pairs are generated with a pair search on a kd-tree of LRG locations. For our simulation boxes, this HOD typically generates around $5.3\times 10^{5}$ LRGs per box. For maximum pair separation $d_{p,\rm{max}}=10$ Mpc in real space, we generate around $8.7\times 10^{5}$ pairs per box, of which around $2.9\times 10^{5}$ are identified as 1-halo pairs, i.e. the two LRGs forming the pair occupy the same halo. We clarify that all pairs within maximum pair separation $d_{p,\rm{max}}$ are counted in our pair catalogues, e.g. we include all 3 pairs generated by a close triplet. Also for the purpose of this paper, we do not consider effects of fiber collisions or survey boundaries.

Figure 2 shows the distribution of pair separation in all 16 simulation boxes in real space. We see that at $d_{p}<2$ Mpc, 1-halo pairs dominate, while at $d_{p}>2$ Mpc, 2-halo pairs dominate.

In Figure 3, we plot the halo mass corresponding to the LRG 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 the host halo, and for 2-halo pairs, $M$ is the sum of the mass of the two halos. The width of the shaded region shown on the plot is twice the sample standard deviation at that $d_{p}$. We see a distinct peak at around $d_{p}\approx 2$ Mpc, suggesting that close LRG pairs at approximately 2 Mpc apart sample the most massive halos, whereas at separation less or greater than $d_{p}\sim 2$ Mpc, the halo mass of the pair decreases.

Figure 4 shows the halo mass probability density distribution for galaxies (in red) and for pairs (in blue), with the average halo mass denoted by the dashed curves. The pair sample shown on this figure has a maximum pair separation $d_{p,\rm{max}}=10$ Mpc in real space. We see that the galaxies sample halos of approximately $10^{13.5}M_{\odot}$ whereas the pairs sample more massive halos at around $10^{14.2}M_{\odot}$. This plot motivates the point in Section 1 that close pairs are biased tracers of more massive halos and that statistics on the pair field can be used to constrain the high mass regime of the HOD.

3.4 2PCF

By the pair-counting definition of 2PCF, we can calculate the 2PCF as

[TABLE]

where $N_{p}(d_{p})$ is the number of pairs with separation $d_{p}$ and $\Delta V_{\rm{shell}}=4\pi d_{p}^{2}\Delta d_{p}$ is the shell volume. $\bar{n}_{g}$ denotes the LRG number density, and $N_{g}$ denotes the total number of LRGs in the sample.

If we consider RSD, the LOS separation of the pair $d_{\parallel}$ becomes largely meaningless due to velocity dispersion along the LOS. Thus, we calculate the LOS projected 2PCF, $w(d_{\perp})$, as a function of the perpendicular (to LOS) component of pair separation $d_{\perp}$,

[TABLE]

where $d_{\parallel,\rm{max}}=\sqrt{d_{p,\rm{max}}^{2}-d_{\perp}^{2}}$ is the maximum pair separation along the LOS given a maximum pair separation $d_{p,\rm{max}}$. $\sigma_{p}(d_{\perp})$ is the integrated column pair density, and $\Delta V_{\perp}=2\pi\Delta d_{\perp}\times 2d_{\parallel,\rm{max}}$. Unlike the real space 2PCF given by equation 35, the LOS projected 2PCF can be inferred from observed LRG distributions without redshift information. For this reason, we use projected 2PCF instead of the full 2PCF in Section 5.1 as well.

Figure 5 shows $w(d_{\perp})$ calculated from our LRG pair catalogue. We also plot in dashed blue and green curves the projected 2PCFs from the CMASS luminous red galaxies sample of Sloan Digital Sky Survey III (SDSS-III; Eisenstein et al. 2011; Dawson et al. 2013) within the redshift range $0.55<z<0.7$ and $0.43<z<0.55$, respectively. These 2PCFs are directly adapted from Guo et al. (2013). Our 2PCF is roughly consistent with the observed 2PCFs from CMASS red galaxies. The inconsistency in shape can be attributed to a combination of differences in cosmology, problems with halo finders, and oversimplifications in our HOD model itself. However, for the purpose of this methodology paper, the exact form of the HOD is not essential, thus a detailed fit to the observed 2PCF is not important for our analysis. Applications of our methodology to more sophisticated halo occupation models are deferred to future papers.

3.5 Mass assignment function

To be able to utilize FFTs to efficiently compute the power spectra of the LRG and LRG pair fields, we require a mass assignment function that assigns LRGs (and pairs) onto a regular 3D grid in such a way that generates a smooth LRG (and pair) density field. For the purpose of this paper, we assign LRGs and pairs onto an evenly spaced $512^{3}$ grid spanning the full simulation box using the Triangular Shaped Cloud (TSC) method (Jeong, 2010). For a simulation box size $L_{\rm{box}}=1100h^{-1}$ Mpc with 5123, the cell size is $\Delta x\approx 2.15h^{-1}$ Mpc. Thus, a TSC of base length $2\Delta x$ spans 27 cells on the grid, with 3 cells along each dimension. The fact that each TSC can span at most 3 cells along each direction becomes important in the next section as we design our custom weighting functions.

3.6 Weighting function

In Section 2, we present our method of efficiently computing the cross-correlation between galaxy pairs and galaxies, from which we infer the squeezed 3PCF. This routine requires 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 considered as the third galaxy when performing the cross-correlation. This means that the weighting function should have a hollow spherical cavity in the centre. The radius of the cavity should be $r_{\rm{cut}}\approx d_{p,\rm{max}}/2+2\Delta x$, where $d_{p,\rm{max}}$ is the maximum pair separation and $\Delta x$ is the cell size. The extra $2\Delta x$ comes from avoiding the TSC of the pair, which can extend a maximum of 2 cells outwards from either galaxy along the axis of the pair.

(b) The weighting function should minimize the FoG effect when we apply our methodology to redshift space data, where the correlation along the LOS becomes largely meaningless. To disregard correlation between pairs and galaxies along the LOS, we “drill" a cylindrical hole through the centre of the weighting function parallel to the LOS. The radius of the cylindrical hole should be $r_{\rm{cut}}\approx d_{p,\rm{max}}/2+2\Delta x$ for the same reason as explained in (a).

(c) The weighting function should sample scales approximately in the range $k=0.01\sim 0.1h$ Mpc*-1*. The scale range is limited on the large scale by the simulation box size, where we cannot meaningfully interpret scales larger than the box size, which translates to a minimum frequency of $k_{\rm{min}}=2\pi/L_{\rm{box}}\approx 0.006h$ Mpc*-1*. The scale range is limited on the small scale by self-counting, so we do not want the third galaxy to be too close to the pair that their TSCs overlap. We also avoid small scale contribution since we want to be computing bias in the linear regime.

In addition to these requirements, we also want the weighting function to be mostly smooth and analytic for ease of performing FFTs and to avoid high $k$ modes in the weighting function. Specifically, our weighting function is given in real space by,

[TABLE]

where $\mathcal{N}$ is a normalization constant. $r_{\perp}$ and $r_{\parallel}$ are distance between the third galaxy and the centre of the pair, projected perpendicular and parallel to the LOS respectively. The function is parameterized by $\sigma$, $r_{s}$, $s_{0}$ and $r_{\rm{cut}}$, which are determined by the maximum pair separtion$d_{p,\rm{max}}$ that we choose. Although the function is divided into two segment and thus not fully analytical, such engineering is necessary to remove contamination from pair self counting. We choose the parameters such that the jump in $W(r_{\perp},r_{\parallel})$ at $r_{\rm{cut}}$ is small (refer to Figure 6).

For the rest of the analysis in this paper, we use $d_{p,\rm{max}}=10$ Mpc when working in real space and $d_{p,\rm{max}}=30$ Mpc when working in redshift space. We choose $d_{p,\rm{max}}=30$ Mpc in redshift space because the typical velocity dispersion along the LOS translates to $\sim 10$ Mpc displacement, so we extend the pair search in both direction by $10$ Mpc to capture most of the pairs. The parameters of the weighting function for both cases are summarized in Table 1.

Figure 6 shows a 2D cutaway of the weighting functions at $r_{\parallel}=0$ in the positive half of $r_{\perp}$. The solid red curve shows the weighting function for $d_{p,\rm{max}}=10$ Mpc and the dashed blue curve shows the weighting function for $d_{p,\rm{max}}=30$ Mpc. The dotted black lines denote the location of $r_{\rm{cut}}$ to show that the jump at $r_{\rm{cut}}$ is small.

According to equation 11, the weighting function is multiplied with the pair density field in Fourier space as a sampling weight for different k modes. The Fourier Transforms of the weighting functions are plotted in Figure 7, where the top panel shows the weighting function for $d_{p,\rm{max}}=10$ Mpc whereas the middle panel shows the weighting function for $d_{p,\rm{max}}=30$ Mpc. The bottom panel shows cutaways of the two cases at $k_{\parallel}=0$. The colorbar gives $\tilde{W}\ k_{\perp}/k$, which is proportional to the number of modes sampled by the weighting function in Fourier space. We see that both weighting functions are indeed heavily sampling scales between $0.01\sim 0.1h$ Mpc*-1*. The first and second order ringing seen at higher frequency are due to the steep drop-off in the real space weighting function towards small $r_{\perp}$. However, the ringings have much lower integrated weight in $k$-space and are not expected to affect our analysis significantly.

Now that we have defined the weighting functions, we can calculate $d_{w}$ (equation 30) and $\xi(d_{w})$ (equation 29), which we need to compute the squeezed 3PCF coefficient $Q_{\rm{eff}}$ from the pair-galaxy bias $b_{p,g}$. To do so, we employ the 2PCF given by Figure 2 and Figure 3 in Eisenstein et al. (2005b), which is constructed from the linear power spectra calculated using CMBFAST (Seljak & Zaldarriaga, 1996; Zaldarriaga et al., 1998; Zaldarriaga & Seljak, 2000). For any of our weighting functions, we fit the 2PCF with a power law in the scale range sampled by the weighting function. Then we use the best fit exponent $p$ to calculate $d_{w}$ and $\xi(d_{w})$. Specifically, the $d_{w}$ and $\xi(d_{w})$ for our $d_{p,\rm{max}}=10$ Mpc and 30 Mpc are given in the last two columns of Table 1. It should be pointed out that the analysis in Eisenstein et al. (2005b) uses a slightly different cosmology, where $h=0.7$, $\Omega_{m}h^{2}=0.13$ and $\Omega_{b}h^{2}=0.024$. However, the correction is small and does not meaningfully change our calculations.

Now that we have defined our mass assignment function and our weighting functions, we apply the TSC to the LRGs and LRG pairs in all 16 simulation boxes to generate smooth density fields of the LRGs and pairs on our grid. The resulting density fields are then convolved with the appropriate weighting function to compute the pair-galaxy bias $b_{p,g}$ (equation 11) and the reduced 3PCF $Q_{\rm{eff}}$ (equation 32).

4 Results and Discussion

In this section, we present the results of implementing our efficient method of calculating $b_{p,g}$ and $Q_{\rm{eff}}$ using the N-body simulation and the Z09 HOD. The results shown are averaged over the 16 cosmological boxes to reduce noise. Section 4.1 presents the results in real space, and Section 4.2 presents the results in redshift space.

4.1 Results in real space

We average the $b_{p,g}$ calculated from all 16 simulation boxes at $z=0.5$ and plot the averaged real space $b_{p,g}$ in the top panel of Figure 8. We have chosen maximum pair separation $d_{p,\rm{max}}=10$ Mpc and a Gauss-Sigmoid weighting function with the following parameters: $\sigma=50$ Mpc, $r_{s}=30$ Mpc, $s_{0}=0.5$ Mpc*-1* and $r_{\rm{cut}}=15$ Mpc (Refer to Section 3.6 for details). The associated uncertainty map is given in middle panel. The bottom panel shows the reduced 3PCF $Q_{\rm eff}$ in real space, as calculated from $b_{p,g}$ using equation 32.

The top panel of Figure 8 shows that the pair-galaxy bias is close to 1.3 at pair separation less than $\sim 1$ Mpc, but dramatically increases by a factor of 2 and peaks at $\sim 2$ Mpc, before decreasing again towards greater pair separation. The bottom panel shows that the squeezed 3PCF coefficient $Q_{\rm{eff}}$ recovers the same behavior. We can analytically estimate the value of $b_{p,g}$ by computing the halo bias as a function of halo mass and then take the ratio of the halo bias averaged over all pairs and the halo bias averaged over all galaxies (using the halo mass function of the pairs and galaxies shown in Figure 4). Our estimate yields $b_{p,g}\approx 1.7$, which is consistent with the values we see in the top panel of Figure 8.

Referring to Figure 2 and Figure 3, we see that pairs with separation $<1$ Mpc are mostly 1-halo pairs with lower halo mass. $b_{p,g}$ is approximately 1.3 at $d_{p}<1$ Mpc. The fact that the pair-galaxy bias $b_{p,g}$ is greater than 1 shows that close galaxy pairs are more biased than individual galaxies. Pairs with separation $\sim 2$ Mpc consist of $\sim 70\%$ larger 1-halo pairs and $\sim 30\%$ 2-halo pairs. Figure 3 shows that these pairs trace the most massive halos, thus giving them the highest bias.

Figure 8 also shows a strong anisotropic signal in $b_{p,g}$ and $Q_{\rm{eff}}$ at large pair separation $d_{p}>4$ Mpc. The anisotropy is due to the fact that the weighting function samples favorably volume in the transverse direction and avoids volumes along the LOS direction. When a pair is oriented transverse to the LOS, it establishes the transverse direction as the preferred direction for the large-scale filament. Thus, when a pair is perpendicular to the LOS, the weighting function overlaps the most with the large-scale filament, giving rise to much stronger clustering. Conversely, when a pair is oriented along the direction of LOS, the filament is preferentially aligned with the LOS, which does not overlap significantly with the weighting function, thus giving rise to weaker clustering. The large $b_{p,g}$ value at $d_{p}>4$ Mpc suggests that the filament, not just the halo, is also strongly biased against the galaxy field.

The middle panel of Figure 8 shows that most of the $b_{p,g}$ values have percent level uncertainties, except the ones towards the top left corner of the plot, where the small sample size contributes to the high noise.

The bottom panel of Figure 8 shows that the reduced 3PCF recovers the same behavior as seen in the pair-galaxy bias $b_{p,g}$. This is consistent with our conclusion that in the strictly squeezed limit, the 3PCF would basically reduce to $Q_{\rm{eff}}=b_{p,g}/2$. At pair separation $d_{p}>5$ Mpc, we continue to see the anisotropy in bias due to the filamentary structure aligning with the preferred direction of the weighting function. However, we also see suppression of a few percent in $Q_{\rm eff}$ due to the $\xi(d_{w})$ term. The associated uncertainty map of $Q_{\rm{eff}}$ is omitted, but shares the same properties as the uncertainty in $b_{p,g}$, but with approximately half the amplitude.

4.2 Results in redshift space

Similar to the last section, we can incorporate RSD and calculate the $b_{p,g}$ and $Q_{\rm{eff}}$ in redshift space, as shown in Figure 9. The top panel shows the the pair-galaxy bias $b_{p,g}$ as a function of pair separation. The middle panel shows the corresponding uncertainty on $b_{p,g}$. The bottom panel shows the reduced 3PCF $Q_{\rm{eff}}$ in redshift space. For these plots, the weighting function is chosen to be the $d_{p,\rm{max}}=30$ Mpc. Figure 10 shows a zoom-in version of the $Q_{\rm{eff}}$ signal in redshift space with finer bins on pair separation to be comparable to the real space $Q_{\rm{eff}}$ signal plotted in the bottom panel of Figure 8. To compensate for the loss in signal-to-noise by choosing finer bins, we use 16 different seeds on top of 16 different simulation boxes, as discussed in Section 3.3. One should disregard the high values along the outer circular edge of plots as these values are calculated from small sample of pairs and are thus dominated by noise.

The top and bottom panel of Figure 9 show a vertical “ridge" at $d_{\perp}\approx 2$ Mpc in $b_{p,g}$ and $Q_{\rm{eff}}$ that corresponds to the $\sim 2$ Mpc peak in the real space $b_{p,g}$ and $Q_{\rm{eff}}$ signal (Figure 8). This ridge tracks the halo mass peak at $d_{p}\sim 2$ Mpc in Figure 3, but smeared out along the LOS by FoG. The strong clustering along the ridge all the way to $d_{\parallel}\sim 30$ Mpc suggests that the galaxy pair-wise peculiar velocities reach $\sim 2000$ km/s in the most massive halos. The zoom-in in Figure 10 shows that the ridge is more clustered at larger LOS separation, which is due to the fact that larger LOS pair separation corresponds to larger peculiar velocities of the galaxies, which in turn corresponds to more massive halos. We also point out that the infall velocity of merging dark matter halos may also contribute to this ridge. However, we do not expect the infall velocity to frequently reach $\sim 2000$ km/s and speculate that the strongest $Q_{\rm{eff}}$ signal towards the top of the ridge should be completely dominated by FoG.

We also see weaker anisotropy at large pair separation in redshift space. This is due to the RSD smearing out the filamentary structures, thus also randomizing the preferred direction given an LRG pair. The rise in bias value towards larger pair separation perpendicular to the LOS can be explained by the fact that more distant pairs along the same filament are better at picking out the preferred direction of the filament, which gives rise to higher bias given the anisotropy of the weighting function.

Both the real space and redshift space $b_{p,g}$ signal show a factor of 2 variation. This suggests that the hierarchical ansatz on the 3PCF does not hold in the squeezed limit. However, as we have discussed in this section, the variations in the squeezed 3PCF signal provides rich information on the galaxy-halo connection on the small scale and on the filamentary structure on the large scale. In the following section, we explore using the squeezed 3PCF as a diagnostics on the galaxy-halo connection.

5 HOD dependency of the 3PCF

We next discuss the relationship between the squeezed 3PCF and the HOD. Specifically, we explore the possibility of using the squeezed 3PCF to break degeneracies in the 2PCF and constrain the high mass end of the HOD.

5.1 Controlling the 2PCF with the HOD

To investigate how the squeezed 3PCF breaks degeneracies in the 2PCF, we first need to understand how the HOD controls the 2PCF and devise HOD recipes that would change the 3PCF while maintaining the same 2PCF. First we attempt a linear emulator to model the 2PCF as a function of HOD parameters. As discussed in Section 3.4, we study the projected 2PCF $w(d_{\perp})$ as given by equation 36. To the first order, we can expand the projected 2PCF as

[TABLE]

where $\textbf{p}=\left[p_{1},p_{2},p_{3},p_{4}\right]=\left[\log_{10}M_{\rm{cut}},\log_{10}M_{1},\sigma,\alpha\right]$ are the HOD parameters that we use to control the 2PCF. Note that we omitted parameter $\kappa$ because our tests show that $\kappa$ perturbs the 2PCF insignificantly for the variation amplitude that we are interested in. $\textbf{p}_{0}$ denotes the Z09 design of the HOD, as quoted in Section 3.3. $w_{i}(\bf{p})$ denotes the value of the projected 2PCF integrated over the $i$-th bin in $d_{\perp}$. By choosing an uniform set of bins on the projected 2PCF, we form a vector $\left[w_{i}(\bf{p})\right]$ that holds the shape information of the projected 2PCF. $\left[w_{i,p_{j}}\right]$ records the matrix of the derivatives of the 2PCF integrated over the $i$-th bin, $w_{i}(\bf{p})$, against the $j$-th HOD parameter, $p_{j}$. We can determine each derivative $w_{i,p_{j}}$ by independently varying each parameter of the HOD and tracking the variation incurred onto $w_{i}(\bf{p})$. Once we have $\left[w_{i,p_{j}}\right]$, we solve for the set of HOD parameters $\bf{p}$ that minimizes $|\left[w_{i}(\bf{p})\right]-\left[w_{i}(\bf{p}_{0})\right]|$. However, it turns out that the linear model does not produce a good emulator for the projected 2PCF. Our tests show that the linear approximation breaks down when varying $\alpha$ by more than $\sim 1\%$.

Thus, we develop an emulator that incorporates second-order derivatives. To formulate a second order emulator of $\left[w_{i}(\bf{p})\right]$, we expand equation 38 to the 2nd order

[TABLE]

where $w_{i}$ corresponds to the integrated value of the projected 2PCF over the $i$-th bin. $w_{i,p_{j}}$ denotes the first derivatives while $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 of 14 different and independent evaluations of $\left[w_{i}(\bf{p})\right]$ to solve for all 14 derivatives.

Once we have solved for the first and second derivatives on the projected 2PCF, we have determined the full parametrized form of our $2^{\rm{nd}}$ order emulator for $\left[w_{i}(\bf{p})\right]$. Then we vary the HOD parameters by first varying $\alpha$ (we choose $\alpha$ because it controls the number of satellites) by $+5\%$ and $-5\%$ and then perturbing the remaining HOD parameters to find the $\left[w_{i}(\bf{p})\right]$ that best emulates $\left[w_{i}(\bf{p}_{0})\right]$. Specifically, we minimize $\chi^{2}$ with an error on $\delta w_{i}$ that scales with inverse sample volume. Considering that our full sample live in a sphere of diameter $d_{p,\rm{max}}$ in pair separation space, the error is given by $\sigma_{i}^{2}=\sqrt{d_{p,\rm{max}}^{2}-d_{\perp,i}^{2}}/d_{\perp,i}$ for an even binning in $d_{\perp,i}$.

For both $\delta\alpha=0.05$ and $\delta\alpha=-0.05$, we find local minima of $\chi^{2}$ from an arbitrarily selected set of initial locations in the HOD parameter space. Then we visually select the set of HOD parameters that best emulates the Z09 2PCF. Table 2 shows the resulting parameters of these best fit HODs. The minimization is done using the Powell algorithm (Powell, 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 Z09 design in Figure 11. The number of centrals and satellites are plotted in dotted red curve. The HOD corresponding to $\alpha=0.95$ is plotted in blue and the HOD corresponding to $\alpha=1.05$ is plotted in magenta. The solid curves show $n_{\rm{cent}}$ and the dashed curves show $n_{\rm{sat}}$. The bottom panel shows the absolute difference between the $\alpha=1.05$ HOD and the $\alpha=0.95$ HOD, where again the solid line shows the difference in $n_{\rm{cent}}$ and the dashed line shows the difference in $n_{\rm{sat}}$. We see that by increasing $\alpha$, we are increasing the number of satellites at the high-mass end of the HOD as expected. To compensate for such changes so as to maintain the same 2PCF, we see a decrease in the number of centrals at halo mass $M<10^{13}M_{\odot}$. Thus, while we are increasing the number of LRG pairs in the largest halos, our algorithm is removing central-central pairs in the more common but smaller halos to compensate.

Figure 12 top panel shows the projected 2PCF corresponding to the Z09 design. The bottom panel shows the relative change in the projected 2PCFs after perturbing HOD parameters. The dashed curves give the relative change in projected 2PCF after only changing $\alpha$, with blue corresponding to $\alpha=0.95$ and magenta corresponding to $\alpha=1.05$. The solid curves show the relative change in projected 2PCF after also perturbing the other HOD parameters to suppress the deviations in the projected 2PCF from that of the Z09 design. We see that for a $10\%$ change in $\alpha$, the best-fit 2PCFs shows an offset of $\sim 3\%$ at most, with the largest fractional deviations occurring between $d_{\perp}\approx 2$ Mpc. The fractional deviation is below $1\%$ beyond 4 Mpc. This change is considerably smaller than the $\sim 16\%$ change that occurs in the projected 2PCF before we adjust the other HOD parameters.

5.2 3PCF dependence on the HOD

Now we compute the change in the reduced 3PCF $Q_{\rm{eff}}$ corresponding to the two best-fit HODs that we found in the last section.

We plot the relative change in $Q_{\rm{eff}}$ in Figure 13, where $\delta Q_{\rm{eff}}$ is defined as

[TABLE]

The $\delta Q_{\rm{eff}}$ signal is averaged over 16 simulation boxes, each with 16 different seeds to suppress shot noise. We see a relative change in $Q_{\rm{eff}}$ not greater than $4\%$ at $d_{\perp}\sim$ 1–3 Mpc, with the relative change increasing towards large $d_{\parallel}$. This suggests that $Q_{\rm{eff}}$ is most sensitive to $\delta\alpha$ along the 2 Mpc ridge (refer to the bottom panel of Figure 9) that corresponds to the most massive halos. This makes sense since $\alpha$ controls the number of galaxies in massive halos.

In detail, the morphology of $\delta Q_{\rm{eff}}$ shown in Figure 13 recalls the variations in the projected 2PCF (Figure 12): an excess of $3\%$ at 2 Mpc scale, dominating at large $d_{\parallel}$ where the most massive halo pairs occur, followed by $<1\%$ differences at larger $d_{\perp}$. We therefore wanted to check how the $Q_{\rm{eff}}$ variations in Figure 13 compared to the variations in the full anisotropic 2PCF. This is shown in Figure 14, where we have computed $\xi(d_{\perp},d_{\parallel})$ for both the $\alpha=0.95$ and $\alpha=1.05$ HODs and then shown the fractional difference. One sees that while the broadest patterns coincide, the details vary significantly. Specifically, the maximum variation in $\xi$ is $\sim 6\%$, twice as large as that of $Q_{\rm{eff}}$. We see a strong signal in $\delta Q_{\rm{eff}}$ at the smallest $d_{\perp}$ but not in $\delta\xi$. The ridge in $\delta Q_{\rm{eff}}$ extends all the way down to the smallest $d_{\parallel}$ whereas the $\delta\xi$ signal only extends down to about $d_{\parallel}\sim 6$ Mpc. There is also a consistently negative region in $\delta\xi$ at $d_{\perp}\sim$ 5–15 Mpc, which is not reflected in the $\delta Q_{\rm{eff}}$ signal. The angular gradient at large pair-separation is also different between $\delta Q_{\rm{eff}}$ and $\delta\xi$. Thus, we conclude that the reduced 3PCF $Q_{\rm{eff}}(d_{\perp},d_{\parallel})$ is not a simple rescaling of the anisotropic 2PCF $\xi(d_{\perp},d_{\parallel})$. However, as the variations in the reduced 3PCF $Q_{\rm{eff}}$ ($<4\%$) are similar in size and scale to those in the projected 2PCF, it is unclear whether the squeezed 3PCF is breaking a degeneracy within this particular space of HOD models.

While this test was inconclusive within this model space, we expect that there could be reasonable extensions of the HOD where the variation to $Q_{\rm{eff}}$ is disentangled from the variation to the projected 2PCF. For these scenarios, $Q_{\rm{eff}}$ would be a new diagnostic to distinguish among different HOD models. However, it is also possible that the squeezed 3PCF is well predicted by the 2PCF for certain classes of HOD models. Such a predicative relationship would be interesting in itself, as any discrepancy between the observed squeezed 3PCF and the 2PCF could then be used to falsify such classes of HOD models.

In future works, we intend to apply our methodology to more sophisticated halo occupation models that incorporate dependence on more parameters than the one we used for this paper, such as models that incorporate halo assembly bias (Wechsler et al., 2002; Gao et al., 2005; Wechsler et al., 2006; Croton et al., 2007; Zentner et al., 2014; Lehmann et al., 2017). In particular, Hearin et al. (2016) parameterize the assembly bias in an HOD-like model (“Decorated HOD"). Recently, there has also been a series of studies exploring the dependence of the HOD on halo environment (e.g. Mandelbaum et al., 2006; Gil-Marín et al., 2011; Croft et al., 2012; Tonnesen & Cen, 2015) and halo geometry (van Daalen et al., 2012). It is plausible that the pair-galaxy bias and the squeezed 3PCF can offer extra constraints for these sophisticated halo occupation models.

6 Conclusions

To summarize, this paper presents the methodology of cross-correlating the galaxy pair field with the galaxy field as a squeezed 3PCF. We develop a fast and efficient method using FFTs to compute the pair-galaxy bias $b_{p,g}$ and its associated reduced squeezed 3PCF $Q_{\rm{eff}}$. We implement our method using a series of N-body cosmology simulations and show the $b_{p,g}$ and $Q_{\rm{eff}}$ as functions of pair separation in real space and in redshift space. In real space, a significant peak in $b_{p,g}$ and $Q_{\rm{eff}}$ at pair separation $\sim$2 Mpc is observed, which is due to the fact that galaxy pairs at $\sim$2 Mpc separation trace the most massive dark matter halos. We also see strong anisotropy in the $b_{p,g}$ and $Q_{\rm{eff}}$ signals that tracks the large-scale filamentary structure.

In redshift space, the FoG effect smears out the 2 Mpc peak into a ridge along the LOS. The anisotropy signal is also smeared out along the LOS, but persists at large pair separation. Both in real space and redshift, we see a factor of 2 variation in the squeezed 3PCF, contradicting the hierarchical ansatz. However, the variations in the squeezed 3PCF offers rich information on the galaxy-halo connection and the large-scale filamentary structure.

Thus, we investigate using the reduced squeezed 3PCF $Q_{\rm{eff}}$ to constrain the high mass end of the HOD while breaking degeneracies in the 2PCF. We develop a second order emulator to model the projected 2PCF around the Z09 design HOD. Using the emulator, we find two simple HOD prescriptions that constrain the variations in the projected 2PCF to $<3\%$ when handed a $10\%$ perturbation to HOD parameter $\alpha$. The resulting perturbation to the squeezed 3PCF is no greater than $4\%$, which cannot be immediately disentangled from the variations in the 2PCF. Nevertheless, we propose that more sophisticated HOD models or different variations to the HOD could further disentangle the squeezed 3PCF from the 2PCF.

In future papers, we intend to construct a more sophisticated emulator of the 2PCF as a function of HOD parameters. Another key objective is to extend our methodology to more sophisticated halo occupation models, such as ones that incorporate assembly bias, and explore how the pair-galaxy bias and the squeezed 3PCF can be used to distinguish among these more complex models. Eventually we want to apply our methodology to observations to detect these signals.

Acknowledgements

We would also like to thank Dr. Zachary Slepian and Joshua Speagle for useful discussions. DJE and LG are supported by National Science Foundation grant AST-1313285. DJE is additionally supported by U.S. Department of Energy grant DE-SC0013718 and as a Simons Foundation Investigator.