How to Optimally Constrain Galaxy Assembly Bias: Supplement Projected Correlation Functions with Count-in-cells Statistics

TL;DR

This paper investigates how combining traditional galaxy clustering and lensing measurements with count-in-cells statistics can more effectively constrain galaxy assembly bias, improving understanding of galaxy-halo connections.

Contribution

It demonstrates that count-in-cells statistics provide more efficient constraints on galaxy assembly bias than traditional methods when combined with correlation functions.

Findings

Count-in-cells statistics outperform traditional observables in constraining assembly bias.

Combining $w_p(r_p)$ with count statistics reduces parameter degeneracies.

Count statistics should be integrated into future galaxy-halo connection studies.

Abstract

Most models for the connection between galaxies and their haloes ignore the possibility that galaxy properties may be correlated with halo properties other than mass, a phenomenon known as galaxy assembly bias. Yet, it is known that such correlations can lead to systematic errors in the interpretation of survey data. At present, the degree to which galaxy assembly bias may be present in the real Universe, and the best strategies for constraining it remain uncertain. We study the ability of several observables to constrain galaxy assembly bias from redshift survey data using the decorated halo occupation distribution (dHOD), an empirical model of the galaxy--halo connection that incorporates assembly bias. We cover an expansive set of observables, including the projected two-point correlation function $w_{\mathrm{p}}(r_{\mathrm{p}})$, the galaxy--galaxy lensing signal $\Delta \Sigma(r_{\mathrm{p}})$, the void probability function $\mathrm{VPF}(r)$, the distributions of counts-in-cylinders $P(N_{\mathrm{CIC}})$, and counts-in-annuli $P(N_{\mathrm{CIA}})$, and the distribution of the ratio of counts in cylinders of different sizes $P(N_2/N_5)$. We find that despite the frequent use of the combination $w_{\mathrm{p}}(r_{\mathrm{p}})+\Delta \Sigma(r_{\mathrm{p}})$ in interpreting galaxy data, the count statistics, $P(N_{\mathrm{CIC}})$ and $P(N_{\mathrm{CIA}})$, are generally more efficient in constraining galaxy assembly bias when combined with $w_{\mathrm{p}}(r_{\mathrm{p}})$. Constraints based upon $w_{\mathrm{p}}(r_{\mathrm{p}})$ and $\Delta \Sigma(r_{\mathrm{p}})$ share common degeneracy directions in the parameter space, while combinations of $w_{\mathrm{p}}(r_{\mathrm{p}})$ with the count statistics are more complementary. Therefore, we strongly suggest that count statistics should be used to complement the canonical observables in future studies of the galaxy--halo connection.

Tables

Table 1. Table 1: In this table, we list the fiducial HOD parameters adopted for each luminosity threshold of galaxies, taken from the fits of Zentner et al. ( 2019 ) . Of the 5 standard HOD parameters, M min subscript 𝑀 min M_{\mathrm{min}} is the mass at which a halo has a 50% probability of hosting a central galaxy; σ log ⁡ M subscript 𝜎 𝑀 \sigma_{\log M} determines the rate that ⟨ N cen | M vir ⟩ inner-product subscript 𝑁 cen subscript 𝑀 vir \langle N_{\mathrm{cen}}|M_{\mathrm{vir}}\rangle transitions from zero to unity; M 0 subscript 𝑀 0 M_{0} is the truncating mass, below which no satellite galaxies are allowed; M 1 subscript 𝑀 1 M_{1} is the halo mass at which the mean satellite number is unity; and α 𝛼 \alpha is the index of the satellite occupation power law. Besides the standard HOD parameters, we also allow A cen subscript 𝐴 cen A_{\mathrm{cen}} and A sat subscript 𝐴 sat A_{\mathrm{sat}} to vary, which control the amount of galaxy assembly bias for central and satellite galaxies respectively. In doing this we treat galaxy assembly bias as a deviation from the standard HOD model to be constrained.

Fiducial HOD Parameters
	$\log M_{\text{min}}$	${\sigma }_{\log M}$	$\alpha$	$\log M_{\text{0}}$	$\log M_{\text{1}}$	$A_{\text{cen}}$	$A_{\text{sat}}$
$M_{\text{r}}$ -19.0	11.64	0.5119	1.040	10.25	12.80	0	0
$M_{\text{r}}$ -19.5	11.75	0.4458	1.116	11.29	13.06	0	0
$M_{\text{r}}$ -20.0	11.97	0.3485	1.144	11.31	13.29	0	0
$M_{\text{r}}$ -20.5	12.25	0.1854	1.197	11.20	13.59	0	0
$M_{\text{r}}$ -21.0	12.82	0.5595	1.337	11.96	13.99	0	0

Table 2. Table 2: Definition of bins for each observable. We show the values of bin centres for the bins in which w p ​ ( r p ) subscript 𝑤 p subscript 𝑟 p w_{\mathrm{p}}(r_{\mathrm{p}}) and Δ ​ Σ ​ ( r p ) Δ Σ subscript 𝑟 p \Delta\Sigma(r_{\mathrm{p}}) are measured, the set of radii of spheres used for evaluating VPF ​ ( r ) VPF 𝑟 \mathrm{VPF}(r) , the intervals defining each bin in the histograms of counts-in-cylinders and annuli, and the ratio N 2 / N 5 subscript 𝑁 2 subscript 𝑁 5 N_{2}/N_{5} . We use the same number of bins for each observable, such that our comparison of constraining power is not sensitive to bin number. * Note that the number density n gal subscript 𝑛 gal n_{\mathrm{gal}} is listed as the first bin of w p ​ ( r p ) subscript 𝑤 p subscript 𝑟 p w_{\mathrm{p}}(r_{\mathrm{p}}) in this table, but it is included in the analysis for all possible combinations of observables.

Observable Bin Definition
Observable	$w_p\left(r_p\right)$	$\Delta \Sigma \left(r_p\right)$	$\mathrm{VPF}\left(r\right)$	$P\left(N_{\mathrm{CIC}}\right)$	$P\left(N_{\mathrm{CIA}}\right)$	$P\left(N_2/N_5\right)$
Bin No.	$r_p\left[h^{-1}\mathrm{Mpc}\right]$	$r_p\left[h^{-1}\mathrm{Mpc}\right]$	$r\left[h^{-1}\mathrm{Mpc}\right]$	$N_{\mathrm{CIC}}$	$N_{\mathrm{CIA}}$	$N_2/N_5$
1	*	0.11	1.00	$\left\{0\right\}$	$\left\{0\right\}$	$[0.000,0.033)$
2	0.11	0.13	1.08	$\left\{1\right\}$	$\left\{1\right\}$	$[0.033,0.067)$
3	0.14	0.16	1.17	$\left\{2\right\}$	$\left\{2\right\}$	$[0.067,0.100)$
4	0.17	0.20	1.27	$\left\{3\right\}$	$\left\{3\right\}$	$[0.100,0.133)$
5	0.20	0.24	1.37	$\left\{4\right\}$	$\left\{4\right\}$	$[0.133,0.167)$
6	0.25	0.29	1.49	$\left\{5\right\}$	$\left\{5\right\}$	$[0.167,0.200)$
7	0.30	0.35	1.61	$\left\{6\right\}$	$\left\{6\right\}$	$[0.200,0.233)$
8	0.37	0.43	1.74	$\left\{7\right\}$	$\left\{7\right\}$	$[0.233,0.267)$
9	0.45	0.52	1.89	$\left\{8\right\}$	$\left\{8\right\}$	$[0.267,0.300)$
10	0.54	0.62	2.04	$\left\{9\right\}$	$\left\{9\right\}$	$[0.300,0.333)$
11	0.66	0.76	2.21	$[10,12)$	$[10,12)$	$[0.333,0.367)$
12	0.81	0.92	2.40	$[12,13)$	$[12,14)$	$[0.367,0.400)$
13	0.99	1.11	2.59	$[13,15)$	$[14,16)$	$[0.400,0.433)$
14	1.20	1.35	2.81	$[15,18)$	$[16,19)$	$[0.433,0.467)$
15	1.47	1.63	3.04	$[18,20)$	$[19,22)$	$[0.467,0.500)$
16	1.79	1.98	3.29	$[20,24)$	$[22,26)$	$[0.500,0.533)$
17	2.18	2.39	3.56	$[24,27)$	$[26,30)$	$[0.533,0.567)$
18	2.66	2.90	3.86	$[27,31)$	$[30,35)$	$[0.567,0.600)$
19	3.24	3.51	4.18	$[31,36)$	$[35,41)$	$[0.600,0.633)$
20	3.95	4.26	4.52	$[36,42)$	$[41,48)$	$[0.633,0.667)$
21	4.82	5.16	4.89	$[42,48)$	$[48,57)$	$[0.667,0.700)$
22	5.88	6.25	5.30	$[48,55)$	$[57,66)$	$[0.700,0.733)$
23	7.17	7.57	5.74	$[55,64)$	$[66,78)$	$[0.733,0.767)$
24	8.75	9.17	6.21	$[64,74)$	$[78,91)$	$[0.767,0.800)$
25	10.6	11.11	6.72	$[74,85)$	$[91,106)$	$[0.800,0.833)$
26	13.0	13.46	7.28	$[85,98)$	$[106,125)$	$[0.833,0.867)$
27	15.8	16.30	7.89	$[98,113)$	$[125,146)$	$[0.867,0.900)$
28	19.3	19.75	8.53	$[113,130)$	$[146,171)$	$[0.900,0.933)$
29	23.6	23.93	9.24	$[130,150)$	$[171,200)$	$[0.933,0.967)$
30	28.7	28.99	10.00	$[150,+\infty )$	$[200,+\infty )$	$[0.967,1.000]$

Table 3. Table 3: M r < − 19.0 subscript 𝑀 𝑟 19.0 M_{r}<-19.0

Observable combination	$\sigma \left(\log M_{\min }\right)$	$\sigma \left({\sigma }_{\log M}\right)$	$\sigma \left(\alpha \right)$	$\sigma \left(\log M_1\right)$	$\sigma \left(A_{\mathrm{cen}}\right)$	$\sigma \left(A_{\mathrm{sat}}\right)$
$w_p\left(r_p\right)$	0.30	0.645	0.038	0.05	0.71	0.34
$\Delta \Sigma \left(r_p\right)$	0.10	0.157	0.066	0.09	0.79	0.32
$\mathrm{VPF}\left(r\right)$	0.22	0.129	0.689	0.37	0.63	3.99
$P\left(N_{\mathrm{CIC}}\right)$	0.18	0.179	0.091	0.18	1.34	0.20
$P\left(N_{\mathrm{CIA}}\right)$	0.06	0.151	0.049	0.07	0.46	0.36
$P\left(N_2/N_5\right)$	0.05	0.073	0.040	0.06	0.57	0.21
$w_p\left(r_p\right)+\Delta \Sigma \left(r_p\right)$	0.06	0.121	0.016	0.02	0.41	0.20
$w_p\left(r_p\right)+\mathrm{VPF}\left(r\right)$	0.05	0.086	0.019	0.03	0.31	0.25
$w_p\left(r_p\right)+P\left(N_{\mathrm{CIC}}\right)$	0.03	0.048	0.013	0.02	0.15	0.15
$w_p\left(r_p\right)+P\left(N_{\mathrm{CIA}}\right)$	0.05	0.111	0.014	0.02	0.16	0.19
$w_p\left(r_p\right)+P\left(N_2/N_5\right)$	0.03	0.062	0.013	0.02	0.27	0.15
$\Delta \Sigma \left(r_p\right)+\mathrm{VPF}\left(r\right)$	0.06	0.084	0.036	0.04	0.31	0.19
$\Delta \Sigma \left(r_p\right)+P\left(N_{\mathrm{CIC}}\right)$	0.04	0.054	0.022	0.03	0.19	0.12
$\Delta \Sigma \left(r_p\right)+P\left(N_{\mathrm{CIA}}\right)$	0.05	0.082	0.030	0.03	0.17	0.19
$\Delta \Sigma \left(r_p\right)+P\left(N_2/N_5\right)$	0.04	0.064	0.022	0.03	0.29	0.14
$\mathrm{VPF}\left(r\right)+P\left(N_{\mathrm{CIC}}\right)$	0.05	0.057	0.030	0.06	0.43	0.19
$\mathrm{VPF}\left(r\right)+P\left(N_{\mathrm{CIA}}\right)$	0.05	0.061	0.046	0.04	0.23	0.35
$\mathrm{VPF}\left(r\right)+P\left(N_2/N_5\right)$	0.03	0.050	0.030	0.04	0.38	0.20
$P\left(N_{\mathrm{CIC}}\right)+P\left(N_{\mathrm{CIA}}\right)$	0.04	0.055	0.017	0.03	0.24	0.14
$P\left(N_{\mathrm{CIC}}\right)+P\left(N_2/N_5\right)$	0.04	0.054	0.019	0.03	0.23	0.13
$P\left(N_{\mathrm{CIA}}\right)+P\left(N_2/N_5\right)$	0.03	0.056	0.023	0.03	0.23	0.17
$w_p\left(r_p\right)+\Delta \Sigma \left(r_p\right)+\mathrm{VPF}\left(r\right)+P\left(N_{\mathrm{CIC}}\right)+P\left(N_{\mathrm{CIA}}\right)+P\left(N_2/N_5\right)$	0.02	0.031	0.009	0.01	0.10	0.07

Table 4. Table 4: M r < − 19.5 subscript 𝑀 𝑟 19.5 M_{r}<-19.5

Observable combination	$\sigma \left(\log M_{\min }\right)$	$\sigma \left({\sigma }_{\log M}\right)$	$\sigma \left(\alpha \right)$	$\sigma \left(\log M_1\right)$	$\sigma \left(A_{\mathrm{cen}}\right)$	$\sigma \left(A_{\mathrm{sat}}\right)$
$w_p\left(r_p\right)$	0.24	0.558	0.029	0.03	0.63	0.35
$\Delta \Sigma \left(r_p\right)$	0.07	0.138	0.080	0.08	0.76	0.33
$\mathrm{VPF}\left(r\right)$	0.13	0.179	0.657	0.30	0.62	3.10
$P\left(N_{\mathrm{CIC}}\right)$	0.09	0.119	0.088	0.14	1.06	0.16
$P\left(N_{\mathrm{CIA}}\right)$	0.05	0.141	0.075	0.11	0.72	0.37
$P\left(N_2/N_5\right)$	0.04	0.076	0.052	0.05	0.58	0.28
$w_p\left(r_p\right)+\Delta \Sigma \left(r_p\right)$	0.05	0.104	0.022	0.02	0.40	0.23
$w_p\left(r_p\right)+\mathrm{VPF}\left(r\right)$	0.04	0.077	0.021	0.02	0.32	0.27
$w_p\left(r_p\right)+P\left(N_{\mathrm{CIC}}\right)$	0.02	0.046	0.016	0.02	0.17	0.12
$w_p\left(r_p\right)+P\left(N_{\mathrm{CIA}}\right)$	0.04	0.091	0.018	0.02	0.18	0.22
$w_p\left(r_p\right)+P\left(N_2/N_5\right)$	0.03	0.057	0.017	0.02	0.31	0.17
$\Delta \Sigma \left(r_p\right)+\mathrm{VPF}\left(r\right)$	0.05	0.079	0.047	0.04	0.30	0.21
$\Delta \Sigma \left(r_p\right)+P\left(N_{\mathrm{CIC}}\right)$	0.03	0.050	0.026	0.03	0.20	0.12
$\Delta \Sigma \left(r_p\right)+P\left(N_{\mathrm{CIA}}\right)$	0.04	0.071	0.038	0.03	0.21	0.21
$\Delta \Sigma \left(r_p\right)+P\left(N_2/N_5\right)$	0.03	0.060	0.028	0.03	0.31	0.16
$\mathrm{VPF}\left(r\right)+P\left(N_{\mathrm{CIC}}\right)$	0.04	0.052	0.039	0.06	0.44	0.14
$\mathrm{VPF}\left(r\right)+P\left(N_{\mathrm{CIA}}\right)$	0.04	0.058	0.062	0.04	0.26	0.35
$\mathrm{VPF}\left(r\right)+P\left(N_2/N_5\right)$	0.03	0.052	0.038	0.04	0.34	0.27
$P\left(N_{\mathrm{CIC}}\right)+P\left(N_{\mathrm{CIA}}\right)$	0.03	0.052	0.023	0.03	0.27	0.12
$P\left(N_{\mathrm{CIC}}\right)+P\left(N_2/N_5\right)$	0.03	0.052	0.025	0.03	0.26	0.13
$P\left(N_{\mathrm{CIA}}\right)+P\left(N_2/N_5\right)$	0.03	0.054	0.032	0.03	0.26	0.22
$w_p\left(r_p\right)+\Delta \Sigma \left(r_p\right)+\mathrm{VPF}\left(r\right)+P\left(N_{\mathrm{CIC}}\right)+P\left(N_{\mathrm{CIA}}\right)+P\left(N_2/N_5\right)$	0.01	0.030	0.012	0.01	0.11	0.07

Table 5. Table 5: M r < − 20.0 subscript 𝑀 𝑟 20.0 M_{r}<-20.0

Observable combination	$\sigma \left(\log M_{\min }\right)$	$\sigma \left({\sigma }_{\log M}\right)$	$\sigma \left(\alpha \right)$	$\sigma \left(\log M_1\right)$	$\sigma \left(A_{\mathrm{cen}}\right)$	$\sigma \left(A_{\mathrm{sat}}\right)$
$w_p\left(r_p\right)$	0.17	0.475	0.031	0.03	0.78	0.45
$\Delta \Sigma \left(r_p\right)$	0.06	0.149	0.081	0.07	0.86	0.43
$\mathrm{VPF}\left(r\right)$	0.10	0.155	0.783	0.21	0.56	2.73
$P\left(N_{\mathrm{CIC}}\right)$	0.09	0.157	0.081	0.11	1.32	0.16
$P\left(N_{\mathrm{CIA}}\right)$	0.04	0.137	0.088	0.11	0.91	0.33
$P\left(N_2/N_5\right)$	0.05	0.112	0.064	0.05	0.80	0.32
$w_p\left(r_p\right)+\Delta \Sigma \left(r_p\right)$	0.05	0.123	0.024	0.02	0.53	0.29
$w_p\left(r_p\right)+\mathrm{VPF}\left(r\right)$	0.03	0.082	0.024	0.02	0.39	0.32
$w_p\left(r_p\right)+P\left(N_{\mathrm{CIC}}\right)$	0.03	0.062	0.020	0.01	0.25	0.12
$w_p\left(r_p\right)+P\left(N_{\mathrm{CIA}}\right)$	0.04	0.100	0.021	0.02	0.26	0.22
$w_p\left(r_p\right)+P\left(N_2/N_5\right)$	0.03	0.069	0.021	0.02	0.41	0.19
$\Delta \Sigma \left(r_p\right)+\mathrm{VPF}\left(r\right)$	0.04	0.083	0.054	0.04	0.39	0.25
$\Delta \Sigma \left(r_p\right)+P\left(N_{\mathrm{CIC}}\right)$	0.03	0.065	0.031	0.02	0.30	0.12
$\Delta \Sigma \left(r_p\right)+P\left(N_{\mathrm{CIA}}\right)$	0.03	0.076	0.042	0.03	0.30	0.22
$\Delta \Sigma \left(r_p\right)+P\left(N_2/N_5\right)$	0.03	0.070	0.033	0.02	0.43	0.19
$\mathrm{VPF}\left(r\right)+P\left(N_{\mathrm{CIC}}\right)$	0.03	0.059	0.033	0.04	0.44	0.13
$\mathrm{VPF}\left(r\right)+P\left(N_{\mathrm{CIA}}\right)$	0.03	0.069	0.071	0.04	0.32	0.32
$\mathrm{VPF}\left(r\right)+P\left(N_2/N_5\right)$	0.03	0.069	0.040	0.03	0.38	0.30
$P\left(N_{\mathrm{CIC}}\right)+P\left(N_{\mathrm{CIA}}\right)$	0.03	0.067	0.027	0.03	0.38	0.11
$P\left(N_{\mathrm{CIC}}\right)+P\left(N_2/N_5\right)$	0.03	0.074	0.030	0.03	0.36	0.12
$P\left(N_{\mathrm{CIA}}\right)+P\left(N_2/N_5\right)$	0.03	0.073	0.036	0.03	0.34	0.24
$w_p\left(r_p\right)+\Delta \Sigma \left(r_p\right)+\mathrm{VPF}\left(r\right)+P\left(N_{\mathrm{CIC}}\right)+P\left(N_{\mathrm{CIA}}\right)+P\left(N_2/N_5\right)$	0.01	0.037	0.014	0.01	0.15	0.08

Table 6. Table 6: M r < − 20.5 subscript 𝑀 𝑟 20.5 M_{r}<-20.5

Observable combination	$\sigma \left(\log M_{\min }\right)$	$\sigma \left({\sigma }_{\log M}\right)$	$\sigma \left(\alpha \right)$	$\sigma \left(\log M_1\right)$	$\sigma \left(A_{\mathrm{cen}}\right)$	$\sigma \left(A_{\mathrm{sat}}\right)$
$w_p\left(r_p\right)$	0.05	0.244	0.038	0.03	1.07	0.60
$\Delta \Sigma \left(r_p\right)$	0.05	0.197	0.098	0.06	1.05	0.52
$\mathrm{VPF}\left(r\right)$	0.07	0.117	1.023	0.11	0.71	2.57
$P\left(N_{\mathrm{CIC}}\right)$	0.07	0.282	0.070	0.06	1.53	0.25
$P\left(N_{\mathrm{CIA}}\right)$	0.04	0.209	0.126	0.12	1.88	0.40
$P\left(N_2/N_5\right)$	0.04	0.186	0.060	0.03	1.42	0.41
$w_p\left(r_p\right)+\Delta \Sigma \left(r_p\right)$	0.03	0.132	0.036	0.02	0.70	0.34
$w_p\left(r_p\right)+\mathrm{VPF}\left(r\right)$	0.02	0.084	0.034	0.02	0.53	0.45
$w_p\left(r_p\right)+P\left(N_{\mathrm{CIC}}\right)$	0.03	0.148	0.031	0.02	0.60	0.19
$w_p\left(r_p\right)+P\left(N_{\mathrm{CIA}}\right)$	0.03	0.136	0.030	0.02	0.60	0.30
$w_p\left(r_p\right)+P\left(N_2/N_5\right)$	0.03	0.134	0.029	0.02	0.79	0.25
$\Delta \Sigma \left(r_p\right)+\mathrm{VPF}\left(r\right)$	0.03	0.083	0.065	0.03	0.51	0.32
$\Delta \Sigma \left(r_p\right)+P\left(N_{\mathrm{CIC}}\right)$	0.04	0.140	0.039	0.02	0.57	0.20
$\Delta \Sigma \left(r_p\right)+P\left(N_{\mathrm{CIA}}\right)$	0.02	0.106	0.054	0.03	0.53	0.27
$\Delta \Sigma \left(r_p\right)+P\left(N_2/N_5\right)$	0.03	0.118	0.039	0.02	0.67	0.24
$\mathrm{VPF}\left(r\right)+P\left(N_{\mathrm{CIC}}\right)$	0.02	0.084	0.036	0.03	0.60	0.22
$\mathrm{VPF}\left(r\right)+P\left(N_{\mathrm{CIA}}\right)$	0.02	0.083	0.073	0.03	0.53	0.37
$\mathrm{VPF}\left(r\right)+P\left(N_2/N_5\right)$	0.02	0.088	0.042	0.03	0.58	0.38
$P\left(N_{\mathrm{CIC}}\right)+P\left(N_{\mathrm{CIA}}\right)$	0.03	0.130	0.034	0.03	0.78	0.18
$P\left(N_{\mathrm{CIC}}\right)+P\left(N_2/N_5\right)$	0.03	0.142	0.037	0.02	0.87	0.20
$P\left(N_{\mathrm{CIA}}\right)+P\left(N_2/N_5\right)$	0.03	0.125	0.042	0.03	0.77	0.27
$w_p\left(r_p\right)+\Delta \Sigma \left(r_p\right)+\mathrm{VPF}\left(r\right)+P\left(N_{\mathrm{CIC}}\right)+P\left(N_{\mathrm{CIA}}\right)+P\left(N_2/N_5\right)$	0.01	0.052	0.020	0.01	0.29	0.11

Table 7. Table 7: M r < − 21.0 subscript 𝑀 𝑟 21.0 M_{r}<-21.0

Observable combination	$\sigma \left(\log M_{\min }\right)$	$\sigma \left({\sigma }_{\log M}\right)$	$\sigma \left(\alpha \right)$	$\sigma \left(\log M_1\right)$	$\sigma \left(A_{\mathrm{cen}}\right)$	$\sigma \left(A_{\mathrm{sat}}\right)$
$w_p\left(r_p\right)$	0.19	0.308	0.097	0.03	0.84	0.90
$\Delta \Sigma \left(r_p\right)$	0.08	0.101	0.228	0.07	0.67	0.81
$\mathrm{VPF}\left(r\right)$	0.12	0.149	1.558	0.10	0.53	4.10
$P\left(N_{\mathrm{CIC}}\right)$	0.19	0.273	0.148	0.11	2.46	0.36
$P\left(N_{\mathrm{CIA}}\right)$	0.06	0.070	0.277	0.15	1.48	0.88
$P\left(N_2/N_5\right)$	0.05	0.081	0.091	0.03	0.80	0.62
$w_p\left(r_p\right)+\Delta \Sigma \left(r_p\right)$	0.05	0.074	0.081	0.02	0.38	0.57
$w_p\left(r_p\right)+\mathrm{VPF}\left(r\right)$	0.05	0.073	0.084	0.03	0.41	0.70
$w_p\left(r_p\right)+P\left(N_{\mathrm{CIC}}\right)$	0.06	0.083	0.065	0.02	0.56	0.26
$w_p\left(r_p\right)+P\left(N_{\mathrm{CIA}}\right)$	0.04	0.061	0.076	0.02	0.40	0.50
$w_p\left(r_p\right)+P\left(N_2/N_5\right)$	0.04	0.068	0.063	0.02	0.49	0.37
$\Delta \Sigma \left(r_p\right)+\mathrm{VPF}\left(r\right)$	0.05	0.062	0.159	0.04	0.36	0.56
$\Delta \Sigma \left(r_p\right)+P\left(N_{\mathrm{CIC}}\right)$	0.05	0.064	0.073	0.02	0.41	0.30
$\Delta \Sigma \left(r_p\right)+P\left(N_{\mathrm{CIA}}\right)$	0.03	0.044	0.119	0.03	0.30	0.47
$\Delta \Sigma \left(r_p\right)+P\left(N_2/N_5\right)$	0.04	0.055	0.077	0.02	0.37	0.41
$\mathrm{VPF}\left(r\right)+P\left(N_{\mathrm{CIC}}\right)$	0.06	0.082	0.073	0.03	0.51	0.30
$\mathrm{VPF}\left(r\right)+P\left(N_{\mathrm{CIA}}\right)$	0.04	0.059	0.205	0.04	0.38	0.78
$\mathrm{VPF}\left(r\right)+P\left(N_2/N_5\right)$	0.05	0.067	0.083	0.03	0.40	0.57
$P\left(N_{\mathrm{CIC}}\right)+P\left(N_{\mathrm{CIA}}\right)$	0.04	0.064	0.062	0.03	0.61	0.25
$P\left(N_{\mathrm{CIC}}\right)+P\left(N_2/N_5\right)$	0.05	0.072	0.065	0.03	0.68	0.27
$P\left(N_{\mathrm{CIA}}\right)+P\left(N_2/N_5\right)$	0.04	0.060	0.077	0.03	0.52	0.42
$w_p\left(r_p\right)+\Delta \Sigma \left(r_p\right)+\mathrm{VPF}\left(r\right)+P\left(N_{\mathrm{CIC}}\right)+P\left(N_{\mathrm{CIA}}\right)+P\left(N_2/N_5\right)$	0.02	0.030	0.042	0.01	0.21	0.18

Equations

$$⟨N_{\mathrm{cen}}|M_{\mathrm{vir}}⟩$$

$$⟨N_{\mathrm{sat}}|M_{\mathrm{vir}}⟩$$

$$⟨N_{\mathrm{gal}}|M_{\mathrm{vir}},x>x_{\mathrm{piv}}⟩$$

$$⟨N_{\mathrm{gal}}|M_{\mathrm{vir}},x\le x_{\mathrm{piv}}⟩$$

$$⟨N_{\mathrm{gal}}|M_{\mathrm{vir}}⟩=\int ⟨N_{\mathrm{gal}}|M_{\mathrm{vir}},x⟩P(x|M_{\mathrm{vir}})\mathrm{d}x,$$

$$\delta N_{\mathrm{cen}}$$

$$\delta N_{\mathrm{sat}}$$

$$w_{\mathrm{p}}(r_{\mathrm{p}})=2{\int }_0^{{\pi }_{\max }}d\pi \xi (r_{\mathrm{p}},\pi )$$

$$\Delta \Sigma (r_{\mathrm{p}})=\bar{\Sigma }(<r_{\mathrm{p}})-\Sigma (r_{\mathrm{p}})$$

$$\mathrm{VPF}(r)=\frac{N_{\text{void}}(r)}{N_{\text{sphere}}(r)},$$

$$\mathbf{I}=\frac{\partial {\mathbf{f}}^{\text{T}}}{\partial \mathbf{p}}{\mathbf{C}}^{-1}\frac{\partial \mathbf{f}}{\partial \mathbf{p}},$$

$$\mathbf{\sigma }=\sqrt[]{\text{Diag}(\mathbf{\Sigma })}=\sqrt[]{\text{Diag}({\mathbf{I}}^{-1})}.$$

$${\mathbf{C}}_{\mathrm{total}}={\mathbf{C}}_{\mathrm{jackknife}}+{\mathbf{C}}_{\mathrm{realization}}+{\mathbf{C}}_{\mathrm{random}}.$$

Full text

††pubyear: 2019††pagerange: How to Optimally Constrain Galaxy Assembly Bias: Supplement Projected Correlation Functions with Count-in-cells Statistics–7

1 Introduction

In the concordance $\Lambda$CDM model of the Universe (Komatsu et al., 2011; Betoule et al., 2014; Planck Collaboration et al., 2014, 2018a; DES Collaboration et al., 2018), galaxies reside in dark matter haloes (White & Rees, 1978; Blumenthal et al., 1984), which form around peaks in the primordial dark matter density field (Bardeen et al., 1986; Bond et al., 1991; Sheth et al., 2001a; Sheth et al., 2001b; Zentner, 2007). In practice, the abundance, clustering, and structure of dark matter haloes have been precisely documented by high-resolution, gravity-only $N$-body simulations of cosmological structure growth (e.g., Navarro et al., 1997; Mo & White, 1996; Tinker et al., 2008b; Tinker et al., 2010). Halo occupation models use empirical data to link galaxies to haloes in a statistical sense. They are useful because they provide a convenient means to compare observed galaxy clustering statistics with theoretical predictions without a complete theory of galaxy formation and evolution. This, in turn, is useful because one can use such models to test cosmological models using data on non-linear scales, and because empirical models distil the formidable amount of information available in survey data into a relatively simpler galaxy–halo relationship that can be used to inform models of galaxy formation and evolution. Wechsler & Tinker (2018) provide a contemporary review of these models.

The term assembly bias has, unfortunately, taken on several related but distinct meanings in the literature. The clustering of dark matter haloes is a strong function of halo mass (Kaiser, 1984; Mo & White, 1996), but it has become clear over the last decade that haloes cluster as a function of a number of other properties (Gao et al., 2005; Wechsler et al., 2006; Gao & White, 2007; Li et al., 2008; Mao et al., 2018; Xu & Zheng, 2018b). Gao et al. (2005) first studied the age-dependence of halo clustering, which led to the term assembly bias, but because many halo properties are correlated with formation history, the dependence of halo clustering on many other properties (e.g., concentration, spin, and so on) is often loosely referred to as assembly bias or halo assembly bias as well. Mao et al. (2018) advocate referring to these dependences as secondary biases111The primary bias is the strong dependence of halo clustering on halo mass.. This nomenclature is clearer because the secondary biases do not necessarily have a clear origin in the correlations of halo properties, such as concentration, with conventional measures of halo formation history.

Assembly bias of observed galaxies as well as simulated halo populations has received significant attention in the recent literature regarding the analysis and interpretation of galaxy survey data. Assembly bias may challenge survey analyses because it may induce (1) systematic errors in the inferred galaxy–halo relationship inferred from survey data (e.g., Zentner et al., 2014) and/or (2) biases in inferred cosmological parameters(e.g., Croton et al., 2007; Wu et al., 2008; McCarthy et al., 2018). On the other hand, an unambiguous detection of assembly bias in survey data may pave the way to a richer and more complete understanding of the connection between galaxies and their host dark matter haloes.

In the context of the interpretation of galaxy surveys, assembly bias can lead to the following possibility. Galaxies in a certain luminosity sample may form in haloes with a probability that depends upon not only the mass of the halo, but on any number of halo properties. If this happens, then the resultant clustering of any galaxy sample must be interpreted within the context of a model that incorporates not only the mass dependence of halo clustering, but also the secondary biases, as is done in e.g., Chaves-Montero et al. (2016); Zentner et al. (2019); Lehmann et al. (2017); Romano-Díaz et al. (2017); Desjacques et al. (2018). This case is sometimes loosely referred to as galaxy assembly bias. Signals of galaxy assembly bias are also present in hydrodynamical simulations (e.g., Xu & Zheng, 2018a; Artale et al., 2018; Bose et al., 2019). The issue facing survey data analysis is that the vast majority of studies treat survey data using models that assume that galaxies of a particular luminosity (or other galaxy properties that determine the sample selection) form within haloes with a probability that depends upon only the mass of the halo (and no other halo property). Thus, these analyses account only for the mass dependence of halo clustering (Wechsler & Tinker, 2018).

The most widely used empirical models for interpreting survey data include the halo occupation distribution (HOD) (e.g., Berlind & Weinberg, 2002) and the conditional luminosity function (CLF) (e.g. Yang et al., 2003). Both of these models in the original (and standard) forms assume that galaxies of a particular type reside in haloes with a probability that depends only on the masses of the haloes. Both the HOD (Zheng et al., 2007; Zehavi et al., 2011; Guo et al., 2012a, b; Leauthaud et al., 2012; Guo et al., 2014; Zentner et al., 2019; Vakili et al., 2016) and the CLF (van den Bosch et al., 2003; Yang et al., 2004; Cooray, 2006; van den Bosch et al., 2007; Cacciato et al., 2013; Lange et al., 2019) have been used successfully to interpret a variety of observational samples. On the other hand, subhalo abundance matching (AM or SHAM, e.g., Kravtsov et al., 2004; Tasitsiomi et al., 2004; Vale & Ostriker, 2004; Trujillo-Gomez et al., 2011; Hearin et al., 2013b) has the power to naturally incorporate assembly bias, by matching galaxy properties to halo properties (e.g., $V_{\rm{peak}}$) that can assume varying values at a fixed halo mass. The success of SHAM in interpreting galaxy survey data (Behroozi et al., 2010, 2013c; Reddick et al., 2013; Lehmann et al., 2017) is striking considering the simplicity of its assumptions.

In this paper, we will use an expanded form of the HOD, known as the decorated HOD (dHOD, Hearin et al., 2016). The dHOD builds upon the traditional HOD by adding parameters that enable tunable levels of galaxy assembly bias. We limit our treatment to the HOD and dHOD for specificity and simplicity.

We explore the utility of several galaxy survey observables to constrain assembly bias within the context of simple dHOD (Hearin et al., 2016) models. Aside from the overall galaxy number density $n_{\mathrm{gal}}$, which we use as a basic constraining observable throughout our analyses, the observables that we explore are the projected galaxy correlation function $w_{\mathrm{p}}(r_{\mathrm{p}})$, the excess surface density inferred from galaxy–galaxy lensing $\Delta\Sigma(r_{\mathrm{p}})$, the void probability function $\mathrm{VPF}(r)$, galaxy counts-in-cylinders $P(N_{\mathrm{CIC}})$, galaxy counts-in-annuli $P(N_{\mathrm{CIA}})$, and the probability distribution of the ratio between $N_{\mathrm{CIC}}$ of different cylinder sizes $P(N_{2}/N_{5})$. In this first study of the subject, we intentionally avoid utilizing satellite kinematics (SK), redshift space distortions (RSD), or other observables that require a detailed model of galaxy velocities relative to haloes. Treating such statistics requires additional modelling and additional assumptions that can greatly complicate such a study. For the observables we study, we examine and compare their effectiveness at constraining not only assembly bias within the context of a dHOD model, but nearly all HOD parameters.

The remainder of this paper gives the details necessary to support the summary of our findings stated in the previous paragraph. In Section 2, we describe the simulation that we use, the dHOD models within which we work, the observables we consider, and our approach to estimating parameter constraints. In Section 3, we present our results in detail. We discuss our results in the context of the contemporary literature, draw broad conclusions, and propose future steps in Section 4.

2 Methods

In this section, we give the details of our analysis procedures. This includes a discussion of the simulation that we use, the dHOD models that we explore, the observables that we consider, our methods for estimating parameter constraints, and our estimates of the covariance matrices used in our analyses.

2.1 Simulation

In order to mitigate the limitations of analytic estimates of clustering and lensing statistics (see e.g., van den Bosch et al., 2013), the calculations that we perform in this paper are based upon $N$-body simulations of the formation of structure in a concordance cosmological model. These simulations evolve dark matter particles under the influence of gravity from initial over-densities in the early universe to the present day. In particular, the analysis in this work utilizes the Bolshoi Planck simulation222Available at www.cosmosim.org/cms/simulations/bolshoip (Klypin et al., 2016; Behroozi et al., 2013a; Behroozi et al., 2013b; Rodríguez-Puebla et al., 2016). Bolshoi Planck is a dark matter only simulation within a cubic box of length $250\,h^{-1}\mathrm{Mpc}$, which adopts values of cosmological parameters from Planck Collaboration et al. (2014), namely $\Omega_{\Lambda}$ = 0.693, $\Omega_{\mathrm{m}}$ = $1-\Omega_{\Lambda}$ = 0.307, $\Omega_{\mathrm{b}}$ = 0.048, h = 0.7, $n_{\mathrm{s}}$ = 0.96, and $\sigma_{8}$ = 0.82. The simulation contains 20483 particles, implying a particle mass of $m_{\mathrm{p}}=1.55\times 10{{}^{8}}\,h^{-1}\mathrm{M}_{\odot}$.

We use the bolplanck halo catalogue included with the Halotools software package333halotools.readthedocs.io (Hearin et al., 2017), which also provides an implementation of customizable dHOD models.

The bolplanck halo catalogue was produced from the Bolshoi Planck simulation using the ROCKSTAR halo-finder (Behroozi et al., 2013a). To compute lensing observables, we use the particle catalogue included with Halotools which contains $10^{6}$ randomly-selected particles from the Bolshoi Planck volume, and make a downsampled catalogue containing $\sim 10^{5}$ particles with an acceptance rate of 0.1 for runtime considerations. We have tested that the measured $\Delta\Sigma(r_{\mathrm{p}})$ is not sensitive to the downsampling, and the noise introduced in this process is accounted for in our covariance matrix (see Section 2.5). The catalogues are included in Halotools version halotools_v0p4, which adopts the virial definition of haloes, and we work at $z=0$, corresponding to an overdensity parameter $\Delta_{\textrm{vir}}=333$ with respect to the mean matter density of the Universe.

2.2 Halo Occupation Model

We describe the galaxy–halo connection using the HOD and the dHOD. Both of these models specify the probability for a halo of mass $M_{\mathrm{vir}}$ to host $N_{\mathrm{cen}}$ central galaxies and $N_{\mathrm{sat}}$ satellite galaxies above a certain threshold stellar mass, $P(N_{\mathrm{cen}}|M_{\mathrm{vir}})$ and $P(N_{\mathrm{sat}}|M_{\mathrm{vir}})$ respectively. Central and satellite galaxies are considered separately because central galaxies reside in the potential wells of host haloes while satellite galaxies are associated with subhaloes and experience different physics of formation and evolution. It is well known that subhaloes experience very distinct evolution from host haloes and thus have demographics that are distinct from host haloes (Kravtsov et al., 2004; Zheng et al., 2005; Zentner et al., 2005). Moreover, numerous observations, using many different approaches, have established that central galaxies and satellites have different properties and can be described as two distinct populations (Oemler, 1974; Dressler, 1980; Postman & Geller, 1984; Zehavi et al., 2005; Weinmann et al., 2006; van den Bosch et al., 2008).

The central galaxy occupation is modelled as a Bernoulli random variable, which takes the value 1 with probability $p$ and the value 0 with probability $1-p$, with $p=\langle N_{\mathrm{cen}}|M_{\mathrm{vir}}\rangle$. Satellite galaxies follow a Poisson distribution with first moment $\langle N_{\mathrm{sat}}|M_{\mathrm{vir}}\rangle$.

In principle, the formalism we use in this work is the same as in Zentner et al. (2019) which, in turn, was chosen to mimic the analysis of Zehavi et al. (2011). However, we adopt the implementation of the model in Halotools, and introduce some subtle modifications which will be elaborated on in Section 2.2.2.

2.2.1 Standard HOD

In the standard HOD, which does not account for any potential galaxy assembly bias, the mass of a halo solely determines the galaxy occupation. The mean central and satellite galaxy occupations vary with halo mass according to

[TABLE]

(Zheng et al., 2007), where $M_{\mathrm{min}}$ is the mass at which a halo has a 50% probability of hosting a central galaxy; $\sigma_{\log M}$ is a measure for the scatter in the stellar mass–halo mass relation that determines the steepness of the $\langle N_{\mathrm{cen}}|M_{\mathrm{vir}}\rangle$ transition from zero to unity; $M_{0}$ is the truncating mass, below which $\langle N_{\mathrm{sat}}|M_{\mathrm{vir}}\rangle=0$; the mass $M_{1}$ indicates the halo mass at which there is, on average, one satellite444More accurately, this mass is $M_{1}$+$M_{0}$, but $M_{0}$ is typically much smaller than $M_{1}$. if a central is present; and, finally, $\alpha$ is the index of the satellite occupation power law. Note that Eq. (2) expresses the probability of having a satellite galaxy for a halo with mass $M_{\mathrm{vir}}$, after marginalizing over the central occupation. The first term on the right hand side indicates the mean satellite occupation in haloes with a central galaxy, while the second term modulates this occupation by the probability for a halo to contain such a central. Hence, the presence of a central boosts the probability for a halo to host satellite galaxies. Note, though, that for individual haloes a central galaxy is not strictly required for satellites to be present. Although this modulation with $\langle N_{\mathrm{cen}}|M_{\mathrm{vir}}\rangle$ is fairly common (e.g., Zheng et al., 2007; Zehavi et al., 2011; Zentner et al., 2019), we emphasize that it is not used by all authors.

These specifications, along with the assumptions that the central galaxy HOD is a Bernoulli distribution and the satellite galaxy HOD is a Poisson distribution, suffice to specify fully the halo occupation statistics of dark matter haloes in a standard HOD model without assembly bias.

2.2.2 Decorated HOD

Galaxy assembly bias can be incorporated into the HOD formalism in any number of ways. For a secondary halo property $x$, (e.g., concentration, spin, etc.), one can specify a functional form for the probability distributions $P(N_{\mathrm{cen}}|M_{\mathrm{vir}},x)$ and $P(N_{\mathrm{sat}}|M_{\mathrm{vir}},x)$. In such a generalized HOD, the clustering of galaxies can be altered if halo clustering depends upon secondary property $x$. The decorated HOD (dHOD, Hearin et al., 2016) is one way of incorporating assembly bias into the HOD formalism such that integrating the dHOD probability distributions over the secondary properties of interest yields the standard HOD.

In the present paper, we use a simple variation of the dHOD as an illustrative model. In particular, we divide haloes into two categories based upon secondary halo property $x$. Haloes with higher values of $x$ are assigned distinct HODs compared to haloes with lower values of $x$, with a pivot value of $x_{\mathrm{piv}}$. This is the “discrete halo subpopulations” example discussed in Section 4.2 of Hearin et al. (2016) and used to analyze SDSS data (Abazajian et al., 2009) in Zentner et al. (2019). To specify completely the dHOD, we assume that $P(N_{\mathrm{cen}}|M_{\mathrm{vir}},x)$ is a Bernoulli distribution and that $P(N_{\mathrm{sat}}|M_{\mathrm{vir}},x)$ is a Poisson distribution, but that these distributions have first moments of

[TABLE]

where we use the notation $N_{\textrm{gal}}$ because this modification applies equally well to both the central and satellite occupations. We choose $x_{\mathrm{piv}}$ to be the median value of $x$ at a given halo mass, so that each population contains 50% of all the haloes. In this toy model, assembly bias manifests itself as a step function in the secondary property $x$, though we expect that any assembly bias realized in nature would be represented by a smooth function of $x$. This simple model is practical in the sense that current data are not sufficient to constrain more complex models (Zentner et al., 2019); however, richer models of assembly bias can naturally be accommodated within the dHOD framework (Hearin et al., 2016) and future data sets are likely to enable constraints on richer models.

The differences $\delta N_{\mathrm{gal}}$ above are characterized by two assembly bias parameters, $A_{\mathrm{cen}}$ and $A_{\mathrm{sat}}$, both constructed so that they range between 1 and -1, in addition to the five standard HOD parameters. A list of the 7 dHOD parameters can be found in Table 1. Positive values of $A_{\mathrm{gal}}$ indicate a positive correlation between galaxy number and halo property $x$ (i.e., haloes with $x>x_{\mathrm{piv}}$ contain more galaxies, on average, than those with $x<x_{\mathrm{piv}}$), while negative values represent anti-correlation. When $A_{\mathrm{gal}}=0$, the model reduces to the traditional standard HOD. Note that $A_{\mathrm{cen}}$ and $A_{\mathrm{sat}}$ vary independently of one another and do not necessarily have the same sign. The stipulations that the occupation of a halo never be negative, and the requirement that

[TABLE]

with $P(x|M_{\mathrm{vir}})$ the probability distribution for $x$ given $M_{\mathrm{vir}}$, implies that

[TABLE]

It should be noted that, when populating a mock galaxy catalogue using HOD or dHOD, the actual number of galaxies in each halo is a random variable: the number of central galaxies follows the Bernoulli distribution, and the number of satellite galaxies follows the Poisson distribution. Since we will be conducting a Fisher analysis, the random fluctuation in realizations can masquerade as a dependence of galaxy number density on (d)HOD parameters, and yield artificially tight constraints. Hence, we need to reduce the random fluctuation in realizations as much as possible so that a small change in one or more (d)HOD parameters results in a small change in the total number of galaxies. We achieve this by assigning to each halo two random variates, $p_{\text{cen}}$ and $p_{\text{sat}}$, both drawn from the uniform distribution $U(0,1)$, independently from the (d)HOD parameter values. We then find the number that corresponds to these $p$-values in the cumulative distribution of a Bernoulli distribution (for $p_{\text{cen}}$) or a Poisson distribution (for $p_{\text{sat}}$). This minimizes the random fluctuations among realizations that only differ slightly in their corresponding (d)HOD parameters.

In the case of the dHOD, the mean number density of galaxies is strictly independent of the dHOD parameters $A_{\mathrm{cen}}$ and $A_{\mathrm{sat}}$, and so the problem of preserving the total number density from one mock realization to another is particularly acute. In the dHOD, changes to $A_{\mathrm{cen}}$ or $A_{\mathrm{sat}}$ result in changes to the mean occupations of individual haloes, but should result in no change to the total number density. If the galaxy occupation for each halo is realized independently, then the total number of galaxies can vary from mock realization to another as $A_{\mathrm{cen}}$ and/or $A_{\mathrm{sat}}$ are varied. The result of such a variation would be to infer additional constraining power on $A_{\mathrm{cen}}$ and $A_{\mathrm{sat}}$ where there should be none. To mitigate this possibility, we slightly modified the dHOD implementation in Halotools555Our implementation is called PreservingNgalHeavisideAssembias in Halotools. to ensure the total number density of galaxies is preserved among mock catalogs that differ only in their values of $A_{\mathrm{cen}}$ and $A_{\mathrm{sat}}$. We achieved this by conditioning the dHOD on the total number of galaxies before realizing the occupation of each individual halo. It should be noted that once the total number of central galaxies is fixed, the number of central galaxies in each halo would no longer be strictly a Bernoulli distribution. However, for satellite galaxies, both the total number of galaxies and the number of galaxies in an individual halo follow Poisson distributions.

In this work, we choose the NFW concentration parameter (Navarro et al., 1997) as our secondary property [so $x=c_{\rm{NFW}}$ in Eqs. (3) and (4)] when studying constraints on the parameters $A_{\mathrm{cen}}$ and $A_{\mathrm{sat}}$. As has been shown in Croton et al. (2007) and Villarreal et al. (2017), concentration only partially accounts for galaxy assembly bias, and other halo properties (e.g., halo age, spin, environment density) may also contribute to assembly bias. Nevertheless, we choose concentration for several physically motivated reasons. First, concentration is known to correlate with assembly history (Wechsler et al., 2002)666Though this does not guarantee that concentration and assembly history metrics will lead to similar secondary biases (Mao et al., 2018)., and has the advantage that it can be measured in a single snapshot of a simulation. Second, the success of abundance matching suggests that the HODs realized by nature may, indeed, have some dependence upon halo structure (Conroy et al., 2006; Reddick et al., 2013; Hearin & Watson, 2013; Zentner et al., 2014; Mao et al., 2015). Indeed, Lehmann et al. (2017) showed that abundance matching in a manner that does not include any concentration dependence is excluded by galaxy clustering. Third, haloes are known to exhibit large concentration-dependent clustering in the mass range of interest to us ($\sim$ a few $\times 10^{12}\ h^{-1}\mathrm{M}_{\odot}$). Consequently, concentration-dependent clustering is an excellent test case with which to study methods to constrain assembly bias. Fourth, concentration-dependent clustering has already been studied in Hearin et al. (2016) for the dHOD and for observational samples by Zentner et al. (2019) and Vakili et al. (2016), providing a baseline for comparison. Given the above reasons, we believe the concentration parameter is the most reasonable choice for this study, yet we note that our findings may not be trivially generalized to the assembly bias induced by other secondary halo properties, as other properties may induce different assembly bias behaviours.

2.2.3 Spatial and Velocity Distribution

The detailed predictions of an empirical model depend not only on the model for halo occupation, but also upon the positions and velocities, relative to the host halo, that are assigned to the galaxies. We place the central galaxy at the halo centre and the central galaxy inherits the host halo’s peculiar velocity. Satellite galaxies are distributed within the virial radius of the host halo according to a spherically symmetric NFW profile characterized by the same concentration as the dark matter distribution. This assumption is supported by various works (e.g., van der Marel et al., 2000; Lin et al., 2004; van den Bosch et al., 2005), though other authors find that the distribution of satellite galaxies are described by a concentration different from that of dark matter particles, depending on the satellite population (Chen, 2008; More et al., 2009; Watson et al., 2012; Tal et al., 2012; Lange et al., 2018). The radial velocity distribution of satellite galaxies is modelled as a Gaussian distribution with the host halo velocity as the first moment and the solution of the isotropic Jeans equation for an NFW profile (Klypin et al., 1999) as the second moment. We assume velocities to be isotropic, and draw the peculiar velocities in each Cartesian direction independently from this distribution. In practice, the statistics that we examine are quite insensitive to moderate alterations to the treatment of galaxy peculiar velocities (this is by design), though it would be interesting to explore statistics that are sensitive to peculiar velocities as a follow-up study. To examine the effect of alternative velocity models, we have tested the velocity bias model in Guo et al. (2015) with $\alpha_{c}=0.3$ and $\alpha_{s}=1$, and find that the systematic change in our observables is negligible (within $1.5\%$ in all cases).

2.3 Observables

In search of effective ways of utilising existing and future galaxy surveys to constrain the dHOD, we consider a number of observables that are sensitive to halo occupation. In particular, while including the overall galaxy number density of the simulation volume, $n_{\mathrm{gal}}$, as a constraining observable in all of our analyses, we examine

the projected two-point correlation function, $w_{\mathrm{p}}(r_{\mathrm{p}})$; 2. 2.

the galaxy–galaxy lensing signal, $\Delta\Sigma(r_{\mathrm{p}})$; 3. 3.

the void probability function $\mathrm{VPF}(r)$; 4. 4.

the distribution of counts-in-cylinders, $P(N_{\mathrm{CIC}})$; 5. 5.

the distribution of counts-in-annuli, $P(N_{\mathrm{CIA}})$ (analogous to counts-in-cylinders, but with an excised inner region); 6. 6.

the distribution of the ratio of counts in cylinders of different sizes $P(N_{2}/N_{5})$.

We discuss these observables in more detail in the remainder of this subsection.

We compute all observables numerically, by generating mock galaxy catalogues and subsequently measuring each observable from the mock catalogue. This forward-modelling approach enables us to mitigate modelling uncertainty associated with analytic approaches to galaxy clustering and to incorporate possible systematic errors into our calculations. All observables are computed in redshift space, as they would be from observational data, namely, the coordinates of galaxies $(x,y,z)$ are mapped onto $(x,y,z+v_{z}/aH(a))$. We show examples of the measured values of the observables and their uncertainties from jackknife subsampling (see Section 2.5) in Fig. 1 for our fiducial HOD models. Our fiducial models are taken from the fits of Zentner et al. (2019), the parameters of which are listed in Table 1.

Each of the observables is binned in a particular manner. We have selected the binning scheme to ensure that our binning does not significantly degrade the constraining power of any individual observable. We do this by performing a series of analyses in which the bin sizes are reduced in each analysis. We choose bin sizes for each observable such that further refinement of the bins would not yield significant improvement in parameter constraints. We specify the range of the independent variable for each observable observable (for example, in the case of $w_{\mathrm{p}}(r_{\mathrm{p}})$, we take $0.1\leq r_{\mathrm{p}}/h^{-1}\mathrm{Mpc}\leq 31.6$) and increase the number of bins until parameter constraints saturate. This process has been described in detail in Hearin et al. (2013a). We find that the constraining power of all observables saturates at fewer than 30 bins, so we take 30 bins for all observables for simplicity. The binning scheme for which our main results are obtained is shown in Table 2. 777Our results are insensitive to the largest length scales included in our analysis because statistics on these scales are measured with relatively low signal-to-noise. We have verified that excluding the few largest bins of $w_{\mathrm{p}}(r_{\mathrm{p}})$ and $\Delta\Sigma(r_{\mathrm{p}})$ from our analyses results in negligible quantitative change in our constraints (typically below $1\%$, and as large as $\sim 2\%$ in the most extreme cases), and no qualitative change to our conclusions.

2.3.1 Projected Two-Point Correlation Function, $w_{\mathrm{p}}(r_{\mathrm{p}})$

The projected two-point correlation function, $w_{\mathrm{p}}(r_{\mathrm{p}})$, is a canonical observable that has been considered in numerous previous analyses to inform halo occupation (e.g., van den Bosch et al., 2007; Zheng et al., 2009; Zehavi et al., 2011). It is defined by

[TABLE]

where $\xi(r_{\mathrm{p}},\pi)$ is the excess probability of finding galaxy pairs with projected and line-of-sight separations $r_{\mathrm{p}}$ and $\pi$, respectively. We estimate $w_{\mathrm{p}}(r_{\mathrm{p}})$ from our mock catalogues by counting galaxy pairs that have a projected separation in a bin of $r_{\mathrm{p}}$ within a perpendicular distance of $\pi_{\mathrm{max}}$ in redshift space. We choose $\pi_{\mathrm{max}}=60h^{-1}\mathrm{Mpc}$, as is done by Zehavi et al. (2011), according to whom this integration limit is large enough to include most correlated pairs and minimize the impact of the details of peculiar velocity models, yet sufficiently small to give a stable result by suppressing noise from very distant, uncorrelated pairs. We compute $w_{\mathrm{p}}(r_{\mathrm{p}})$ in 29 logarithmically spaced radial bins from $r_{\mathrm{p}}=0.1\,h^{-1}\mathrm{Mpc}$ to $r_{\mathrm{p}}=31.6\,h^{-1}\mathrm{Mpc}$.

The projected two-point clustering of our fiducial models are shown for two luminosity threshold samples in the upper, left panel of Fig. 1. The figure exhibits several well-known characteristics of galaxy clustering. First, brighter galaxies cluster more strongly. Second, the galaxy two-point correlation function can roughly be described as a power law, $w_{\mathrm{p}}=(r_{\mathrm{p}}/r_{0})^{\alpha}$, with index $\alpha\approx-0.8$. Third, in more detail, the correlation function exhibits a small deviation from a power law near $r_{\mathrm{p}}\sim 2\,h^{-1}\mathrm{Mpc}$ which is due to the transition from galaxy pairs that reside in distinct haloes (the “two-halo” term) on large scales ($r_{\mathrm{p}}\gtrsim 2h^{-1}\mathrm{Mpc}$) and pairs of galaxies that reside in a common halo (the “one-halo term”) on scales $r_{\mathrm{p}}\lesssim 1h^{-1}\mathrm{Mpc}$.

2.3.2 Galaxy–Galaxy Weak Lensing, $\Delta\Sigma(r_{\mathrm{p}})$

In addition to the projected two-point clustering, galaxy–galaxy weak lensing is another observable statistic that has been used by many previous authors to constrain halo occupation from observational data (e.g., Mandelbaum et al., 2006a; Cacciato et al., 2009; van Uitert et al., 2016). The canonical observable, $\Delta\Sigma(r_{\mathrm{p}})$, is the excess surface density of mass around galaxies projected along the line-of-sight, and averaged over all potential lens galaxies in the sample,

[TABLE]

where $\Sigma(\textit{r}_{\textrm{p}})$ is the projected surface density evaluated at position $r_{\mathrm{p}}$ relative to the centre of the lens galaxy, and $\bar{\Sigma}(<r_{\mathrm{p}})$ is the mean projected, two-dimensional, surface mass density within a projected distance or $r_{\mathrm{p}}$ from the lens galaxy. We compute $\Delta\Sigma(r_{\mathrm{p}})$ in 30 logarithmically-spaced radial bins from $r_{\mathrm{p}}=0.1\,h^{-1}\mathrm{Mpc}$ to $r_{\mathrm{p}}=31.6\,h^{-1}\mathrm{Mpc}$. The simulations that we use are gravity-only N-body simulations, so our estimates of $\Delta\Sigma(r_{\mathrm{p}})$ include neither baryonic mass nor any influences of baryons on the dark matter distribution (e.g., Rudd et al., 2008).

The galaxy–galaxy lensing signal in our fiducial catalogues is depicted in the upper, right panel of Fig. 1. It is evident from this panel that galaxies in the higher-luminosity samples are more strongly correlated with mass, indicating the well-known fact that more luminous galaxies tend to reside in more massive dark matter haloes. The feature due to the transition between the one-halo and two-halo terms in the galaxy matter correlation function is evident near $r_{\mathrm{p}}\sim 2\,h^{-1}\mathrm{Mpc}$ as well.

2.3.3 Void Probability Function $\mathrm{VPF}(r)$

We examine a number of options for observables in addition to the canonical $w_{\mathrm{p}}(r_{\mathrm{p}})$ and $\Delta\Sigma(r_{\mathrm{p}})$, among them the VPF. The VPF has been examined in previous studies on assembly bias with mixed conclusions (Tinker et al., 2013; Zentner et al., 2014) and, in principle, depends upon all of the n-point functions (e.g., Peebles, 1980). To estimate $\mathrm{VPF}(r)$, we randomly place spheres of radius $r$ throughout our simulation volume and enumerate the probability of the spheres containing zero galaxies (and thus being classified as voids) as

[TABLE]

where $N_{\mathrm{sphere}}$ is the total number of spheres that we use for the estimate ($N_{\mathrm{sphere}}=10^{5}$ in this work) and $N_{\textrm{void}}$ is the number of spheres that are found to enclose zero galaxies. We compute $\mathrm{VPF}(r)$ at 30 logarithmically-spaced radii from $r=1\,h^{-1}\mathrm{Mpc}$ to $r=10\,h^{-1}\mathrm{Mpc}$. We remind the reader that these calculations are performed in redshift space, by mapping the coordinates of galaxies according to their line-of-sight velocities.

The $\mathrm{VPF}(r)$ of our fiducial models are depicted in the left, middle panel of Fig. 1. As with $w_{\mathrm{p}}(r_{\mathrm{p}})$ and $\Delta\Sigma(r_{\mathrm{p}})$, several expected features of the $\mathrm{VPF}(r)$ are evident. The VPF drops from nearly unity on small scales to well below unity beyond a scale of $r\sim 10\,h^{-1}\mathrm{Mpc}$ and voids are more likely for higher luminosity galaxy samples, due largely to their overall lower number density.

2.3.4 Counts-in-cylinders (CIC) Statistic $P(N_{\mathrm{CIC}})$

Galaxy counts, particularly counts of galaxies within cylindrical volumes in redshift space, have been studied for decades (e.g., Fry & Peebles, 1978; Balian & Schaeffer, 1988; Alimi et al., 1990; Baugh et al., 1995; Colombi et al., 1995; Szapudi et al., 1996; Kim & Strauss, 1998; Hogg et al., 2004; Kauffmann et al., 2004; Blanton et al., 2006; Barton et al., 2007; Reid & Spergel, 2009; Berrier et al., 2011; Oguri & Lin, 2015; Gruen et al., 2018). The average number of companions that a galaxy will have within a particular cylinder can be computed from the two-point correlation function; however, the distribution of counts-in-cylinders depends, at least in principle, upon all of the higher n-point functions (Peebles, 1980) and can complement the two-point function as a study of the galaxy halo relationship.

We compute counts in cylinders (CIC) from our galaxy catalogues as follows. We centre a cylinder of transverse radius $r_{\mathrm{CIC}}$ and depth $\pm\Delta v$ (in redshift space) on each galaxy in the sample and count the number of companion galaxies that fall within the cylinder. This procedure enables us to estimate a probability distribution of companion number, $P(N_{\mathrm{CIC}})$, which is the probability that any galaxy has $N_{\mathrm{CIC}}$ companions within the cylinder. For the primary results that we present in this paper, we use $r_{\mathrm{CIC}}=2h^{-1}\mathrm{Mpc}$ and a maximum relative velocity of $\Delta v=1000\,\mathrm{km\ s^{-1}}$, corresponding to a half-length of $L=10h^{-1}\mathrm{Mpc}$, assuming velocities are only due to the Hubble flow. We choose cylinders of a transverse radius $r_{\mathrm{CIC}}$ on the order of a few $h^{-1}\mathrm{Mpc}$ in order to include galaxy companions separated by a scale on which assembly bias is known to introduce a distinct feature in halo clustering (Hearin et al., 2016; Zentner et al., 2019; Sunayama et al., 2016). We have experimented with a variety of alternative cylinder radii and depths, finding that our results remain qualitatively similar. When characterizing the PDF $P(N_{\mathrm{CIC}})$, we count the first few values of $N_{\mathrm{CIC}}$ individually and group larger $N_{\mathrm{CIC}}$ values into logarithmically spaced bins, as indicated in Table 2.

Examples of our counts-in-cylinders distributions, $P(N_{\mathrm{CIC}})$, for two of our luminosity threshold samples, can be seen in the middle, right-hand panel of Fig. 1. As is expected, the probability of having a large number of companions in a cylindrical cell increases dramatically with decreasing galaxy luminosity due to the higher number density of galaxies with lower luminosities.

2.3.5 Counts-in-annuli (CIA) Statistic $P(N_{\mathrm{CIA}})$

To complement counts-in-cylinders, we also examine counts of neighbour galaxies in annuli. The counts-in-annuli (CIA) enable one to get a sense of clustering as a function of scale and to compare smaller-scale, intra-halo clustering to larger-scale clustering. The statistic $P(N_{\mathrm{CIA}})$ is the probability that the number of companions within the annulus is equal to $N_{\mathrm{CIA}}$, analogous to $P(N_{\mathrm{CIC}})$. Unlike the counts-in-cylinders statistic, which roughly probes the halo-occupation statistics on the “one-halo” scale, $P(N_{\mathrm{CIA}})$ is a novel statistic, introduced here, and specifically designed to probe the immediate, supra-halo environments of galaxies.

As with CIC, we choose fixed dimensions for the annuli that we use and explore the constraining power of the distribution of counts around galaxies in our catalogues. Our annuli have inner radii of $r_{\mathrm{inner}}=2\,h^{-1}\mathrm{Mpc}$ and outer radii of $r_{\mathrm{outer}}=5\,h^{-1}\mathrm{Mpc}$. As with CIC, the annuli have a depth in the redshift dimension of $10\,h^{-1}\mathrm{Mpc}$, corresponding to a velocity difference of $\Delta v=1000\,\mathrm{km\ s^{-1}}$. As with CIC, this geometry is chosen in order to probe the immediate environments of haloes, particularly on scales where assembly bias has already been shown to induce a feature in galaxy clustering (e.g., Hearin et al., 2016; Sunayama et al., 2016). We have experimented with moderately different annular dimensions and obtained qualitatively similar results in all cases. We group values of $N_{\mathrm{CIA}}$ in a similar way to $N_{\mathrm{CIC}}$, as detailed in Table 2.

Examples of $P(N_{\mathrm{CIA}})$ for our fiducial catalogues are given in the lower, left-hand panel of Fig. 1. The CIA distribution shares most of the qualitative features of the CIC distribution, though the counts are generally higher because the volumes of our annuli exceed the volumes of our cylinders by a factor of $\sim 5$.

2.3.6 Distribution of Cylinder Count Ratios, $P(N_{2}/N_{5})$

As a distinct way of characterising the clustering environments of galaxies, we also consider the distribution of the ratio of two cylinder counts. The first count is within a cylinder with a radius of $r_{\mathrm{CIC}}=2\,h^{-1}\mathrm{Mpc}$ and the second, larger cylinder has a radius $r_{\mathrm{CIC}}=5\,h^{-1}\mathrm{Mpc}$ for each galaxy. For both cylinders we adopt the same depth, $\Delta v=1000\,\mathrm{km\ s^{-1}}$ as for the CIC and CIA statistics discussed above. For each galaxy in our catalogues, we compute the companion counts $N_{2}$ and $N_{5}$ within each of these cylinders, and take the probability distribution of the ratio of these two numbers as the statistic of interest (notice that the inner cylinder is the same cylinder used in our CIC calculations, so that $N_{5}=N_{2}+N_{\mathrm{CIA}}$. Similar to $N_{\mathrm{CIA}}$, this is a novel statistic to probe the large scale distribution of galaxies that, to the best of our knowledge, has not been utilized before. The intention of this statistic is to probe the relative clustering within a halo (the “one-halo term”) to that in its immediate vicinity. We measure the probability distribution of this ratio, $P(N_{2}/N_{5})$ in 30 linearly-spaced bins from 0 to 1.

Probability distributions of the cylinder count ratio are shown in the lower, right-hand panel of Fig. 1. It is evident that the higher luminosity sample has a distribution with more discreteness noise, as a result of the low number density. The more luminous sample has a significantly higher probability of having low values of $N_{2}/N_{5}$ than the lower luminosity sample. This is due to the fact that satellite galaxies are increasingly rare in the higher-luminosity samples, so that $N_{2}$ is increasingly likely to be either small or zero in such samples compared to lower-luminosity samples (see the right, middle panel of Fig. 1).

2.4 Fisher Analysis

We use a Fisher matrix analysis (Bond et al., 1998; Dodelson, 2003) to forecast the constraining power of each of the observables described above and combinations thereof. Despite its approximate nature, we elect to use a Fisher matrix due to the computational expense of utilizing other techniques. Employing a technique that directly quantifies uncertainty by averaging the posteriors found from each mock catalogue, either via Markov Chain Monte Carlo (Goodman & Weare, 2010; Foreman-Mackey et al., 2013) or Approximate Bayesian Computation (Weyant et al., 2013), to perform the very large number of analyses that we undertake using mock catalogues is substantially more computationally intensive.

The Fisher matrix is a measure of the ideal amount of information that can be obtained from an experiment. The Fisher matrix is defined as

[TABLE]

where $\mathbf{f}$ is the set of observables, $\mathbf{C}$ is the covariance of the observables, and $\mathbf{p}$ is the parameter set to be constrained. The set of observables $\mathbf{f}$ includes all bins of each observable quantity (e.g., 29 bins of $w_{\mathrm{p}}(r_{\mathrm{p}})$, 30 bins of $\Delta\Sigma(r_{\mathrm{p}})$, etc.), while $\mathbf{p}$ represents the set of all model parameters. The notation $\partial\mathbf{f}/\partial\mathbf{p}$ represents the matrix of values constructed by differentiating each observable with respect to each of the individual parameters, so that the matrix element $\partial f_{\mathrm{i}}/\partial p_{\mathrm{j}}$ represents the derivative of the $i^{\mathrm{th}}$ observable with respect to the $j^{\mathrm{th}}$ model parameter. Both the derivatives and the covariance are evaluated at a single, fiducial point in the parameter space, which is assumed to be the true underlying model.

The expected $1\sigma$ error on any inferred parameter, marginalized over all other parameters, can be obtained by taking the square root of the corresponding diagonal term of the posterior covariance matrix $\bm{\Sigma}$, which is the inverse of the Fisher matrix. Hence, the forecasted marginalized uncertainties in the parameters are

[TABLE]

In our study of the constraining power of various observables, we explore the 6 dimensional parameter space, spanned by $\alpha$, $\log{M_{\mathrm{1}}}$, $\sigma_{\log{M}}$, $\log{M_{\rm{min}}}$, $A_{\rm{cen}}$, and $A_{\rm{sat}}$. The parameter $\log{M_{\mathrm{0}}}$ is part of both the standard HOD and the dHOD models. However, we set $\log{M_{0}}$ to its fiducial value in our analyses and do not allow it to vary. We do this because $\log{M_{0}}$ is poorly constrained by these data (see Zentner et al., 2019).

The fiducial points about which we evaluate our Fisher matrices are given by the parameters listed in Table 1, and differ for each luminosity threshold sample. In the next two subsections, we discuss the computation of the covariance matrix, $\mathbf{C}$, and the derivatives of the observables, $\partial\mathbf{f}/\partial\mathbf{p}$, respectively.

2.5 Covariance

In order to implement the Fisher approximation for the marginalized constraints on model parameters, we must compute a covariance matrix about the fiducial point in the parameter space. The covariance matrix that we calculate has three contributions. The first, and dominant, contribution is from sample variance (sometimes called “cosmic” variance in this context). We estimate the sample variance contribution using jackknife resampling of the simulation volume, while recognizing the caveat that jackknife resampling is known to underestimate covariances. We will refer to this component of the covariance as $\mathbf{C}_{\mathrm{jackknife}}$. The second contribution to the covariance matrix is due to the stochasticity of populating a simulation with galaxies drawn from the probability distribution functions of the (d)HOD. Multiple realizations of the same underlying model in identical, finite volumes will lead to mildly different predictions due to this stochasticity. We refer to this contribution to the covariance as $\mathbf{C}_{\mathrm{realization}}$. Third, we use a fixed set of randomly distributed centres of spheres in the calculation of $\mathrm{VPF}(r)$ as well as a fixed subsample of dark matter particles in the calculation of $\Delta\Sigma(r_{\mathrm{p}})$. These choices introduce a small contribution to the covariance that we denote $\mathbf{C}_{\mathrm{random}}$. The total covariance matrix that we use is the sum of each of these contributions

[TABLE]

As an example, Fig. 2 depicts the covariance matrices $\mathbf{C}_{\mathrm{jackknife}}$, $\mathbf{C}_{\mathrm{random}}$, $\mathbf{C}_{\mathrm{realization}}$ and their linear combination $\mathbf{C}_{\mathrm{total}}$, for the $M_{r}<-19.0$ threshold sample, as normalized correlation matrices.

The contributions from $\mathbf{C}_{\mathrm{realization}}$ and $\mathbf{C}_{\mathrm{random}}$ are straightforward to compute. To estimate $\mathbf{C}_{\mathrm{realization}}$, we populate the halo catalogue with the fiducial HOD multiple times, each time using a new random seed, and compute the covariance across the measurements from the resultant mock galaxy catalogues. $\mathbf{C}_{\mathrm{realization}}$ is displayed in the upper right panel of Fig. 2. For $\mathbf{C}_{\mathrm{random}}$, we repeat measurements of $\Delta\Sigma(r_{\mathrm{p}})$ and $\mathrm{VPF}(r)$ on the same mock galaxy catalogue, but with different sets of particle subsamples and sphere centres, respectively, and calculate the covariance. Note that this matrix only consists of the two corresponding blocks, as shown in the lower left panel of Fig. 2.

Finally, to estimate $\mathbf{C}_{\mathrm{jackknife}}$, which is shown in the upper left panel of Fig. 2, we divide the simulation box into 10$\times$10 cuboid cells, each of size $25\,h^{-1}\mathrm{Mpc}\times 25\,h^{-1}\mathrm{Mpc}\times 250\,h^{-1}\mathrm{Mpc}$. The long axes of each cuboid are the same as the length of the simulation volume and are assumed to lie along the line of sight. For each mock catalogue, we construct three such sets of jackknife samples by choosing, in turn, the $x$, $y$, and $z$ dimensions of the simulation cube as the line-of-sight direction. Our final covariances are the averages of the three covariances computed for each of the three projections. We construct this average to minimize the contributions from any significant variations that may, by chance, fall along any individual projection. For each set of jackknife samples coming from each of the three projections of the mock catalogues, we exclude individual jackknife cells in turn, and compute the jackknife contribution in the usual manner (Jones, 1956; Quenouille, 1956; Tukey, 1958; Wall & Jenkins, 2003).

For the purposes of computing jackknife covariances only, the mock catalogues that we use are not based on our fiducial HODs. Our jackknife covariance mock catalogues are based upon abundance matching with zero scatter (Kravtsov et al., 2004). We construct these catalogues by populating haloes that have the highest values of $V_{\mathrm{peak}}$, with galaxy number densities consistent with HOD realizations. This modification is necessary for the following reason. Each mock HOD-based catalogue is a realization of the underlying HODs. Therefore, there is inherent stochasticity in the covariance matrix estimates. Moreover, in the HOD formalism, each luminosity threshold must be treated independently, which, in turn, means that the covariances in different threshold samples can fluctuate independently. This makes comparing covariances across luminosity thresholds challenging because to do this using the HOD approach requires marginalizing the stochasticity over a very large number of mock catalogues. The abundance matching approach that we have adopted allows us to circumvent this difficulty because there is no stochasticity in the mock catalogues. Therefore, the stochasticity associated with building mock catalogues does not contribute to our $\mathbf{C}_{\mathrm{jackknife}}$ estimates. This ensures that our $\mathbf{C}_{\mathrm{jackknife}}$ estimates vary smoothly with the luminosity threshold of the sample. We have found that this procedure reduces the noisiness of our forecasts, yet does not alter our qualitative results.

All six of our candidate observables are based on pair or neighbour counting, which reduces the choice of algorithm to determining which counts to exclude for each jackknife subsample. For $w_{\mathrm{p}}(r_{\mathrm{p}})$, we discard a pair if either or both of the galaxies reside in the excluded cell. For $\Delta\Sigma(r_{\mathrm{p}})$, we only calculate the dark matter density profile around galaxies that live outside the excluded cell. Note, though, that in doing so we include dark matter particles that lie in that cell. Excluding such particles would lead to anomalous density profiles that are not easily corrected because the subsampling procedure violates the periodicity of the simulation volume. Similarly, for $\mathrm{VPF}(r)$, we place random spheres about points outside of the excluded cell; however, for the purposes of determining whether or not a particular sphere is a void region, galaxies within the excluded cell are taken into account. And for the count statistics, $P(N_{\mathrm{CIC}})$, $P(N_{\mathrm{CIA}})$, and $P(N_{2}/N_{5})$, we only centre cylinders on galaxies outside the excluded cell, but include companion galaxies within the excluded cell in our counts.

Of the three contributions, $\mathbf{C}_{\mathrm{jackknife}}$ is the dominant component, and $\mathbf{C}_{\mathrm{random}}$ is negligibly small, suggesting that we have used sufficiently large samples of VPF centres and dark matter particle positions to render the noisiness induced by finite sampling of these distributions negligibly small.

Direct inversion of the covariance matrix $\mathbf{C}$ is problematic numerically. Briefly stated, the uncertainty in the covariance will lead to the smallest eigenvalues of $\mathbf{C}$ being dominated by noise. When a matrix is inverted, its eigenvalues are inverted, which implies that the small, noisy eigenvalues of $\mathbf{C}$ become the large and noisy eigenvalues of its inverse. The inverse hence becomes dominated by this noise. This problem is further compounded by the large differences among the matrix elements inherited from the differences between the natural scales of the different observables, leading to extremely large differences in the sizes of the matrix eigenvalues. For example, the natural scale of the two-point function, $w_{\mathrm{p}}(r_{\mathrm{p}})$, is $\sim 10^{2}$, whereas the natural scale of $P(N_{\mathrm{CIC}})$ is on the order of $\sim 10^{-2}$ (see Fig. 1) and this difference leads to very different covariance matrix elements.

A common approach when faced with this problem is to truncate the smallest eigenvalues of $\mathbf{C}$, and calculate the Moore-Penrose pseudo-inversion (Penrose, 1955). We therefore normalize $\mathbf{C}$, and perform pseudo-inversion, excluding the eigenvalues of the correlation matrix smaller than $10^{-5}$ times the largest eigenvalue. Choosing other reasonable values of this cutoff does not substantially impact our results. In the future, when larger data sets are available, and the noise level of the covariance is sufficiently low, pseudo-inversion may no longer be necessary.

Our covariance matrices have several noteworthy features. Firstly, from an inspection of Fig. 2, we see strong self correlation between the bins of $w_{\mathrm{p}}(r_{\mathrm{p}})$, $\Delta\Sigma(r_{\mathrm{p}})$, and $\mathrm{VPF}(r)$ over a wide range of scales. On the contrary, $P(N_{2}/N_{5})$ shows weak correlation among its bins and with other observables, as it measures the distribution of the dimensionless ratio $N_{2}/N_{5}$, and is insensitive to the cosmic variance of galaxy number density. Secondly, for $P(N_{\mathrm{CIC}})$ and $P(N_{\mathrm{CIA}})$, the probability of smaller counts and larger counts sum up to unity, and are anti-correlated by construction, producing the sign reversal in the corresponding matrix blocks. Additionally, when comparing $\mathbf{C}_{\mathrm{jackknife}}$ and $\mathbf{C}_{\mathrm{realization}}$, it is obvious that the observable values are more correlated among jackknife subsamples than stochastic realizations. The sign of correlation coefficients approximately coincide between the two contributions, with the exception of blocks involving $\mathrm{VPF}(r)$. In $\mathbf{C}_{\mathrm{realization}}$, $\mathrm{VPF}(r)$ has a weak positive correlation with $w_{\mathrm{p}}(r_{\mathrm{p}})$ and $\Delta\Sigma(r_{\mathrm{p}})$ in most of the bins, while in $\mathbf{C}_{\mathrm{jackknife}}$, $\mathrm{VPF}(r)$ is anti-correlated with $w_{\mathrm{p}}(r_{\mathrm{p}})$ at larger $r_{\mathrm{p}}$ and $\Delta\Sigma(r_{\mathrm{p}})$. These are non-trivial effects, as $\mathrm{VPF}(r)$ is dependent on multiple moments of the galaxy number density field. The jackknife subsamples probe different regions of the box, with denser regions corresponding to stronger galaxy–galaxy and galaxy–matter correlation as well as fewer voids, leading to the anti-correlation in $\mathbf{C}_{\mathrm{jackknife}}$. On the other hand, among different realizations, higher values of $w_{\mathrm{p}}(r_{\mathrm{p}})$ and $\Delta\Sigma(r_{\mathrm{p}})$ result not from higher galaxy number densities but when more galaxies are concentrated in clusters, allowing more voids to exist in the rest of the space, giving rise to a positive correlation with $\mathrm{VPF}(r)$. Covariances for other luminosity samples have qualitatively similar features.

When applying the same analyses to galaxy survey data, observational uncertainties need to be taken into account, here we discuss how our covariances compare to observational covariances from SDSS data. For statistics that depend only on galaxy distribution, e.g., all of our observables except $\Delta\Sigma(r_{\mathrm{p}})$, the observational uncertainties depend on the survey volume and target number density to first order, both of which we mimic in our analyses. In comparison with the covariances from SDSS, which has a similar volume for the $M_{r}<-20.0$ sample to the Bolshoi Planck simulation, indeed our jackknife covariance matrix for $w_{\mathrm{p}}(r_{\mathrm{p}})$ is comparable to the measurement in Zehavi et al. (2011). For the fainter samples, SDSS has smaller volumes and hence larger covariances than ours, and vice versa for the brighter samples. Similarly, we expect the other observables to have covariances comparable to SDSS observation, with the exception of $\Delta\Sigma(r_{\mathrm{p}})$. Aside from the first-order comparability, there are other factors that need to be accounted for in observation. Fibre collision (see, e.g., O’Mill et al., 2012) affects pairs with small separations, i.e., the smallest bins of $w_{\mathrm{p}}(r_{\mathrm{p}})$ and $\Delta\Sigma(r_{\mathrm{p}})$, as well as $P(N_{\mathrm{CIC}})$ and therefore $P(N_{2}/N_{5})$. This effect needs to be forward modelled in survey data analyses. However, $P(N_{\mathrm{CIA}})$ does not depend on the closest galaxy companions, and is robust to fibre collision. As we will show in Section 3, $P(N_{\mathrm{CIA}})$ typically provides tight constraints comparable to $P(N_{\mathrm{CIC}})$. Some of the other second-order observational uncertainties are blending and saturation in target selection, redshift measurement errors, and geometric features, all of which we expect to be subdominant in this context.

On the other hand, the covariances for $\Delta\Sigma(r_{\mathrm{p}})$ measured for SDSS data (Mandelbaum et al., 2006b) are significantly larger than our covariances. This is expected because the survey data is dominated by shape noise, while we neglect shape noise in our study, assuming infinite source densities. Taking non-zero shape noise into account will result in weaker constraints from lensing than those that we find here, and our forecasts therefore must be regarded as the upper limit of constraining power that can be achieved with $\Delta\Sigma(r_{\mathrm{p}})$.

2.6 Derivative Fitting

In order to compute the elements of the Fisher matrix, it is necessary to estimate the partial derivatives of the observables with respect to the parameters in the neighbourhood of the fiducial point in the parameter space. We designated these derivatives as $\partial\mathbf{f}/\partial\mathbf{p}$ in Eq. (11) above. Assessing these derivatives from realizations of the perturbed models is non-trivial due to the inherent stochasticity of using mock catalogues based upon the direct population of N-body simulations. Therefore, we give a detailed description of our approach to estimating derivatives in Appendix A.

3 Results

We perform Fisher matrix analyses in order to forecast the constraints on the dHOD model that can be extracted from combinations of the galaxy observables described in Section 2.3. In Table 1, we list the fiducial HOD parameters corresponding to 5 galaxy samples selected by luminosity. In this section, we present our primary results in terms of estimated posterior 1$\sigma$ constraints for the four lower luminosity samples that we have studied. We exclude the brightest, $M_{r}<-21.0$ sample from our primary results because the results from this sample are subject to excessive statistical fluctuations due to the relatively small number of galaxies above this luminosity threshold within the volume of the Bolshoi Planck simulation. For completeness, the results from the $M_{r}<-21.0$ sample are included in our comprehensive list of results in Appendix B.

Before proceeding to our results, we note that the dominant contribution to our errors are from sample variance due to the finite volume of the Bolshoi Planck simulation (see Section 2.5 above). The absolute constraints on parameters will decrease with increasing volume. As such, we focus on relative constraints on parameters from different combinations of observables, rather than on the absolute values. Our study requires that haloes be resolved with a very large number of particles which, in turn, stipulates the use of a high-resolution, relatively smaller volume simulation such as Bolshoi Planck.

3.1 Assessing the Complementarity of Observables

We begin with a discussion of our forecast constraints from individual observables and combinations of any two observables. We use the $M_{r}<-20.0$ sample as an example in Fig. 3, to compare the constraints on $A_{\mathrm{cen}}$ and $A_{\mathrm{sat}}$ from the individual observables and all the possible combinations of two observables. We caution that each of the four panels has a different $y$-axis range, which the reader must account for when comparing results among different panels.

We display constraints from the individual observables on $A_{\mathrm{cen}}$ in the upper left panel of Fig. 3 and $A_{\mathrm{sat}}$ in the lower left panel of Fig. 3. The columns are ordered according to increasing constraining power, and the filled circles indicate the observables from which the constraints were derived. For example, the leftmost column in the upper left panel of Fig. 3 shows the constraint on $A_{\mathrm{cen}}$ derived from $P(N_{\mathrm{CIC}})$.

The constraints from individual observables, displayed in the left panels of Fig. 3 convey several points. The observable $P(N_{\mathrm{CIC}})$ more strongly constrains the assembly bias of satellite galaxies than that of central galaxies. This is expected because $P(N_{\mathrm{CIC}})$ primarily probes the “one-halo term”, and is sensitive to the satellite population that accompany centrals. For the brighter samples with higher satellite fractions, $P(N_{\mathrm{CIA}})$ is more dependent on the satellite population, and constrains $A_{\mathrm{sat}}$ more strongly, while its constraining power decreases for $A_{\mathrm{cen}}$. The constraints from $w_{\mathrm{p}}(r_{\mathrm{p}})$ are dominated by the smaller radial bins, which have higher signal to noise ratios than the measurements at larger scales. This causes the absolute constraints on $A_{\mathrm{cen}}$ to be weaker than those on $A_{\mathrm{sat}}$. However, with larger volumes, e.g., DESI (Levi et al., 2013), large-scale clustering will be measured with higher precision, enabling better constraints on $A_{\mathrm{cen}}$. We also find that $\mathrm{VPF}(r)$ gives strong constraints on $A_{\mathrm{cen}}$, but is extremely inefficient in constraining $A_{\mathrm{sat}}$. This can be explained by the fact that a single galaxy suffices to eliminate the possibility that a region could be a void. The vast majority of satellite galaxies reside in haloes where there are central galaxies, which already eliminate the void, therefore the void probability function is largely insensitive to the abundances of satellite galaxies. For this reason, $\mathrm{VPF}(r)$ is a poor probe of not only $A_{\mathrm{sat}}$, but all of the parameters that determine satellite populations (see Tables 3-6 in Appendix B).

The constraints from the combinations of two observables are shown in the right-hand panels of Fig. 3. The columns are again ordered from least constraining to most constraining, and the observables used in each analysis are marked by filled circles. For example, the leftmost column in the upper right panel of Fig. 3 shows constraints derived from the combination of $w_{\mathrm{p}}(r_{\mathrm{p}})$ and $\Delta\Sigma(r_{\mathrm{p}})$. The combination of $w_{\mathrm{p}}(r_{\mathrm{p}})$ and $\Delta\Sigma(r_{\mathrm{p}})$ has been used in a number of previous studies; however, we find that this is one of the least constraining of the combinations that we have considered for both $A_{\mathrm{cen}}$ and $A_{\mathrm{sat}}$. As we discuss further below, this is chiefly because $w_{\mathrm{p}}(r_{\mathrm{p}})$ and $\Delta\Sigma(r_{\mathrm{p}})$ share largely common degeneracies among the (d)HOD parameters, so that combining the two does not result in significant decreases in marginalized constraints that can be expected when combining highly complementary data.

The primary result to be gleaned from Fig. 3 is the overall efficacy of the count-based observables, i.e., $P(N_{\mathrm{CIC}})$ and $P(N_{\mathrm{CIA}})$, to complement either $w_{\mathrm{p}}(r_{\mathrm{p}})$ or $\Delta\Sigma(r_{\mathrm{p}})$ to constrain the galaxy assembly bias parameters. We find that when used in combination with count-based statistics, $w_{\mathrm{p}}(r_{\mathrm{p}})$ typically outperforms $\Delta\Sigma(r_{\mathrm{p}})$. Moreover, it is worth noting that our lensing covariance assumes an infinite density of lensing sources, so lensing constraints realized from a real survey analysis will be further diluted by shape noise contributions to the covariance, as we have discussed in Section 2.5. Therefore, we suspect that this general result will be robust to actual survey analyses. While Fig. 3 displays only constraints on $A_{\mathrm{cen}}$ and $A_{\mathrm{sat}}$ from a single luminosity threshold sample, we find that these qualitative results hold for all thresholds samples considered here (see Fig. 4 and Appendix B).

3.2 Complementarity with Clustering

We now turn to a more detailed exploration of the complementarity of various observables with galaxy clustering, as quantified by $w_{\mathrm{p}}(r_{\mathrm{p}})$. Figure 4 displays the constraining power of different combinations of observables on the two galaxy assembly bias parameters of our dHOD model: $A_{\mathrm{cen}}$, the central galaxy dHOD assembly bias parameter (top panel), and $A_{\mathrm{sat}}$, the satellite galaxy dHOD assembly bias parameter (bottom panel). We include in Fig 4 constraints from $w_{\mathrm{p}}(r_{\mathrm{p}})$ individually along with constraints from combining $w_{\mathrm{p}}(r_{\mathrm{p}})$ with each of the other observables that we study. For completeness, we also show the constraints from all the observables combined as an illustration of the maximal constraining power that can be achieved using the complete set of observables considered in our study. The constraints in Fig. 4 are depicted as bar plots, with the bars grouped by combination of observables. Bars of different colours within each group correspond to different luminosity threshold samples, as indicated. Finally, the heights of the bars represent the fully marginalized 1$\sigma$ constraints, with smaller values corresponding to tighter, more restrictive constraints. Similar plots for the other dHOD parameters can be found in Appendix B.

Examining the bars in Fig. 4, several general trends are apparent. Most prominently, similar to what we find for the $M_{r}<-20.0$ sample in Fig. 3, for all the luminosity samples we study, when combined with $w_{\mathrm{p}}(r_{\mathrm{p}})$, $P(N_{\mathrm{CIC}})$ is generally more effective at constraining $A_{\mathrm{cen}}$ and $A_{\mathrm{sat}}$ than the commonly-used $\Delta\Sigma(r_{\mathrm{p}})$, despite our assumption of infinite lensing source density. In fact, in most cases the count statistics are the most effective observables to combine with clustering in an effort to constrain dHOD models of assembly bias. Furthermore, as is shown in Appendix B, this statement is typical of the constraints on most of the HOD parameters, especially for the fainter samples.

When comparing results for galaxy samples defined by different luminosity thresholds, the constraints are typically tighter for the fainter samples. This mainly reflects the fact that brighter samples have lower number densities, resulting in higher levels of noise. However, there are some exceptions. For example, the constraints on $A_{\mathrm{cen}}$ from $w_{\mathrm{p}}(r_{\mathrm{p}})$ alone are tighter for the $M_{r}<-19.5$ sample than for the $M_{r}<-19.0$ sample. In these cases, the degeneracies among different parameters depend on luminosity, such that after marginalization over all other parameters the noiser, high luminosity sample yields tighter constraints. The unmarginalized constraints are all monotonically increasing functions of luminosity threshold, as they must be.

3.3 Marginalized Two-Dimensional Constraints

In addition to fully marginalized constraints, it is interesting to examine parameter constraints in two-dimensional subspaces of the full parameter space. In Figure 5, we plot the marginalized 1$\sigma$ contours in each of the 2D projections of our 6-dimensional parameter space. As we utilize a Fisher matrix to estimate parameter constraints, all contours are elliptical and are centred around the fiducial point in the parameter space. Fig. 5 corresponds to the $M_{r}<-20.0$ threshold sample, but the other samples have qualitatively similar features. Different contours correspond to different observables, or combinations thereof, as indicated, and we have highlighted the results for $w_{\mathrm{p}}(r_{\mathrm{p}})+\Delta\Sigma(r_{\mathrm{p}})$ and $w_{\mathrm{p}}(r_{\mathrm{p}})+P(N_{\mathrm{CIC}})$ using thicker contours. To avoid crowding, we use $P(N_{\mathrm{CIC}})$ as the representative case for the various count statistics. For comparison, the gray shaded ellipse shows the constraints derived from exploiting all of our observables simultaneously.

From Fig. 5, it is apparent that in all projections, $w_{\mathrm{p}}(r_{\mathrm{p}})+P(N_{\mathrm{CIC}})$ (thick red lines) is superior to the other combinations of observables, particularly the commonly-used combination of $w_{\mathrm{p}}(r_{\mathrm{p}})+\Delta\Sigma(r_{\mathrm{p}})$. This result remains strictly true for the $M_{r}<-19.0$ and $M_{r}<-19.5$ samples, while for the brighter, $M_{r}<-20.5$ sample, a few parameter combinations are more tightly constrained using parameter combinations other than $w_{\mathrm{p}}(r_{\mathrm{p}})+P(N_{\mathrm{CIC}})$.

Figure 6 displays confidence contours for our model parameters constrained by $w_{\mathrm{p}}(r_{\mathrm{p}})$, $\Delta\Sigma(r_{\mathrm{p}})$ and $P(N_{\mathrm{CIC}})$ individually (rather than combinations of observables). This visualization can aid in the qualitative understanding of our results. Combining $w_{\mathrm{p}}(r_{\mathrm{p}})$ with $\Delta\Sigma(r_{\mathrm{p}})$ yields limited improvement because both of these observables share similar degeneracy directions in multiple dimensions of the parameter space. The combination of $w_{\mathrm{p}}(r_{\mathrm{p}})$ with $P(N_{\mathrm{CIC}})$ is superior because these observables have largely complementary degeneracy directions in the parameter space and combining these observables leads to the simultaneous breaking of multiple degeneracies. Notice that constraints from $P(N_{\mathrm{CIC}})$ on any single parameter are not particularly restrictive; however, the constraints from $P(N_{\mathrm{CIC}})$ restrict parameter values to exceedingly narrow degeneracy regions, which, in turn, leads to significant improvements in constraining power when combined with $w_{\mathrm{p}}(r_{\mathrm{p}})$.

3.4 Constraints on Parameters of the Standard HOD

In the previous subsections, we focused on constraints on assembly bias parameters in the dHOD model and showed that counts-in-cylinders is an effective comlement to the galaxy projected two-point function for diagnosing and constraining assembly bias. However, it is also interesting to study constraints on the standard HOD parameters in a standard HOD model that does not include assembly bias.

Figure 7 depicts 2D marginalized, projected constraint contours on the standard HOD parameters from an analysis to constrain a standard HOD model using the various observables that we consider. It is clear that the complementarity of $P(N_{\mathrm{CIC}})$ extends to the parameters of the standard HOD, as the combination of $w_{\mathrm{p}}(r_{\mathrm{p}})$ with $P(N_{\mathrm{CIC}})$ outperforms the combination of $w_{\mathrm{p}}(r_{\mathrm{p}})$ with $\Delta\Sigma(r_{\mathrm{p}})$ in all projections. This strongly suggests that $P(N_{\mathrm{CIC}})$ is a favourable observable even in studies using standard HOD.

3.5 Limitations and Caveats

Our results are subject to several limitations and caveats that we describe in this subsection. First, all Fisher matrix analyses are based on a linear expansion of the likelihood: The observable–parameter relation is treated as a linear function of the observable on the parameter. This linearity is not an accurate model of the observable–parameter relation over the entirety of the relevant domain of parameters for all of the observable–parameter combinations that we explore. As a result, the derivative values $\partial\mathbf{f}/\partial\mathbf{p}$, and therefore the constraints, depend on the choice of the fiducial point in the parameter space. For the results presented above, the fiducial point is motivated by previous data analyses using a standard HOD, specifically, the study of Zentner et al. (2019). The fiducial values of the assembly bias parameters, $A_{\mathrm{cen}}$ and $A_{\mathrm{sat}}$ are set to zero. We address this particular caveat in the following subsection, pointing out the dependence of our conclusions upon the fiducial location in the parameter space.

Second, the decorated HOD parameters $A_{\mathrm{cen}}$ and $A_{\mathrm{sat}}$ can only vary over the interval $[-1,1]$ because the degree of galaxy assembly bias that is possible is limited (see Hearin et al., 2016, for details). The Fisher formalism assumes a multivariate Gaussian posterior distribution, so such hard boundaries on the parameter space can lead to gross violations of this assumption. In particular, any time that $\sigma_{A_{\mathrm{cen}}}$ or $\sigma_{A_{\mathrm{sat}}}$ approach unity (or even exceed it, see Appendix B), our estimates will not be a reliable, quantitative estimate of the constraining power of the observables. However, the qualitative comparisons among observables should not be impacted by this shortcoming and the Fisher matrix will still give a reliable ranking of the relative utility of different combinations of observables.

An additional caveat to our results is associated with the particular model that we explore. While we phrase our results qualitatively in terms of constraints on galaxy assembly bias, it is important to realize that our calculations pertain only to a specific model, namely the dHOD with a binary split on galaxy populations. It is possible that our conclusions would change significantly if a different halo occupation model is used. Examples of different models might include a standard HOD with an augmented set of parameters or a wholly different model for the galaxy–halo relationship, such as the conditional luminosity function (e.g., Yang et al., 2003) or a parameterized form of abundance matching (e.g., Lehmann et al., 2017). There are a limited number of models that include tunable galaxy assembly bias based on halo properties. Aside from the dHOD, Lehmann et al. (2017) parametrizes the dependence of galaxy luminosity on halo concentration at a fixed halo mass in the abundance matching model with an interpolation scheme, and Tinker et al. (2008a), Wibking et al. (2019) and McEwen & Weinberg (2018) allow for local density-dependent variations of the HOD. We limit our results to the dHOD model and assembly bias due to concentration, and relegate more comprehensive studies of galaxy assembly bias to future work.

Our constraints from the galaxy–galaxy weak lensing signal $\Delta\Sigma$ represent the upper limit of information that can be gained in a real data analysis, because we use a weak lensing covariance that includes only sample variance. This is equivalent to assuming an infinite background source galaxy density and thus a shape-noise-free measurement of the galaxy–galaxy lensing signal. Since the galaxy–galaxy lensing signal presented is optimistic, our primary qualitative result, namely that the combination of $w_{\mathrm{p}}(r_{\mathrm{p}})$+$P(N_{\mathrm{CIC}})$ is superior to combinations that include the galaxy–galaxy lensing signal, will not be affected when complete galaxy–galaxy lensing covariances are used.

As a final caveat, we emphasize that our work uses only a single simulation and thus, we work in the context of a single set of cosmological parameters. In particular, we consider the best-fit Planck cosmology used as the modelling framework in the Bolshoi Planck simulation. This limitation is difficult to circumvent at this time due to the need for simulations that are both large volume (to model clustering) and high resolution (to measure the internal properties of haloes). The computational costs of such simulations prohibit simultaneous explorations of assembly bias and cosmology within the scope of this paper. However, important steps are being taken in precisely this direction (e.g., Zhai et al., 2019). Nonetheless, the Planck constraints on cosmological parameters are quite restrictive (Planck Collaboration et al., 2018a, b) and we do not expect modifications to the cosmological model to have a significant impact on our qualitative results.

3.6 Dependence on Fiducial Parameters

One of the caveats mentioned in the previous section is that Fisher analyses yield results that may depend upon the fiducial point in the parameter space about which the likelihood is expanded. Another way to say this is that the constraints depend upon the point in parameter space that corresponds to the true underlying model. In the results we presented above, we assumed that the true fiducial model corresponded to zero galaxy assembly bias ($A_{\mathrm{cen}}=A_{\mathrm{sat}}=0$).

To examine how our results depend on this choice of the fiducial model, we now repeat our analysis for two alternative assumptions for the fiducial values of the galaxy assembly bias parameters.

In the first, we adopt a fiducial central galaxy assembly bias of $A_{\mathrm{cen}}=0.5$, which is motivated by the recent analysis of galaxy clustering by Zentner et al. (2019), while keeping $A_{\mathrm{sat}}=0$. In this case, we find results that are both qualitatively and quantitatively similar to our main model with $A_{\mathrm{cen}}=A_{\mathrm{sat}}=0$.

We have also explored the dependence of our forecast constraints on the underlying amount of satellite assembly bias. To do so, we repeated our analyses with the fiducial satellite assembly bias parameter set to $A_{\mathrm{sat}}=-0.6$. This value of $A_{\mathrm{sat}}$ has several motivations. First, Zentner et al. (2019) showed that clustering of galaxies in the SDSS $M_{r}<-19.5$ threshold sample is consistent with significantly negative values of $A_{\mathrm{sat}}$. Furthermore, it is known that the abundance of dark matter subhaloes is anti-correlated with host halo concentration (e.g., Zentner et al., 2005; Mao et al., 2015; Jiang & van den Bosch, 2017). Since subhaloes are believed to host satellite galaxies, this anti-correlation implies a negative value for $A_{\mathrm{sat}}$.

The results of this experiment are shown in Fig. 8. While the results for $A_{\mathrm{cen}}$ (upper panel) are similar to the case of our main model with $A_{\mathrm{cen}}=A_{\mathrm{sat}}=0$ (i.e., $A_{\mathrm{cen}}$ is always best constrained by the combination of $w_{\mathrm{p}}(r_{\mathrm{p}})$ plus $P(N_{\mathrm{CIC}})$), the results for $A_{\mathrm{sat}}$ (lower panel) are notably different. In particular, the combination of $w_{\mathrm{p}}(r_{\mathrm{p}})$ and $\Delta\Sigma(r_{\mathrm{p}})$ now yields the tightest constraints on $A_{\mathrm{sat}}$, rather than the weakest. This very tight constraint stems from two things. First, $w_{\mathrm{p}}(r_{\mathrm{p}})$ on small scales ($\lesssim 1\,h^{-1}\mathrm{Mpc}$) has a much stronger dependence on $A_{\mathrm{sat}}$ near $A_{\mathrm{sat}}=-0.6$. This improves constraints from $w_{\mathrm{p}}(r_{\mathrm{p}})$ alone, reducing parameter degeneracy. Second, the observables $w_{\mathrm{p}}(r_{\mathrm{p}})$ and $\Delta\Sigma(r_{\mathrm{p}})$ are more complementary to one another near $A_{\mathrm{sat}}=-0.6$ because the degeneracy directions selected by the $w_{\mathrm{p}}(r_{\mathrm{p}})$ constraints change their orientation slightly in the parameter space. These improvements jointly boost the constraining power of the combination of $w_{\mathrm{p}}(r_{\mathrm{p}})$ and $\Delta\Sigma(r_{\mathrm{p}})$. However, we must note that several observational and theoretical factors that will likely impact the constraining power of very-small-scale clustering statistics are not accounted for in our analysis. Chief among these omissions is our neglect of shape noise in galaxy–galaxy lensing covariances. Including shape noise is likely to reduce significantly the complementarity of $\Delta\Sigma(r_{\mathrm{p}})$ with other probes, including $w_{\mathrm{p}}(r_{\mathrm{p}})$. Modelling uncertainties, such as the choice of the radial distributions of satellite galaxies will also reduce the constraining power of observables on small-scales with similar result. Hence, the tight constraints on $A_{\mathrm{sat}}$ we observed here are likely to be quite optimistic.

To summarize, the results that we have shown so far point toward a clear conclusion: Count-based galaxy clustering statistics, such as counts-in-cells distributions ($P(N_{\mathrm{CIC}})$), can be instrumental in constraining galaxy assembly bias, and are particularly powerful in constraining the galaxy assembly bias of central galaxies. For reference, we tabulate the forecast constraints on all dHOD parameters from all of the observable combinations we study, and for all five luminosity threshold samples in Appendix B.

4 Discussion and Conclusions

Constraining galaxy assembly bias is important to the study of the connection between galaxies and haloes and for extracting the maximum possible information on both galaxy evolution and cosmology from survey data. Numerous studies use galaxy clustering to constrain either the galaxy–halo connection or cosmology or both (e.g., Hawkins et al., 2003; van den Bosch et al., 2003; Yang et al., 2004; Seljak et al., 2005, 2006; Cooray, 2006; van den Bosch et al., 2007; Zheng et al., 2007; Moster et al., 2010; Zehavi et al., 2011; Anderson et al., 2012; Guo et al., 2012a, b; Leauthaud et al., 2012; Cacciato et al., 2013; Reddick et al., 2013; Guo et al., 2014; Reid et al., 2014; Coupon et al., 2015; More et al., 2015; Saito et al., 2016; Lange et al., 2019; Cowley et al., 2018; Sinha et al., 2018; Xu et al., 2018). Several of these works combine clustering with either weak galaxy–galaxy lensing measurements or with measurements of redshift space distortions in order to constrain the galaxy–halo connection and/or cosmology, and this use of complementary variables is becoming increasingly common. However, these results may suffer from systematic bias when galaxy assembly bias is not properly included in the model (e.g., Zentner et al., 2014). To date, there are only a small number of studies constraining assembly bias with galaxy clustering data (Lehmann et al., 2017; Zentner et al., 2019; Vakili & Hahn, 2019).

We have studied the ability of various galaxy clustering statistics to constrain assembly bias in an effort to determine which combination(s) of observables are most informative. In particular, we have estimated the relative constraining power of several spatial galaxy clustering statistics to constrain the assembly bias parameters of the decorated halo occupation distribution (dHOD). In this first study of its kind, we have restricted our attention to statistics that are not particularly sensitive to galaxy peculiar velocities (however, we do work in redshift space so our results are not completely immune to peculiar velocities). We have chosen to do this because including peculiar velocities (such as redshift space distortions) opens up many distinct modeling questions. Exploring redshift space distortions will be the subject of future work.

In general, we have found that the combination of $w_{\mathrm{p}}(r_{\mathrm{p}})$ with the counts-in-cylinders distribution ($P(N_{\mathrm{CIC}})$) is the most effective combination of two observables for constraining assembly bias within the dHOD modeling framework. This combination outperforms the commonly-used combination of $w_{\mathrm{p}}(r_{\mathrm{p}})$ and $\Delta\Sigma(r_{\mathrm{p}})$ in constraining central galaxy assembly bias by a factor of $\sim 2$ for lower luminosity samples, and to a slightly lesser extent in constraining satellite galaxy assembly bias. The primary reason for this is that $w_{\mathrm{p}}(r_{\mathrm{p}})$ and $\Delta\Sigma(r_{\mathrm{p}})$ share roughly common degeneracy directions in the dHOD parameter space, while the combination $w_{\mathrm{p}}(r_{\mathrm{p}})$ and $P(N_{\mathrm{CIC}})$ is much more complementary. This implies that the combination of $\Delta\Sigma(r_{\mathrm{p}})$+$P(N_{\mathrm{CIC}})$ is nearly as good as $w_{\mathrm{p}}(r_{\mathrm{p}})$+$P(N_{\mathrm{CIC}})$, which we have confirmed. The complementarity between $\Delta\Sigma(r_{\mathrm{p}})$ and $P(N_{\mathrm{CIC}})$ is not unexpected. The lensing signal traces the matter density contrast around galaxies while the counts-in-cylinders statistics probe the galaxy distribution profiles in approximately the same regions, and are therefore expected to complement each other in constraining cosmology as well as the connection between galaxies and the matter field (Gruen et al., 2018; Friedrich et al., 2018).

We have shown that the count statistics are also effective in constraining the standard HOD parameters, independent of whether the actual halo occupation statistic are affected by galaxy assembly bias or not. In addition, we have tested different fiducial models and find that when strong satellite assembly bias is present, the combination of $w_{\mathrm{p}}(r_{\mathrm{p}})+\Delta\Sigma(r_{\mathrm{p}})$ actually provides the tightest constraints on $A_{\mathrm{sat}}$ among all sets of observables studied here. Therefore we caution that the preferred statistics may depend on the true, underlying relationship between galaxies and dark matter haloes, and the degree to which galaxy assembly bias is realized in nature.

Our results complement recent work studying the information that can be gained from higher-order statistics beyond the two-point function. In particular, in Yuan et al. (2017, 2018) it was shown that the three-point function in the squeezed limit contains significant additional constraining power on HOD parameters that is complementary to $w_{\mathrm{p}}(r_{\mathrm{p}}).$ This is consistent with our findings in the sense that the complementarity of counts statistics to $w_{\mathrm{p}}(r_{\mathrm{p}})$ derives precisely from the extraction of information in higher-order $n$-point moments of the density field. Considering the results in Yuan et al. (2018) together with our findings, statistics beyond two-point clustering and lensing should be seriously considered in future analyses of large-scale structure data that utilize models of the galaxy-halo connection.

In conclusion, our results strongly suggest that galaxy assembly bias may be significantly better constrained and/or better understood by employing simple counts statistics as measured from forthcoming and present-day data sets. There are numerous forthcoming data sets with which can be used to inform assembly bias and/or whose interpretation may be challenged by small levels of assembly bias. These include large redshift surveys, such as may be carried out by DESI (Levi et al., 2013) or WFIRST (Spergel et al., 2015). As we have already mentioned, interesting follow up work includes an exploration of velocity statistics, such as redshift space distortions, in redshift surveys. While our work relates specifically to redshift surveys, it would be interesting to explore possible avenues for studying assembly bias within photometric surveys, such as the DES (The Dark Energy Survey Collaboration, 2005; DES Collaboration et al., 2018) and LSST (LSST Science Collaboration et al., 2009). It is our aim to study and deploy these statistics to constrain the galaxy–halo connection and to encourage others to do the same.

Acknowledgements

We thank Brian Batell, Arthur Kosowsky, Rachel Mandelbaum, Dazhuang Mao, Jeffrey Newman, Troy Raen, Li-Cheng Tsai, and Rongpu Zhou for useful discussions. KW, ARZ and ASV are supported by the US National Science Foundation (NSF) through grant AST 1517563 and by the Pittsburgh Particle Physics Astrophysics and Cosmology Center (PITT PACC). ASV is additionally supported by the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357. YYM is supported by PITT PACC through the Samuel P. Langley Postdoctoral Fellowship. FvdB and JUL are supported by the US National Science Foundation (NSF) through grant AST 1516962. FvdB received additional support from the Klaus Tschira foundation, and from the National Aeronautics and Space Administration through Grant No. 17-ATP17-0028 issued as part of the Astrophysics Theory Program. CMS is supported by NSF through grant DMS 1521786. APH is supported under the DOE contract DE-AC02-06CH11357. DC is supported by a McWilliams Postdoctoral Fellowship.

The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) and the Partnership for Advanced Supercomputing in Europe (PRACE, www.prace-ri.eu) for funding the MultiDark simulation project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, www.lrz.de). The Bolshoi simulations have been performed within the Bolshoi project of the University of California High-Performance AstroComputing Center (UC-HiPACC) and were run at the NASA Ames Research Center.

This research made use of Python, along with many community-developed or maintained software packages, including IPython (Pérez & Granger, 2007), Jupyter (jupyter.org), Matplotlib (Hunter, 2007), NumPy (van der Walt et al., 2011), and SciPy (Jones et al., 2001). This research made use of NASA’s Astrophysics Data System for bibliographic information.

Appendix A Derivative Fitting

In this appendix, we describe our approach to obtain the derivative matrix $\partial\mathbf{f}/\partial\mathbf{p}$.

We begin from the fiducial parameter set specified in Table 1, and add perturbations to each dHOD parameter (except $M_{0}$, which we do not vary in our analysis) in turn, keeping all other parameters fixed, to examine the dependence of observables on each individual parameter (i.e., the set of $\mathbf{f}$–$\mathbf{p}$ relations). The process is non-trivial because we construct our observables from mock catalogues based on the population of a simulation of finite size with galaxies drawn from the (d)HOD. The fact that we construct our observables using mock catalogues has the advantage of accuracy compared to analytic approximation methods, but it also has the unfortunate consequence that individual observables can exhibit non-negligible stochasticity from one mock catalogue realization to another. In order to mitigate the impact of this variability, rather than taking $\Delta\mathbf{f}/\Delta\mathbf{p}$ at a single perturbed point, we consider a series of perturbed models in the neighbourhood of the fiducial value, along each dimension of the parameter space. We then fit our set of $\mathbf{f}$–$\mathbf{p}$ relations for a slope, as an estimate of the partial derivative in the neighbourhood of the fiducial point. We will elaborate on these procedures in the remainder of this appendix.

As we have discussed in Section 2.2.2, in the Halotools implementation of the (d)HOD model, the galaxy occupation of each halo is randomly drawn from the probability distribution function determined by the properties of the halo. The mean occupation varies with the (d)HOD parameters, whereas any particular realization of the mock galaxy catalogue is also dependent on the sequence of random numbers used in this process. To generate the mock catalogues that are suitable for Fisher analysis, we use a fixed random seed for each random variable in the (d)HOD model. The values of random seeds are set independently of the (d)HOD parameter values. In this fashion, the number of galaxies in each halo will always have the same random deviate each time the halo is populated from the underyling (d)HOD model, regardless of the parameter values. Had we not implemented this, the $\mathbf{f}$–$\mathbf{p}$ relations we measure would be severely impacted by the significantly greater stochasticity introduced by the random process. It would then be required to construct a very large number of mock catalogues to marginalize over this stochasticity. By employing a common seed, the differences between these catalogues become primarily attributable to parameter differences and not statistical fluctuations from catalogue to catalogue that could be induced by finite sampling of the (d)HOD models. This exercise minimizes the variation of the observables attributable to stochasticity; however, $\mathbf{f}$ is generally not a smooth function of $\mathbf{p}$, due to the intricate nature of the dependence of the observables on the galaxy distribution.

Determining the $\mathbf{f}$–$\mathbf{p}$ relations from a single random number seed is not sufficient because the $\mathbf{f}$–$\mathbf{p}$ relation has a small dependence upon the random number seed. Consequently, the values of the resulting derivatives will vary slightly with different random seeds. For this reason, we repeat this entire process for a large number (or order $\sim 100$) of different random number seeds and take the trimmed mean of $\mathbf{f}$, averaging only the central $68\%$ values. In this manner, we construct smooth $\mathbf{f}$–$\mathbf{p}$ relations that do not dependent on the choice of the random seed. In the example shown in Fig. 9, we assess the relation between the projected two-point correlation function $w_{\mathrm{p}}$ at $r_{\mathrm{p}}=1.74h^{-1}\mathrm{Mpc}$ and the central galaxy assembly bias parameter $A_{\mathrm{cen}}$. Thin coloured lines correspond to different random number seeds, and for clarity we only plot results for a subset of all random seeds used. The thicker black line shows the trimmed mean from the central 68% of the fixed-seed $\mathbf{f}$–$\mathbf{p}$ data points, which serves as our estimate of the observable–parameter relation. This procedure provides us with a set of $\mathbf{f}$–$\mathbf{p}$ relations, one for each observable–parameter pair, that we fit as described below.

We use the R (R Core Team, 2016) package locfit (Loader, 2013) to fit local linear derivatives to the $\mathbf{f}$–$\mathbf{p}$ relation that we have obtained. We choose the degree of local polynomials to be two, which captures the shape of the curve without excessive overfitting. The locfit package provides 5 commonly used weighting kernels. We have confirmed that different kernels yield similar results. We present results obtained with the default tricube kernel. Because some of the $\mathbf{f}$–$\mathbf{p}$ relations are strongly non-linear, the smoothing scale for the local fit needs to be chosen with care. We choose the smoothing scale of fitting following the principle that the range of the parameter considered for the local derivative fit should be comparable with the posterior 1$\sigma$ constraint for each observable combination respectively. Qualitatively, this is motivated by the fact that the 1$\sigma$ constraints defines what it means to be in the “neighbourhood” of the fiducial point. To this criterion, we add two additional restrictions:

For each $\mathbf{f}$-$\mathbf{p}$ relation, we obtain the optimal smoothing parameter from generalized cross validation (GCV), using the loess.as function in the R package fANCOVA (Wang, 2010). We use this smoothing scale as a lower limit, in order to avoid fitting numerical noise. 2. 2.

We apply another lower limit of smoothing scale for each parameter in each threshold, based on the physical interpretation of $\mathbf{f}$-$\mathbf{p}$ relations, as some of the unphysical effects from the mock are not recognized by GCV. We have tested that our conclusions are not sensitive to this choice within a reasonable range.

In the example of Fig. 9, we fit a local linear derivative to the relation at the fiducial parameter set, marked by the vertical grey dashed line, i.e., $[\partial w_{\mathrm{p}}(r_{\mathrm{p}}=1.74h^{-1}\mathrm{Mpc})/\partial A_{\mathrm{cen}}]|_{A_{\mathrm{cen}}=0}$. The smoothing scale is shown as a grey band that is symmetric around the fiducial parameter. The local linear fit is shown by the solid red line. The fitted slopes are the derivatives we use in our forecasts.

Appendix B Forecast Constraints on HOD Parameters

In this appendix, we present a comprehensive list of our results for constraints on dHOD parameters.

In Fig. 10 and Fig. 11, we show bar plots of constraints on the decorated HOD parameters besides $A_{\mathrm{cen}}$ and $A_{\mathrm{sat}}$ (shown in Figure 4), excluding $\log M_{0}$ which is poorly constrained. These include $\alpha$, $\log M_{1}$, $\log M_{\mathrm{min}}$, and $\sigma_{\mathrm{log}M}$, which are also the original parameters of the standard HOD.

We also list the posterior constraint values in Tables 3, 4, 5, 6, and 7, for the 6 parameters that we allow to vary and for all 5 of the luminosity samples we study. The constraints we include are from all the individual observables, all the possible combinations of two observables, and the combination of all 6 observables.