{'fraction_non_alphanumeric': 0.031696525712627835, 'fraction_numerical': 0.016037652216661428, 'mean_word_length': 4.498958564148017, 'pattern_counts': {'":': 1, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Toward Accurate Modeling of Galaxy Clustering on Small Scales: Constraining the Galaxy-Halo Connection with Optimal Statistics

Applying halo models to analyze the small-scale clustering of galaxies is a proven method for characterizing the connection between galaxies and their host halos. Such works are often plagued by systematic errors or are limited to clustering statistics which can be predicted analytically. In this work, we employ a numerical mock-based modeling procedure to examine the clustering of SDSS DR7 galaxies. We apply a standard halo occupation distribution (HOD) model to dark-matter-only simulations with a LCDM cosmology. To constrain the theoretical models, we utilize a combination of galaxy number density and selected scales of the projected correlation function, redshift-space correlation function, group multiplicity function, average group velocity dispersion, mark correlation function, and counts-in-cells statistics. We design an algorithm to choose an optimal combination of measurements that yields tight and accurate constraints on our model parameters. Compared to previous work using fewer clustering statistics, we find significant improvement in the constraints on all parameters of our halo model for two different luminosity-threshold galaxy samples. Most interestingly, we obtain unprecedented high-precision constraints on the scatter in the relationship between galaxy luminosity and halo mass. However, our best-fit model results in significant tension (>4 sigma) for both samples, indicating the need to add second-order features to the standard HOD model. To guarantee the robustness of these results, we perform an extensive analysis of the systematic and statistical errors in our modeling procedure, including a first-of-its-kind study of the sensitivity of our constraints to changes in the halo mass function due to baryonic physics.


INTRODUCTION
One of the great strengths of ΛCDM is in its accurate predictions of the large-scale density field of the Universe today. It is common to quantify this density field through statistical measurements of the clustering of galaxies. ΛCDM is highly successful at predicting galaxy clustering on large physical scales where galaxies are simple tracers of the underlying matter density field. On small scales ( 10 h −1 Mpc), however, the spatial distribution of galaxies is affected not only by our cosmological model but also by all of the poorly understood physics of galaxy formation and evolution.
We can bypass our ignorance of galaxy physics via the use of so called "halo models," which are motivated by our understanding that galaxies form and reside in arXiv:2110.03701v1 [astro-ph.CO] 7 Oct 2021 gravitationally bound, virialized regions of dark matter known as halos (e.g., Neyman & Scott 1952;Peebles 1974;McClelland & Silk 1977;Scherrer & Bertschinger 1991;Jing et al. 1998;Benson et al. 2000;Ma & Fry 2000;Peacock & Smith 2000;Seljak 2000;Scoccimarro et al. 2001;Sheth et al. 2001;White et al. 2001;Cooray & Sheth 2002). These models assume that the clustering of galaxies can be fully described by (i) the clustering of their host halos and (ii) the way in which galaxies occupy these halos. The first piece can be accurately predicted by ΛCDM via cosmological N-body simulations of dark matter only (i.e., gravity only). The second piece can be described by a halo model, which provides a statistical description of how galaxies occupy halos (i.e., the galaxy-halo connection).
A key ingredient of the halo model is the Halo Occupation Distribution (HOD), which specifies via a few parameters the probability that a halo of mass M contains N galaxies (above some luminosity threshold) and provides a prescription for the distribution of galaxies within a halo (Berlind & Weinberg 2002;Berlind et al. 2003). The standard form of the HOD (Zheng et al. 2005) contains at most five free parameters that specify the mean occupation number of galaxies and assumes that galaxies trace the dark matter inside halos. Constraining these parameters when fitting to observational data provides a useful empirical measurement against which we can test competing theories of galaxy formation and evolution.
Several studies which use this procedure produce fits which, if taken at face value, would result in ruling out the ΛCDM + HOD model (e.g., Zehavi et al. 2011). The errors used in these studies are typically derived via the jackknife method, which has been shown to produce biased results (Norberg et al. 2009). Meanwhile, spectroscopic surveys such as the Sloan Digital Sky Survey (SDSS; York et al. 2000) are providing increasingly high precision measurements of galaxy clustering at a wide range of physical scales. If we wish to take advantage of the information present at small scales to constrain the galaxy-halo connection, it is imperative that we carefully understand and minimize the uncertainty in our modeling procedure.
Recently, Sinha et al. (2018) (S18, hereafter) developed a numerical mock-based modeling procedure that significantly improved the accuracy of HOD modeling. They compared the clustering of SDSS galaxies to a ΛCDM + HOD model, measuring the projected correlation function, group multiplicity function, and galaxy number density. Carefully controlling for systematic errors allowed them to interpret the goodness of fit of their model. They found that their best-fit HOD model was unable to jointly fit the clustering statistics, revealing significant tension between SDSS galaxies and their ΛCDM + HOD model. Because this tension did not exist when they considered only measurements of the projected correlation function (as is done in many studies), S18 demonstrated the value of adding additional statistics in small-scale clustering analyses.
Motivated by these findings, in this work we extend the procedure used in S18 in order to maximize the return from spectroscopic surveys. We do so with two main goals in mind: to tighten the constraints on our adopted five-parameter HOD model and to increase the tension found in S18. To accomplish these goals, we include in our analysis the same clustering statistics used in S18 (galaxy number density, projected correlation function, and group multiplicity function) as well as four additional clustering statistics: redshift-space correlation function, group velocity dispersion, mark correlation function, and counts-in-cells statistics. The fully numerical modeling procedure built by S18 allows us to accurately model all of these clustering statistics without having to rely solely on analytic predictions (which would limit us to two-point statistics) as is done in most studies.
This work can be seen as the second in a series of several works attempting to model the clustering of SDSS galaxies to constrain the galaxy-halo connection. Ultimately, our goal is to expand our HOD model to include features like assembly bias (Gao et al. 2005;Wechsler et al. 2006;Croton et al. 2007;Pujol & Gaztañaga 2014;Hearin et al. 2016;Pujol et al. 2017;Artale et al. 2018;Salcedo et al. 2018;Xu & Zheng 2018;Zehavi et al. 2018;Bose et al. 2019;Contreras et al. 2019;Padilla et al. 2019;Wang et al. 2019;Zentner et al. 2019;Salcedo et al. 2020;Hadzhiyska et al. 2020Hadzhiyska et al. , 2021bMontero-Dorta et al. 2021), spatial bias (Watson et al. 2012;Piscionere et al. 2015), and velocity bias (Guo et al. 2015a,b). However, in order to constrain these extra features, it is crucial that we first incorporate clustering statistics into our modeling procedure that are sensitive to these biases. Thus, in this work, as we add clustering statistics to heighten the tension found in S18, we do so keeping in mind the ability of each statistic to constrain extended features of the HOD model.
The projected correlation function has been utilized in many galaxy clustering studies (e.g., Zehavi et al. 2011;Leauthaud et al. 2012;Coupon et al. 2015) and could potentially be used to constrain assembly bias in future work (Zentner et al. 2014Vakili & Hahn 2019). A few studies have used the redshiftspace correlation function (e.g., Tinker et al. 2006b;Parejko et al. 2013;Guo et al. 2015b;Beltz-Mohrmann et al. 2020), which is sensitive to the velocity assumptions in our HOD model and thus could be used to investigate velocity bias. Even fewer studies have incorporated other clustering statistics, like the mark correlation function, which is sensitive to assembly bias (e.g., Zu & Mandelbaum 2018); the group multiplicity function, which is sensitive to both high-mass occupation (e.g., Zheng & Weinberg 2007;Sinha et al. 2018) and spatial bias (e.g., Beltz-Mohrmann et al. 2020); the group velocity dispersion, which must be sensitive to velocity bias; or counts-in-cells statistics, which are sensitive to assembly bias (e.g., Tinker et al. 2006a;McCullagh et al. 2017;Walsh & Tinker 2019;Wang et al. 2019;Beltz-Mohrmann et al. 2020).
Although the HOD model that we use in this work does not include exotic features, these new clustering statistics can still provide significant constraining power for our current model. In this work we use an optimal combination of these statistics on various scales to extract more clustering information from the SDSS survey than in all previous studies. Additionally, we greatly improve the accuracy of our modeling procedure compared to S18 and perform an extensive analysis of both the statistical and systematic errors present in our procedure, further pushing galaxy clustering analysis into the accurate regime. These improvements will allow us to tighten our constraints on the model and more robustly characterize the tension found in S18.
We discuss our observational data in Section 2, our clustering statistics in Section 3, our simulations and mock galaxy catalogs in Section 4, and our modeling procedure in Section 5. In Section 6 we discuss our selection process for the measurements with the most constraining power. We discuss our results in Section 7 and provide conclusions in Section 8. Finally, in the appendices we provide a detailed analysis of the accuracy of several components of our modeling procedure, including our fiber collision correction (Appendix A), our model estimation (Appendix B) and covariance matrix (Appendix C), and the ability of our procedure to recover the correct HOD parameters of mock galaxy catalogs (Appendix D).

OBSERVATIONAL DATASET
In this work, we make use of the seventh data release (DR7; Abazajian et al. 2009) of the Sloan Digital Sky Survey (SDSS; York et al. 2000). Specifically, we utilize the large scale structure samples from the NYU Value Added Galaxy Catalog (NYU-VAGC; Blanton et al. 2005). As in S18, we use a parent sample of over 500,000 galaxies covering the northern footprint of SDSS (see S18 for more details). The absolute magni- Note-The columns list (from left to right): the absolute magnitude threshold of each sample at z = 0.1; the minimum, maximum, and median redshifts, respectively; the effective volume of each sample (corrected for incompleteness); and the galaxy number density of each sample.
tudes of these galaxies have been k-corrected to restframe magnitudes at redshift z = 0.1 (Blanton et al. 2003) but, unlike in S18, have not been corrected for passive luminosity evolution. We find that, when comparing volume-limited samples constructed with and without the passive luminosity correction, the number density of the sample without the correction is closer to being constant in redshift. Therefore, we do not include this correction. From this parent sample, we construct two volumelimited subsamples, each complete down to a specified rband absolute magnitude threshold. To construct these samples, we first choose redshift limits, z min and z max . For both samples, we adopt a z min of 0.02 so that the redshift-space distortions are small compared to cosmological redshifts. We choose a value of z max such that the redshifted spectra of all galaxies with our limiting absolute magnitude would still fall within the redshift survey's apparent magnitude and surface brightness limits if placed at z max . Our low-luminosity sample is complete down to an r-band absolute magnitude of −19 in the redshift range 0.02 − 0.07, while our high-luminosity sample is complete down to an r-band absolute magnitude of −21 in the redshift range 0.02 − 0.158. We refer to our two samples as "−19" and "−21" henceforth. The redshift limits, median redshift, effective volume, and number density of these two samples are listed in Table 1. These values are slightly different than those used in S18.
Throughout this paper, we adopt a flat ΛCDM cosmological model with Ω m = 0.302 and h = 1. The comoving distances of SDSS galaxies are determined based on this assumption and have units of h −1 Mpc. Similarly, reported absolute magnitudes are actually M r + 5logh 1 .

CLUSTERING STATISTICS
Many studies which use clustering statistics 2 to constrain the galaxy-halo connection (e.g., Zehavi et al. 2011) focus on analytic predictions of these measurements. Such studies are limited to employing only those clustering statistics for which analytic predictions are possible. Given the difficulty in writing down the analytic form of an arbitrary clustering statistic, as well as the need to consider the impacts of survey geometry and incompleteness, this limitation unnecessarily excludes many clustering statistics with potentially useful information for the galaxy-halo connection. One of the key advantages of employing a numerical modeling procedure is that clustering statistics can be measured in identical ways on both the SDSS galaxies and on mock galaxy catalogs generated according to our model. Such a procedure thus allows us to make apples-toapples comparisons between the data and our theoretical model, allowing us to employ any arbitrary clustering statistic with the knowledge that any systematic errors (e.g., edge effects) are reflected the same way in both the model and the data. Berlind & Weinberg (2002), and more recently S18 and Hadzhiyska et al. (2021a), demonstrate that different clustering statistics are sensitive to different combinations of HOD parameters, and that analyses involving several different galaxy clustering statistics have the most power to constrain galaxy bias within an HOD framework. In their analysis, S18 employ the galaxy number density, the projected correlation function, and the group multiplicity function. With the aim of increasing the precision of the S18 results, we utilize four additional clustering statistics in this work: the redshiftspace correlation function, the average group velocity dispersion, the mark correlation function, and counts in cells statistics. We choose these new clustering statistics for their potential ability to constrain models beyond the standard HOD (i.e., those that include assembly bias, spatial bias, and velocity bias), which we intend to explore in future work. All of our clustering statistics (other than galaxy number density) are described below.
3.1. The projected correlation function: w p (r p ) In general, a correlation function is the excess number of galaxy pairs above that which is expected for a random distribution of points, as a function of pair separation. Widely used in studies of galaxy clustering (e.g., Zehavi et al. 2011, S18), the projected correlation function aims to remove the effect of redshift-space distortions by first counting pairs of galaxies in bins of their line-of-sight and projected components, π and r p , and then integrating over π: (1) We use the Zehavi et al. (2002) definitions of r p and π, and we count pairs of galaxies in 14 evenly spaced logarithmic bins of projected separation r p between 0.1 and 20 h −1 Mpc. We then integrate out to π max of 40 and 80 h −1 Mpc for the −19 and −21 samples, respectively. We calculate ξ(r p , π) with the Landy-Szalay estimator (DD − 2DR + RR)/RR (Landy & Szalay 1993), where DD, DR, and RR are the normalized numbers of datadata, data-random, and random-random pairs, respectively. We use the publicly available corrfunc (Sinha & Garrison 2019;Sinha & Garrison 2020) to compute our projected correlation function.

The redshift-space correlation function: ξ(s)
The redshift-space correlation function ξ(s) is the excess number of galaxy pairs above that which is expected for a random distribution of points, as a function of redshift-space (i.e., not projected) pair separation s. Because the redshift-space distortions of galaxies depend on their velocities, measuring ξ(s) will allow us to investigate the assumption that galaxies trace the velocity distribution of dark matter within their halo. To compute ξ(s), we once again use corrfunc and the Landy-Szalay estimator, counting pairs in the same bins of separation as those used for w p (r p ) (but not projected).
3.3. The group multiplicity function: n(N ) The group multiplicity function is the number density of galaxy groups as a function of the number of galaxies ("richness") in the group, n(N ) (e.g., Berlind & Weinberg 2002). We use the Berlind et al. (2006) friends-of-friends algorithm to identify groups: galaxies are linked together if their projected and line-of-sight separations are both less than a corresponding linking length. We adopt the Berlind et al. (2006) linking lengths of b ⊥ = 0.14 and b = 0.75, which are given in units of the mean inter-galaxy separation n −1/3 g , where n g is the sample number density. For the −21 sample, we measure groups in the following five (inclusive) bins of N : (5 − 6), (7 − 9), (10 − 13), (14 − 19), (20 − 27). For the −19 sample, we measure groups in the following seven (inclusive) bins of N : (5 − 6), (7 − 9), (10 − 13), (14 − 19), (20 − 32), (33 − 52), (53 − 84). These are the same bins as those used in S18, with the exception of the largest bin in S18 for each sample, which is excluded from this analysis.

Average group velocity dispersion: σ v (N )
We also measure the average velocity dispersion of galaxy groups as a function of richness, σ v (N ). The velocity of each galaxy in the group is first computed by subtracting the velocity of the entire group (we only consider line-of-sight peculiar velocities). For each galaxy group, we then calculate the group velocity dispersion. Lastly we find the average group velocity dispersion of all groups within a particular richness bin (e.g., groups containing 10 − 13 galaxies). We use the same richness binning as in the group multiplicity function. Like ξ(s), measuring σ v (N ) also probes the assumption that galaxies trace the velocities of the underlying dark matter particles within the halo.
3.5. The mark correlation function: mcf(s) The mark correlation function, mcf(s), is given by where DD(s) is the number of galaxy-galaxy pairs as a function of redshift-space separation s, and MM(s) is the sum of a weighted pair-product of the same galaxygalaxy pairs. Both DD and MM are normalized by, respectively, the total number of galaxy-galaxy pairs and the sum of the product of weighted galaxy-galaxy pairs, at all separations.
To compute mcf(s), we must first assign a weight (or "mark") to each galaxy. Motivated by Salcedo et al. (2018), who found that secondary bias of halos can be explained by a halo's distance to a massive neighbor, we choose as our mark the distance from a galaxy to the nearest galaxy cluster. If the galaxy is located in a cluster, then we take the distance to the cluster's center. To choose the minimum number of galaxies to be considered a "cluster," we take the best-fit HOD parameters from S18 and determine the average integer number of galaxies that would be placed in a halo of mass 10 14 M according to Equations 3 and 4 (see Section 4.2). With this definition, to be considered a cluster, a galaxy-group must contain 15 and three galaxies for the −19 and −21 samples, respectively. For our purposes, whether or not this definition accurately weights each galaxy by its distance to the nearest cluster is irrelevant. Ultimately, we only care whether or not this statistic contains information that can be used to constrain our model, a question we explore in Section 6. We again make use of corrfunc to compute mcf(s).

Counts in cells: P N (R)
Counts-in-cells statistics measure the probability of finding a given number of galaxies within a randomly placed finite region (e.g., a sphere) as a function of region size (e.g., radius). One special case of this is the void probability function (VPF), which measures the probability of finding no galaxies in a random region of space. Another variation of counts-in-cells is the "singular probability function," (SPF) or the probability of finding exactly one galaxy in a randomly placed region. Tinker et al. (2006a) and more recently McCullagh et al. (2017) and Wang et al. (2019) found that, in contrast to the projected correlation function, the VPF is sensitive to environmental variations of the HOD, and thus could be used to rule out certain HOD models. We measure counts-in-cells (both the VPF and the SPF) in spheres of evenly spaced bins with radii of 2, 4, 6, 8, and 10 h −1 Mpc. We again make use of corrfunc to compute VPF and SPF.

Simulations and halo catalogs
In this work we make use of dark matter only (DMO) cosmological N-body simulations from the Large Suite of Dark Matter Simulations project (LasDamas; McBride et al. 2009). These simulations were run on the Texas Advanced Computing Center's Stampede supercomputer using the public code gadget-2 (Springel 2005). We generate power spectra and initial conditions with cmbfast (Seljak & Zaldarriaga 1996;Zaldarriaga et al. 1998;Zaldarriaga & Seljak 2000) and 2lptic (Scoccimarro 1998; Crocce et al. 2006), respectively. All simulations were run with the following cosmological parameters (Planck Collaboration et al. 2014): Ω m = 0.302, Ω Λ = 0.698, Ω b = 0.048, h = 0.681, σ 8 = 0.828, and n s = 0.96. In particular, we run two sets of simulations: one to mimic the SDSS −19 sample (which we call Consuelo), and another to mimic the SDSS −21 sample (which we call Carmen). Starting from z = 99, we run Consuelo to z = 0.054, and we run Carmen to z = 0.132, which are roughly the median redshifts of the SDSS −19 and −21 samples, respectively 3 .
For each sample, we run two high-resolution simulations and 100 low-resolution simulations. We use the high-resolution simulations to estimate the model observables and the low-resolution simulations to construct 3 These values are the median redshifts of the SDSS samples in S18.
Our median redshifts (Table 1) differ slightly due to changes in how we process SDSS (see Section 2). Note-The columns list (from left to right): what each simulation is used for, the absolute magnitude threshold of the corresponding SDSS sample, the name of the simulation, the seeds used, the boxsize, number of particles, mass resolution, force softening, and the number of simulations.
a covariance matrix representing cosmic variance (see Section 5). Each of the low-resolution simulations differs in the phases of the density modes of the power spectrum, which is controlled by a seed supplied to 2lptic. The seeds used for the high-resolution simulations were chosen from the 100 low-resolution seeds to minimize the cosmic variance error in our model (see S18 for more details). The details of these simulations are given in Table 2.
We identify halos using a spherical over-density (SO; Lacey & Cole 1994) algorithm using the rockstar phase-space temporal halo finder (Behroozi et al. 2013). We find halos using the M vir halo definition, where the halo density depends on both cosmology and redshift, as given by Bryan & Norman (1998).

Building mock galaxy catalogs
We can construct mock galaxy catalogs by directly populating the dark matter halos from our simulations. This population is performed via the Halo Occupation Distribution (HOD). Specifically, the form of the HOD we use is the "vanilla" model used by S18 and previously . In this form of the HOD, the population statistics of dark matter halos depend only on halo mass, with central and satellite galaxies treated separately (Kravtsov et al. 2004;Zheng et al. 2005).
For a halo of mass M , the average number of central galaxies is given by Within the error function, there are only two parameters, M min and σ logM , which control the population statistics for central galaxies. M min is the mass at which approximately half of the halos host a central galaxy, while σ logM dictates how rapidly the central population goes to zero with decreasing mass. When we consider whether or not to assign a central galaxy to a specific halo of mass M , we draw a random number R from a uniform distribution on the interval [0, 1). If R < N cen , then a galaxy is assigned. This central galaxy is placed at the halo center and given the mean velocity of the halo.
If we have placed a central galaxy in a halo of mass M , then the mean number of satellite galaxies is given by The exact number of satellite galaxies we place in a specific halo of mass M is determined by drawing from a Poisson distribution with a mean N sat . While there is a dependence on N cen (and thus on M min and σ logM ), the satellite population statistics are primarily governed by three parameters: M 1 , α, and M 0 . More intuitively, halos of mass M 1 on average host one satellite galaxy, while α dictates how rapidly halo occupation increases with increasing halo mass. The parameter M 0 technically sets the halo mass below which we do not find any satellite galaxies, but in practice this parameter has not been well-constrained in studies which use luminosity samples similar to our own (e.g., Guo et al. 2015b, S18). After we choose the exact number of satellite galaxies in a specific halo, each satellite is given the position and velocity of a randomly chosen dark matter particle in the halo. Once we populate our dark matter halos with galaxies, we must build realistic mock galaxy catalogs that resemble our SDSS samples of interest. To do this, we first need to transpose the mock galaxies from Cartesian to spherical coordinates by positioning an observer at the center of the box and converting the positions of the galaxies into RA, DEC, and co-moving distances. Due to the smaller volume of the SDSS sample, we are able to carve out four independent mock galaxy catalogs from each simulation box. We then need to incorporate the same systematic effects that plague our observational dataset, such as redshift-space distortions, sample geometry, and incompleteness.
To introduce redshift-space distortions, we determine the line-of-sight peculiar velocities of the galaxies and calculate the redshift as 1 + z = (1 + z cosmological )(1 + z doppler ), where z cosmological is the cosmological redshift and z doppler is the redshift due to the radial peculiar velocity. We then eliminate any mock galaxies outside the redshift limits or sky footprint of our SDSS sample. This procedure ensures that any effect on measured clustering statistics due to sample geometry are present in both the mock galaxies used in our model and in the SDSS data.

Correcting for fiber collisions
In SDSS, spectroscopic fibers cannot be placed closer to each other than their own thickness (55 ). This limitation results in approximately 7% of targeted galaxies lacking a spectroscopically measured redshift due to their proximity to a neighboring galaxy. Modeling these fiber collisions is difficult primarily because it requires a simulation with a large enough volume and a high enough resolution to reproduce the flux-limited footprint of SDSS. Instead of incorporating fiber collisions into our model, we choose to make a correction to the SDSS data. It should be noted that this is the one case where we necessarily diverge from our general modeling philosophy -instead of incorporating this systematic error into our generative model, we attempt to remove it from the SDSS data.
To accomplish this task, we first apply the nearest neighbor correction to galaxies lacking a spectroscopic redshift (Zehavi et al. 2002). We then estimate the error in measured clustering statistics from using only the nearest neighbor correction and adjust our measurements accordingly. Our full procedure is detailed in Appendix A, but we provide a brief summary here. We utilize the SDSS plate overlap regions to estimate the error in our clustering statistics that results from applying this nearest neighbor correction. In overlap regions, we have spectroscopic redshifts for many galaxies that are within 55 of a neighboring galaxy. These regions thus allow us to investigate the impact of applying the nearest neighbor correction. Briefly, we find that the nearest neighbor correction alone is not sufficient for many of the clustering statistics we use. Thus, we choose to adjust the values of our observables to account for this error. We use the adjusted observables in our analysis of SDSS (Section 7). We provide both the original nearest neighbor SDSS observables and the adjusted observables in the machine-readable format.

SUMMARY OF MODELING PROCEDURE
Our main goal in this work is to constrain the galaxyhalo connection. To achieve this end we implement a fully numerical modeling methodology, based on and adapted from S18. We summarize the procedure here, as well as highlight a few key differences from S18.
To employ our numerical modeling procedure, we utilize the Texas Advanced Computing Center's Stam-pede2 supercomputer. We use a Markov Chain Monte Carlo (MCMC) algorithm to explore the HOD parameter space.
In particular, we developed a Cimplementation of the popular affine-invariant sampler emcee (Foreman-Mackey et al. 2013), which we call emcee in c 4 .
An MCMC (or "chain") involves probabilistic sampling of an unknown posterior distribution of parameter values given some assumed prior distribution. We employ flat prior distributions on the same parameter ranges given in S18. At each HOD point we sample, we must evaluate the likelihood that this HOD model could have generated a dataset with the same clustering as SDSS. Given a K-dimensional vector D of observables measured on the SDSS dataset and a corresponding vector M of the same observables for an HOD model ("model observables" hereafter), this likelihood is given by where C is the K-dimensional covariance matrix of these observables representing cosmic variance (see Equation  6). Ignoring the factor of −1/2, the term within the exponential is χ 2 .
In the context of a numerical modeling methodology, ideally we would obtain the observable vector M by first populating a large number of high-resolution DMO simulations with galaxies according to our set of HOD parameters. From these populated simulation boxes, we would then carve out SDSS-like mock galaxy catalogs, measure all observables, and find the average value of each observable across all mocks. In practice, we need to populate halos in a large enough volume in order to yield a cosmic variance error in M that is sufficiently smaller than the uncertainty in D so as not to dominate the overall error budget. As discussed in Section 4.1, S18 found that two boxes were sufficient to accomplish this task for an analysis using n gal , w p (r p ), and n(N ). We continue forward with this approach using the same boxes as S18 and investigate the cosmic variance error in M for the expanded set of clustering statistics we use. emcee in c was made to meet the software restrictions of the Stampede2 supercomputer. We have verified that emcee in c performs identically to emcee.  For each HOD point we explore in the MCMC, we directly populate the halos in these two boxes 5 . We have a choice as to whether to measure each clustering statistic either on the full box(es) or on mock galaxy catalogs carved out from these boxes (see Section 4.2). The goal is to measure each statistic in a way that yields model observables closest to the ideal case described above (i.e., mean of many mocks). Briefly, utilizing the two low-resolution counterparts to our high-resolution boxes, we compare different methods of measurement and choose the method for each clustering statistic that produces model observables closest to the mean measurement across 400 low-resolution mocks. We describe our full procedure in Appendix B. We summarize the results of this test in Table 3, where, for each sample, we indicate whether we choose to measure each clustering statistic on a box (B) or on mocks (M) and which box(es)/mock(s) are used (see table description for more details).
We take one last step when estimating the model observables in our chains. For a given HOD, the process of populating dark matter halos with galaxies is stochastic (see Section 4.2) and is controlled with a "population seed." For a fixed HOD, varying this population seed primarily affects the number of centrals/satellites assigned to each halo and can lead to significant differences in 5 To meet the memory requirements of the Stampede2 compute nodes, we use downsampled versions of these boxes in which we only keep a randomly selected 5% of the dark matter particles in each halo. This still leaves ample particles on which to place satellites (∼ 500 times more particles than satellites on average). Thus, downsampling results in no loss of accuracy. measured clustering statistics (see Appendix B and Figure 14). Even when using a fixed population seed, a slight change in HOD can also lead to significant differences in measured clustering statistics. Such differences do not arise from an actual difference in the likelihoods of these two different models but rather from the noise in how we are estimating the model observables. To reduce this noise, for each HOD point in our chain, we populate our simulation boxes four times, using four fixed population seeds. We then measure the observables (in the way given in Table 3) on each of these four resultant galaxy catalogs. We use the average of these four measurements as our model observables. Similar to the model observables M, ideally we would obtain the covariance matrix C using a large number of high-resolution mock galaxy catalogs. Because of the high computational cost associated with running many high-resolution DMO simulations, we instead construct the covariance matrix from 400 mock galaxy catalogs built from 100 low-resolution boxes, all populated according to a set of fiducial HOD parameters (see Sections 6 and 7). The elements of the covariance matrix are given by Here, the sum is taken over the N = 400 mocks. The values y i and y j are the i th and j th observables measured on each mock and thus vary in the sum. The values y i and y j are the mean values of the i th and j th observables, respectively, and thus are fixed in the sum. Each diagonal element, C ii , of the matrix is simply the variance among 400 mocks for observable i. In this paper, we refer to √ C ii as the "cosmic variance uncertainty" of observable i. For simplicity, we denote the cosmic variance uncertainty of an arbitrary observable as σ obs .
There are three important differences between the covariance matrices we use and those used in S18. First, we use twice the number of mocks as S18, resulting in less noise in our covariance matrix. Second, the lowresolution boxes we use to construct the matrix all have the same cosmology and halo definition as the highresolution boxes we use to obtain our model observables. This is not the case in many of the chains run in S18. Third, instead of using the raw covariance matrix produced by these 400 mocks, S18 rescale the elements of the covariance matrix by the ratio of SDSS to mock measurements (see S18 for more details). When evaluating L(D|M ), however, the covariance matrix represents the variation in clustering statistics calculated on different mock realizations generated by the same set of HOD parameters. Therefore, we instead simply use the raw covariance matrix in our likelihood calculation and do not rescale its elements. While we view all three of these points as improvements to the modeling procedure, we do not quantify their impact.
In addition to these changes, we investigate sources of noise in our covariance matrix in Appendix C. We examine only the noise in the cosmic variance uncertainty of each observable, ignoring off-diagonal elements of the matrix. To summarize, for the −19 sample, the major source of noise is due to the number of mocks we use, which gives an error of ∼3.5% on all values of σ obs . For a few observables, the largest source of noise is due to our choice to use a fiducial matrix instead of re-building the matrix for each set of HOD parameters we evaluate in the chain. Still, this error is at worst only 10%. For the −21 sample, the largest source of noise comes from using low-resolution mocks instead of high-resolution mocks to build our matrices. In the worst case, we find that we overestimate the errors by ∼15% for several small scales of w p (r p ) and ξ(s). We do not make any corrections to our matrix to account for this discrepancy, believing that overestimating our errors will lead to more conservative posterior results.
While the tests in Appendix C quantify the error in the elements of the covariance matrix, it is a much more difficult task to quantify the joint impact of these errors. In particular, we wish to know how these errors affect our constraints on SDSS. In general, using more observables will produce tighter constraints, but this improvement occurs at the cost of more noise in the covariance matrix and thus in the resulting constraints. We therefore seek to obtain a set of observables that provide both tight and reliable constrains, a task we address in the following section.
6. OPTIMIZATION ALGORITHM From our set of N observables, we seek to find the subset of size K ≤ N that gives the tightest constraints on our HOD parameters, while still being robust with respect to the systematic errors in our modeling procedure. Because we do not know what an acceptable value for K is, exploring all combinations quickly becomes an impossible task. With N = 63, as in our −21 sample, N choose K for K : 1 → N represents ∼ 10 19 unique combinations! Therefore, we design and implement an algorithm that is an approximation of this task.
The algorithm we implement orders our N observables by joint constraining power. When choosing the K th observable, we pick the one that, when combined with our K −1 previously chosen observables, produces the tightest constraints on our HOD parameters. Once the order is created, we analyze the reliability of our constraints for different values of K. This analysis leads to a choice of K observables which we utilize in an MCMC to constrain the HOD of our SDSS volume-limited samples. In the following subsections, we describe our algorithm, its raw results, and our choice of observables.

Algorithm design
The first step in our algorithm is to set up a grid of HOD parameters which we will use to explore the constraining power of different combinations of observables. Because some of our parameters are strongly correlated (e.g., logM min and σ logM ), a uniformly-spaced grid contains a large number of unrealistic HOD parameter combinations. Therefore, to set up a grid we perform the following steps: 1. Choose a fiducial set of HOD parameters. 2. Using these parameters, create a high-resolution fiducial mock galaxy catalog. 3. Measure n gal and w p (r p ∼ 0.3 h −1 Mpc) on this mock. 4. Run an MCMC with this mock's observables as the "data." For each HOD point: • Calculate and record all N model observables. (See Section 5 for a description of how model observables are calculated.) • Evaluate likelihood L init using only n gal and w p (r p ∼ 0.3 h −1 Mpc).
This procedure produces an initial non-uniform grid of fairly reasonable HOD points which could have generated our fiducial mock. We choose to perform this analysis on a mock galaxy catalog for two reasons. First, we wish to know whether a chain run with a specific set of observables is able to recover the true HOD, which is not something we can determine with SDSS. Second, we wish to choose observables based on their universal constraining power. By using a mock galaxy catalog, we avoid choosing observables that are over-fit to SDSS. We are, however, potentially over-fitting to this particular mock galaxy catalog. Therefore, we perform steps 2-4 for four different mock galaxy catalogs. These mocks have the same HOD and only differ due to cosmic variance. Thus, we construct four different grids which we use to explore the constraining power of different combinations of observables.
For a given mock/grid, we quantify the constraining power of n gal and w p (r p ∼ 0.3 h −1 Mpc) by measuring the standard deviation of HOD values of each parameter. We call this initial constraint σ par,init . We wish to examine how different combinations of observables can improve upon the initial constraints. To do so, we utilize the fact that the value of L init at a particular point in HOD-space is proportional to the density of points in the grid at that location. Because we measure all N observables when constructing the grid, for each HOD point we can compute a likelihood L K using any arbitrary set of K observables and the appropriate K by K covariance matrix. We can then assign each point in the grid a weight L K /L init . The weight assigned to each HOD point is proportional to the density of points we would have at this location if we had we run a chain on this mock using this set of K observables. Computing a weighted standard deviation σ par,K for each parameter gives an estimate of the constraints from said chain. This process, known as "importance sampling," thus provides a way of estimating the constraining power of an arbitrary combination of K observables.
We use importance sampling to order our N observables by their cumulative constraining power. In short, we choose observables, one at a time, always choosing the observable that, when combined with the previously chosen observables, best constrains all parameters of interest. To choose the K th observable, we perform the following steps: observables not yet chosen, pick a "trial observable." 2. In one of the four grids, using this trial observable and the K − 1 previously chosen observables, compute a likelihood L K for each point. 3. Weight each point in the grid by L K /L init . 4. Compute a weighted standard deviation, σ par,K , for each HOD parameter. 5. Add in quadrature the fractional reduction in the constraints on each parameter achieved by including this observable: M grid = par σ par,K σ par,K−1 2 . 6. Repeat steps 2 -5 for each of the grids. Add the values of M grid to get the metric M for this trial observable. 7. Repeat steps 1 -6 for each of the N − (K − 1) possible trial observables. 8. Choose the trial observable with the lowest value of M .
Once the K th observable is chosen, we move on to choose the (K + 1) th until all N observables have been ordered. The metric M always compares the improvement in constraints when going from K −1 to K observables. This metric rewards observables which constrain parameters that have not already been well-constrained by the first K −1 observables. Additionally, by minimizing the sum of M grid across multiple grids, the choices we make are more robust to the peculiarities of a specific mock due to cosmic variance.
In principle, this procedure could be used to order all observables from K = 1 to K = N . In practice, however, we make the decision to start with n gal and w p (r p ∼ 0.3 h −1 Mpc) and only choose observables from K = 3 to K = N . In general, importance sampling works well when the target distribution is close to or contained within the initial distribution. This criteria is met when choosing the third observable but is not always true when we attempt to choose the first or second. Our choice to build the initial grid using n gal and w p (r p ∼ 0.3 h −1 Mpc) thus affects all successive choices of observables. We assume that n gal will be included in any clustering analysis, given the low computational cost and high information content. Thus, including it in the construction of this grid seems reasonable. We first attempted to create our grid using only n gal , but this grid contained HOD points very unlike SDSS, particularly in the satellite parameters. In preliminary analysis, we found that w p (r p ∼ 0.3 h −1 Mpc) is particularly good at constraining the satellite parameters. Thus, we include it to acquire a more reasonable starting grid of HODs.
Importance sampling works if the grid contains a sufficient number of points in the region of the target distribution. As we go to higher values of K, the target distributions occupy a smaller region in HOD space, and the density of our grid can become an issue. Ideally, each time we choose a new observable, we would construct a new grid by running a new chain with the already chosen observables. This denser grid would then be used to choose the next observable. Such a procedure would minimize the error in importance sampling but has a high computational cost. We make the decision to build a denser grid whenever the importance sampling procedure produces noisy posterior distributions, which we determine by visual inspection. We note at which number of observables we make the decision to construct denser grids in the following section.

Ordering of observables
We apply the algorithm outlined in the previous section to order the observables according to their potential constraining power. We perform this procedure (separately, for each sample) on mocks constructed using HODs appropriate for the −19 and −21 SDSS volumelimited samples. The fiducial HOD parameters we use in building the grids are given in Table 4. These parameters are the best-fit HOD values from a chain run on SDSS 6 using only n gal , w p (r p ), and n(N ) and thus con- Figure 1. Projected constraints (1-σ) on four of our HOD parameters as we increase the number of observables. The constraints for the −19 and −21 mocks are shown in blue and red, respectively. The order of the observables is determined by our algorithm and is shown in Table 5. Each dashed line gives the projected constraints from one of our grids while the solid line shows the average of all four grids. For comparison, the constraints on our SDSS samples from a chain run using n gal , wp(rp), and n(N ) (as was done in S18) are shown with dotted horizontal lines. Note-Unless otherwise noted, we use these fiducial HOD parameters to construct the covariance matrices and mocks we use throughout this paper.
stitute reasonable parameter values for SDSS. For each sample, we construct both the mocks and the covariance matrix from these HODs. For the −21 mocks, the value of logM 0 is too low to have any impact on the mathematical form of our HOD. For the −19 mocks, however, logM 0 is large enough to have a significant impact on halo occupation. Therefore, we decide to simultaneously optimize using all five parameters for the −19 mocks but ignore logM 0 when optimizing the −21 mocks (i.e., in our calculation of M grid ).
In Figure 1, we present the results of our algorithm for four of the five HOD parameters (excluding logM 0 ).
The dashed/solid lines show, for our −19 (blue) and −21 (red) mocks, the estimated posterior constraints (weighted standard deviations) on each parameter as we increase the number of observables in the order chosen according to our algorithm. Each dashed curve gives the estimated constraint obtained from one of the grids (mocks), while the solid curve gives the average constraint. The horizontal dotted lines show the constraints achieved on our SDSS volume-limited sample when we only use n gal , w p (r p ), and n(N ) (NWG, hereafter), as was done in S18 (see Section 7).
In Figure 1, we can see that our average projected constraints are tighter than the NWG constraints for all parameters by the time we reach ∼ 10 observables. Thus, we need only half as many observables as are used in the NWG chain to obtain the same constraining power. For many of the parameters (e.g., −21 central parameters), we are outperforming the NWG constraints with even fewer (∼ 5) observables. These results highlight the strong constraining power that comes from combining various scales of multiple clustering statistics.
Comparing the two samples, we project tighter constraints on three of the parameters for the −21 mocks compared to the −19 mocks, with α as the only pa- rameter better constrained for −19. As can be seen by the horizontal dotted lines, this trend is also true when constraining SDSS using only n gal , w p (r p ), and n(N ). Additionally, the projected constraints are all quite similar among the four mocks for −21, while there is a greater dispersion among the −19 mocks, particularly for logM min and σ logM . Because the fiducial value of σ logM in the −19 mocks (0.22) is closer to zero than in the −21 mocks (0.46), the posterior contours are often cut off by our flat prior lower bound (σ logM cannot be negative). This results in the shapes of the contours having more variety between the four mock realizations, affecting both logM min and σ logM due to their correlation.
Remarkably, we are also able to obtain tight constraints on logM 0 for the −19 mocks using a collection of 10 observables. We show this result in Figure 2, which has the same format as Figure 1 but excludes the −21 results. When using only n gal , w p (r p ), and n(N ), logM 0 is almost entirely unconstrained. With the availability of the clustering statistics we introduce in this paper, our algorithm very quickly selects observables which constrain logM 0 .
At the end of Section 6.1, we describe how we construct denser grids whenever we deem it necessary. For both samples, we choose observables 3 − 5 using our K = 2 grid (constructed with just n gal and w p (r p ∼ 0.3 h −1 Mpc)). After choosing observable 5, the posterior contours arising from importance sampling appear to be noisy. We thus reconstruct denser grids for both samples using the first five observables (K = 5 grids). For −21, we use these K = 5 grids to choose observables 6 − 63. For −19, constraining logM 0 notably has a significant impact on the allowable region of logM min and logM 1 , which are both anti-correlated with logM 0 . Once we begin to obtain tight constraints on logM 0 , the den- Note-The type of clustering statistic and the bin number (1-indexing) for the first 15 observables chosen (in order) for each sample. Note that the first two observables were fixed. We include the full order in the machinereadable format.
sity of our grid is again no longer sufficient. Therefore, for −19, we only use the K = 5 grid to choose observables 6 − 8 and switch to an even denser grid at K = 8. We then use this K = 8 grid to choose observables for all larger values of K.
In Table 5, we show the order in which the first 15 observables are chosen for each sample. We indicate the clustering statistic and the bin number, where Bin 1 is the smallest bin of the corresponding clustering statistic. We include the full order in the machine-readable format. Several small and intermediate ( 1 h −1 Mpc) scales of w p (r p ) and ξ(s) are chosen early for both samples. While these scales are highly correlated and thus provide similar information, their joint constraining power is apparently strong enough to justify their inclusion. For −19, among the first 15 observables are four scales of n(N ), while only one scale is found among the first 15 choices for −21. This statistic is sensitive to the satellite population and therefore may be more relevant for −19, which has a higher satellite fraction. On the other hand, for −21, we find four scales of mcf(s) among the first 15 choices, while finding only one scale for −19.
While commenting on the precise reasons for the choice of each observable is beyond the scope of this paper, we wish to get a general picture of the ability of each individual clustering statistic to constrain each HOD parameter. To estimate the results of running a chain with only one clustering statistic (e.g., mcf(s)), we perform a task very similar to the one described in Section 6.1. Using the sparsest (K = 2) grids, we weight every point in our grid by L stat /L init , where L stat is the likelihood calculated using all scales of one clustering statistic and L init is the likelihood calculated using only n gal and w p (r p ∼ 0.3 h −1 Mpc). We then importance sample the grid to estimate the constraints on each HOD parameter that would be obtained from running a chain with this clustering statistic.
We show the results of this test in Figure 3. Like in Figure 1, each panel shows the estimated constraints on a different HOD parameter, with constraints for −19 and −21 shown in blue and red, respectively. The height of each smaller, filled bar shows the constraints on an individual mock when using the clustering statistic 7 indicated on the horizontal axis, while the empty vertical bar shows the average constraint across all four mocks. We choose to treat VPF(R) and SPF(R) as one statistic (P N (R)) in this figure.
For the central parameters, logM min and σ logM , it is notable that no one clustering statistic does a particularly good job of constraining the −19 mocks. Contrast this with the −21 mocks, where ξ(s) and mcf(s) (as well as w p (r p ) and P N (R) to a lesser extent) each alone provide very tight constraints. In Figure 1, however, we see that the achievable constraints on logM min for −19 and −21 are quite similar. This highlights the advan-tage of using a variety of clustering statistics measured at different scales.
For both samples, we see that w p (r p ) and ξ(s) provide excellent constraints on the satellite parameters, logM 1 and α. We can thus see why so many scales of these observables are chosen early for both −19 and −21. For logM 1 , it is notable that mcf(s) again provides extremely tight constraints for −21 but not for −19, further illustrating why so many mcf(s) scales are chosen early. Similarly for α, n(N ) provides better constraints for −19 than for −21, illustrating why more scales are chosen for that sample.
We should caution that the results presented in Figure 3 are estimates. As we have already mentioned, importance sampling works well if the grid has a sufficient density of points in the neighborhood of the target distribution. Because we are using our sparse K = 2 grids, this condition is not generally met for the results shown in Figure 3. This fact is reflected by the large variation from mock to mock for a given statistic (e.g., ξ(s) for −19). We only show this figure to provide some qualitative insight into the order of observables presented in Table 5.
Looking again at Figures 1 and 2, we can see that by ∼ 20 observables, the constraints on virtually all parameters begin to asymptote. There is still, however, a slight improvement in the constraints beyond 20 observables. From these figures alone, however, we do not have a handle on the noise present in each of these estimates of our constraints. It is thus difficult to gauge how much information we gain as we increase the number of observables. We address this point in the following section.

Choosing an optimal set
In Figure 1, we see that using more observables always improves the constraints on our HOD parameters. Using more observables also adds more elements to the covariance matrix. All elements of the covariance matrix contain some noise because we calculate them using a finite number of mocks. In principle this noise would disappear if we had an infinite number of mocks, but in practice we have only 400. Because calculating a likelihood involves an inversion of the covariance matrix (see Equation 5), any noise in the matrix is propagated into the likelihood. Using more observables introduces more noise into the calculation of a likelihood and thus into the posterior results of an MCMC. Therefore, as we increase the number of observables, there is a trade-off between the increase in constraining power and the added noise in these constraints.
We wish to know, given a particular set of K observables, how our constraints would change if we were to estimate the covariance matrix using a different set of 400 mocks. We design a test to estimate this "error" in our constraints for different values of K. We show the results of this test in Figure 4. As in previous figures, the −19 and −21 samples are plotted in blue and red, respectively. The dark solid lines show the grid-averaged estimated constraint (i.e., the same lines as are in Figure 1), while the shaded region shows our estimate of the error in these constraints. To estimate this error, we perform the following steps for each of our grids and for each ordered subset of K observables: 1. Treating the K by K covariance matrix as a multivariate Gaussian, randomly sample 400 points from this distribution. 2. From these 400 points, construct a new "resampled" covariance matrix. 3. Compute a new likelihood L K for each point in the grid 8 . In this calculation, use the resampled K-dimensional covariance matrix. 4. Importance sample the grid, weighting each point by L K /L init (see Section 6.1 for description of L init ). 5. Re-calculate the weighted posterior standard deviations, σ par,resamp , of each HOD parameter. 6. Repeat steps 1 − 5 100 times, yielding 100 values of σ par,resamp for each HOD parameter. 7. For each parameter, compute an error δ(σ par ) by taking the standard deviation of these 100 values of σ par,resamp .
Performing these steps for each of our four grids yields four estimates of δ(σ par ) for each HOD parameter as we increase K. For each value of K, we choose the largest value of δ(σ par ) among these four estimates as our error 9 .
As expected, we see in Figure 4 that δ(σ par ) increases as we increase the number of observables. We wish to choose a value of K where both the constraint (σ par ) and the error on the constraint (δ(σ par )) are small. To accomplish this task, we employ the "one standard-error rule": we choose the lowest value of K such that the constraint at this value is within one standard error of the constraint at all higher values of K. We have to simultaneously meet this condition for all parameters of interest to us, which pushes us to higher values of K than are necessary for some parameters. Coincidentally, 8 We use the most dense available grid for each value of K. 9 We also explore taking the average error among the four grids.
The choice makes very little difference so we opt for the more conservative approach (i.e., larger errors). for both samples, this happens to be the same number: K = 36. We mark this value of K in Figure 4 with a dot. We also indicate the estimated constraining power for K = 36 with a horizontal dashed line. Although not shown here, logM 0 is also taken into account when choosing K for the −19 sample.
In choosing a value of K, we have focused solely on the constraining power achievable by using a particular combination of observables. The other critical requirement of our set of observables is that, when used in an MCMC, the true HOD parameters used to construct our mock catalogs should lie somewhere inside of our posterior contours. Where exactly these parameters lie depends on the cosmic variance of the dataset on which a chain is run and thus will differ for each of our four mocks. Testing whether any bias exists for a given set of observables, therefore, is a far greater challenge. Instead, we choose to just visually inspect whether or not the true HOD parameters lie within the 95% posterior probability region for each of our HOD parameters. We find that for our K = 36 observables, this condition is met for all mocks in both samples. In Appendix D, we go one step beyond importance sampling and run a full MCMC on one mock from each sample, using our K = 36 observables. We again find that the true HOD parameters are recovered for both samples.
We highlight the variety of observables chosen by our algorithm in Figure 5. We show measurements of all scales of each clustering statistic (excluding n gal ) for our two volume-limited SDSS samples. The observables we ultimately choose are indicated with solid points. We also provide a shaded region to indicate the size of the cosmic variance uncertainty used in our analysis, which represents the variance associated with our process of generating SDSS-like mocks.
The degree of correlation between any two observables i and j is given by the correlation coefficient, where C ij is an element of the covariance matrix as defined in Equation 6. Figure 6 shows the correlation matrices for the −19 and −21 samples. On both axes the observables are placed in the order in which they are chosen according to our algorithm (listed in Table 5). We only show the 36 observables we use in our analysis, excluding all others. The shading of a particular cell represents the degree of correlation between an observ- Figure 5. We show here the SDSS measurements of the observables we use in our analysis, with each panel showing measurements for a different clustering statistic. For the −19 (blue) and −21 (red) samples, the solid lines give the SDSS measurements on all scales we consider using of each clustering statistic. The points mark the scales we ultimately choose to use in our analysis (see Section 6.3). We also use galaxy number density, which is not shown here. For illustrative purposes, we show the cosmic variance uncertainty from our 400 mocks with the shaded region centered on the SDSS measurements.
able on the x-axis and the corresponding observable on the y-axis. The most remarkable feature of Figure 6 is that both matrices are extremely diagonal, particularly when compared to those of S18. By using our algorithm, we generally avoid choosing highly correlated observables which contain redundant information. Looking deeper, the observables for the −19 sample exhibit a higher degree of correlation than the corresponding observables for the −21 sample. This is true in general and not just for the observables chosen here. Still, we can see that several highly correlated observables are chosen early for both samples, particularly for −19. While this is not what we anticipated, it is clear that these correlated observables contain enough joint constraining power to be selected by our algorithm.
In this section, we have established a framework for selecting optimal observables that can be used to constrain the clustering of SDSS galaxies for a given HOD model. While this set of observables is specifically optimized to the HOD model employed in this work, this procedure could be repeated if we make any changes to our HOD model (e.g., adding assembly bias) in the future. In the next section, we use our chosen observables to constrain and test our standard HOD model against our two SDSS volume-limited samples.

SDSS results
With a set of 36 optimal observables chosen to produce tight and reliable constraints, we run an MCMC to constrain the HOD for each of our two volume-limited SDSS samples. We refer to these chains as the "OPT" chains henceforth. We wish to compare the results of these OPT chains to those of S18. Because we have altered the SDSS dataset (Sections 2 and 4.3) and details of the modeling procedure (Section 5), a direct comparison is not appropriate. Therefore, we also run chains using only n gal , w p (r p ), and n(N ), as was done in S18, and refer to these chains as "NWG." For each sample, we use the same HOD to construct the covariance matrix for both the OPT and NWG chains. The parameters we use to construct this matrix are given in Table 4.
All chains were run on the Texas Advanced Computing Center's Stampede2 supercomputer using 272 KNL CPUs spread across four compute-nodes. Each chain was run using 542 walkers for ∼1000 iterations and thus involved ∼500,000 total HOD model evaluations. For  Table 4. Observables are placed in the order in which they are chosen by our algorithm (see Section 6.3). We show only the first 36 observables (i.e., those which we use in our analysis).
the OPT chains, each parameter evaluation took approximately five CPU-minutes, and thus a full chain took ∼45,000 CPU-hours (or ∼650 node-hours). We determine that a chain has converged when the probability distribution of each parameter stabilizes, which occurs after ∼200 iterations.
The joint parameter constraints resulting from these chains are shown in Figure 7. The left and right panels show the constraints on the −19 and −21 samples, respectively. The top panels show σ logM vs. logM min while the bottom panels show α vs. logM 1 . The blue contours show the constraints achieved in the NWG chains, while the red contours show our results when running the OPT chains. The shading shows the 68 and 95% probability regions.
It is clear from Figure 7 that we achieve significantly tighter constraints on the joint distributions of all HOD parameters in the OPT chains compared to the NWG chains. With these tighter constraints, we demonstrate the power of combining an assortment of clustering statistics measured at various physical scales. As a result, we can now detect a significant difference in the values of σ logM between the two samples, with the −21 sample having a higher value than the −19 sample. This result is to be expected, given that the value of σ logM can roughly be interpreted as the scatter in halo mass at fixed luminosity, which is generally greater for more luminous galaxies (Behroozi et al. 2010). However, several previous HOD modeling works that rely on clustering (e.g., Sinha et al. 2018;Zentner et al. 2019), have not successfully detected this distinction in σ logM between different luminosity samples. Moreover, the few studies that have found that σ logM generally increases with luminosity (Zehavi et al. 2011;Guo et al. 2015b) did not have tight enough constraints to detect this difference to the level of significance that we achieve in this work.
In the NWG chains (and the chains in S18), logM 0 is entirely unconstrained for both samples. Using our chosen set of observables, however, we are able to obtain tight constraints on logM 0 for the −19 sample. In Figure 8, we show the joint constraints on logM 0 and logM 1 for the −19 OPT chain. We exclude the results for the NWG chain, where the contours extend all the way down to the lower prior bound on logM 0 at a value of 6. This parameter is anti-correlated with both logM 1 and logM min (not shown). Our ability to constrain logM 0 , thus, affects the allowable values of these two parameters.
In Table 6 we present the marginalized constraints on each parameter for the NWG and OPT chains. For each sample, we give the median parameter values as well as the upper and lower limits corresponding to the 84 and 16 percentile parameter values. To summarize the level of improvement, we divide the range of the inner 68% of parameter values for the NWG chain by the corresponding range for the OPT chain. We report this result in Table 6 in the column labeled "Constraint Ratio." This ratio gives the factor by which we shrink the marginalized constraints when using our set of chosen observables, compared to using just the S18 observables. As was evident in Figure 7, we see improvement in the constraints on every parameter. For both samples, we can see a particularly strong improvement in the constraints on the central parameters, logM min and σ logM , and a weaker improvement for the satellite parameters, logM 1 and α. The greatest improvement is for logM 0 in the −19 sample, again reflecting the fact that this parameter is unconstrained when using only n gal , w p (r p ), and n(N ).
In addition to tightening the constraints on our HOD parameters, we have also significantly heightened the tension between our model and SDSS. In Table 7, for each chain, we provide the best-fit HOD parameters and their associated χ 2 , degrees of freedom, and p-values. For the NWG chains, the p-values of the best-fit points suggest that our model is sufficient to describe the clustering of SDSS 10 . For the OPT chains, however, the pvalues indicate that our model is unable to accurately match the clustering of SDSS. We find significant tension of 4.5 and 4.1 σ for the −19 and −21 samples, respectively.

1.688
Note-Marginalized constraints on SDSS for each chain shown in Figure 7. We present the median parameter values along with upper and lower limits corresponding to the 84 and 16 percentiles respectively. We also provide the ratio of constraints (inner 68 percentile range) of the NWG chain to the OPT chain. These numbers indicate the factor by which we have improved our constraints. Note-Best-fit HOD parameters from all chains shown in Figure 7. We also indicate the goodness of fit of each parameter combination with a χ 2 , the number degrees of freedom (d.o.f), and a p-value. sample, using the set of observables we choose in Section 6.3. Unlike the S18 results, we are able to obtain very tight constraints on logM0, which has a significant anti-correlation with logM1.
To explore this tension further, we wish to identify which observables are contributing the most to the high value of χ 2 for each sample. In Figure 9, for each observable, we show the value χ obs = (D −M )/σ obs , where D is the measurement on SDSS, M is the best-fit model measurement, and σ obs is the cosmic variance uncertainty. For each clustering statistic, indicated with a color-coded label, we show χ obs values for all observables. We indicate with larger, darker points those observables that are actually used in our analysis and thus contribute to the overall χ 2 value. Because many of the observables are correlated, the overall χ 2 value is not just the sum of squares of these values, but we can still make some general conclusions about which observables are contributing the most to our high χ 2 values. For the −19 sample, notably both n gal and small scales ( 1 h −1 Mpc) of ξ(s) are quite poorly fit. With the exceptions of w p (r p ) and mcf(s), some scales of all other clustering statistics seem to contribute significantly to the tension as well. For −21, while n gal is well-fit, all other clustering statistics appear to be poorly fit on at least one scale. Overall, it is clear that the tension arises due to the inability of our five-parameter model to jointly fit all of these observables.
These results suggest that extensions to the HOD may be necessary in order to match the clustering of SDSS. If our model requires additional features (e.g., assembly bias, velocity bias) which are not included in our fiveparameter HOD, then our posterior results may have significant systematic errors (Zentner et al. 2014). In particular, such errors may be exacerbated in the case of the OPT chains where we include clustering statistics which are sensitive to effects beyond those included in our form of the HOD. Attempting to fit to the measurements of clustering statistics on SDSS with an inadequate model may cause a bias in our posterior results. Indeed, looking again at Figure 7, this effect may explain the offsets we observe between several of the parameter constraints of the NWG and OPT chains. Without including these additional features in our HOD, however, we cannot verify this assertion.
We also consider the impact of our decision to adjust our observational data in order to account for the errors due to our use of the nearest neighbor correction (Section 4.3). Particularly, we wish to know how this treatment affects our posterior results and general con-clusions about the success of our five-parameter HOD model. Therefore, using the OPT observables, we also run chains on the unadjusted datasets for each sample and compare the results to our chains on the adjusted datasets. For −19, we find shifts in the median parameter values of 1σ and 2σ for central and satellite parameters, respectively. Additionally the p-value of the best-fit HOD point is even lower for the unadjusted chain. For −21, we find the opposite trends: shifts in the median parameter values of 2σ and 1σ for central and satellite parameters, respectively. The p-value of the best-fit HOD point is slightly larger, but this change is only enough to decrease the tension from 4.1 to 3.5 σ. Thus, our treatment of fiber collisions has little impact on the general conclusions we draw about the ability of our five-parameter HOD to match the clustering of SDSS.
Before continuing, we address one last final detail concerning our results. When running an MCMC, we are assuming that our fiducial matrix is not too dissimilar from the matrices we would construct from each HOD point in our chain. The HOD parameters we use to construct our fiducial matrix, however, are far from the locations of our posterior (2 σ) contours for both the −19 and −21 OPT chains. We consider performing an iterative procedure in which we reconstruct the covariance matrix from the best-fit HOD point of the OPT chain and then re-run the OPT chain with this new matrix. This procedure could be repeated until some sort of convergence criteria is met. To gauge whether or not this method is necessary, in Appendix D, we test whether two chains run on the same mock galaxy catalog, but differing in the HOD used to construct the matrix, produce similar posterior results. We find that the results are largely similar and thus do not perform any iteration for our SDSS results.

Modifying the halo mass function
The use of the HOD framework as we have applied it is predicated on an assumption that our dark matter only simulations produce the correct halo mass function (HMF). Hydrodynamic effects, however, can alter the masses of dark matter halos (e.g., Vogelsberger et al. 2014;Schaller et al. 2015;Springel et al. 2018;Beltz-Mohrmann et al. 2020;Beltz-Mohrmann & Berlind 2021). By not including baryonic physics, our simulations potentially have incorrect HMFs compared to that of the real Universe.
Beltz-Mohrmann & Berlind (2021) (BM21 hereafter) investigate the differences in HMFs between matched DMO and hydrodynamic simulation boxes in EAGLE, Illustris, and IllustrisTNG (TNG hereafter). For halos at z = 0, they find that in each simulation, stellar feedback generally reduces the masses of low mass halos ( 10 11 h −1 M ), while AGN feedback generally reduces the masses of high mass halos (between 10 12 and 10 13 h −1 M ) compared to their DMO counterparts. However, the exact effect that feedback has on the halo masses differs dramatically from one simulation to the next.
Based on each of these simulations, BM21 identify formulae that can be used to "correct" the masses of halos from a DMO simulation in order to reproduce the HMF from a given hydrodynamic simulation. They provide these corrections for several different redshifts and halo definitions. In this work, we apply the z = 0, M vir corrections to modify the masses of the halos in our DMO simulations. We then re-run our OPT SDSS chains three times, using halos modified according to the prescriptions for EAGLE, Illustris, and TNG. In each of these chains we use the same covariance matrix as in the previous section with parameters listed in Table 4.
We show the results of these chains in Figure 10, with the same general format as Figure 7. In each panel, the EAGLE results are shown in orange, Illustris in blue, and TNG in green. "DMO" (red) refers to the chains run using the original halos (replicated from Figure 7). For both samples, we see that some of the HMF modifications produce significant shifts in the median values of the HOD parameters. For the −21 sample, we see significant shifts in logM min and logM 1 and smaller shifts in σ logM and α. For the −19 sample, we see significant shifts in all parameters except σ logM .
In the −21 sample, we can see large decreases in the values of logM min and logM 1 in the Illustris chain. These shifts are due to the fact that the Illustris correction reduces the masses of all halos above 10 12 h −1 M by 15%. The EAGLE and TNG corrections only slightly reduce the masses of halos between 10 12 and 10 13 h −1 M but have little to no impact on the highest mass halos. This explains the slight decrease in logM min and the lack of change in logM 1 for these chains.
In the −19 sample, we can see a large decrease in the value of logM min in the EAGLE chain but little change in the Illustris and TNG chains. This effect is because the EAGLE correction reduces the masses of halos between 10 11 and 10 12 h −1 M by up to 20%, while the Illustris and TNG corrections have little impact in this regime. Like in the −21 sample, the Illustris correction leads to the biggest change in the value of logM 1 . Unlike −21, however, we see some shifts in logM 1 for EAGLE and TNG because we are in a regime in which the halo masses are slightly reduced by these corrections. Perhaps the most significant result of Figure 10 is that there is no shift in the value of σ logM for the -19 sample and very little shift for the -21 sample. Thus, our claim in Section 7.1 concerning the detectable difference in σ logM between the two samples holds even when considering the impact of baryons on the HMF.
While the changes in logM min and logM 1 are fairly straightforward, the small changes in α in both samples (and in σ logM in the −21 sample) are more difficult to intuit. This demonstrates the fact that the halo mass corrections in BM21 are not simple parameter shifts but complex functions of mass that impact all of our HOD parameters.
Regardless of the reasons for the observed shifts, with the increased constraining power afforded by our numerical approach, the differences in recovered HOD parameter values are generally not robust with respect to the uncertainty in the HMF due to baryonic physics. We certainly do not believe that these three hydrodynamic simulations capture the full uncertainty in baryonic physics or that any of them necessarily result in more accurate HMFs than are produced by DMO simulations. This exercise is simply meant to emphasize that our HOD parameter constraints are subject to baryonic effects. Figure 10. Results of chains run on SDSS using halo catalogs with masses modified according to the prescriptions of BM21. The shifts in logMmin and logM1 can be understood by considering the global number density, while the shifts in σ logM and α are in general more complex (see main text for more details).
Despite the observed changes in the posterior contours, the tension that we find between SDSS and our best-fit HOD model remains largely unchanged for both samples. In the -19 sample, the tension that we find between SDSS and our best-fit HOD model applied to DMO simulations (4.5 σ) remains the same after the Illustris correction, and increases slightly (to 4.6 σ) after the EAGLE and TNG corrections. In the -21 sample, the original tension (4.1 σ) remains the same after the Illustris correction, increases slightly (to 4.2 σ) after the EAGLE correction, and decreases slightly (to 3.9 σ) after the TNG correction. In any case, we can conclude that the tension that we find between SDSS and our standard HOD model is not significantly reduced or exacerbated by the halo mass corrections we employ here.
In this section, we have only demonstrated the effects of correcting the mass function for the simple fiveparameter HOD we use in this work. As we look to constrain additional features of the galaxy-halo connection (e.g., assembly bias) in the future, we must continue to consider the impact of baryonic physics on our posterior results.

CONCLUSIONS
In this work, we explore the ability of a variety of clustering statistics to jointly constrain the connection between galaxies and the dark matter halos in which they live. We develop an algorithm to choose a set of "optimal" observables based on estimates of their joint constraining power. We choose these observables from a superset of multiple scales of the projected correlation function, the redshift-space correlation function, the group multiplicity function, the group velocity dispersion, the mark correlation function, and counts-in-cells statistics. In particular, for two different volume-limited SDSS samples, we employ a Markov Chain Monte Carlo (MCMC) method to constrain a simple five-parameter halo occupation distribution model for connecting galaxies to dark matter halos. We also explore how robust our results are with respect to the influence of baryonic physics on the halo mass function. Our main results are summarized below: • Compared to a similar analysis using only n gal , w p (r p ), and n(N ), when using our chosen set of observables, we are able to achieve tighter constraints on all HOD parameters for both the −19 and -21 SDSS volume-limited samples. Specifically, we tighten the constraints on the central parameters (logM min and σ logM ) by factors of ∼5.8 and ∼3.5 for the −19 sample and ∼4.3 and ∼6.1 for the −21 sample. We tighten the constraints on the satellite parameters (logM 1 and α) by factors of ∼1.8 and ∼2.0 for the −19 sample and ∼1.6 and ∼1.7 for the −21 sample. The fifth parameter, logM 0 , is unconstrained when using n gal , w p (r p ), and n(N ) but is well-constrained for the −19 sample when using our chosen observables. (See Figures 7 and 8, as well as Table 6.) • Our tighter constraints have yielded high precision measurements of σ logM for the −19 and −21 samples, an informative parameter which quantifies the scatter in halo mass at fixed luminosity and is related to the scatter in the stellar-to-halo mass relation (see Figure 7 and Table 6). In particular, we can clearly distinguish between the low and high scatter for the -19 and -21 samples, respectively. • For the −19 sample, we find a significant shift in the posterior distributions when using our chosen observables, compared to when we use only n gal , w p (r p ), and n(N ). This shift may be because we incorporate in our modeling procedure new clustering statistics that depend on additional HOD features not included in our five-parameter model, such as assembly bias, velocity bias, etc. (see Figure 7). • We find substantial tension between SDSS and our best-fit HOD model. This tension is slightly stronger for the −19 sample (4.5 σ) than for −21 (4.1 σ). This is in contrast to when we use only n gal , w p (r p ), and n(N ) to constrain the model, which exhibits no tension. Our inability to jointly match the clustering of SDSS suggests the need for a more flexible HOD model. (See Table 7.) • We find that our joint and marginalized constraints are affected by the impact of baryonic physics on the masses of dark matter halos. The prescriptions of Beltz-Mohrmann & Berlind (2021) which we use to modify our halo masses, according to the hydrodynamic physics of EAGLE, Illustris, and IllustrisTNG, have varying degrees of impact on logM min , logM 1 , and α for the two samples, though σ logM appears relatively stable, particularly for -19 (See Figure 10.) These results demonstrate the power of our numerical modeling methodology to constrain the galaxy-halo connection. By directly populating dark matter simulations and carving out realistic mock galaxy catalogs, we are able to measure clustering statistics in the exact same way on both our SDSS dataset and our model, allowing us to employ less widely used clustering statistics. Further, by choosing a combination of different physical scales for each clustering statistic, we avoid using observables which introduce more noise than information into our posterior results. Using this numerical approach, we are also able to quantify and minimize the error associated with both our model observables (Appendix B) and our covariance matrix (Appendix C). We are also able to verify that an MCMC run on a mock galaxy catalog using our chosen set of observables is able to recover the true HOD parameters (Appendix D), suggesting that the remaining error in our methodology has little impact on our posterior results.
Potentially the largest source of error in this work is our treatment of fiber collisions (Appendix A). Lacking an accurate model of fiber collisions, we instead decide to adjust our SDSS measurements to account for the error arising from only using nearest neighbor corrections. While we believe that our estimation of this error and use of the resulting adjusted measurements are an improvement over nearest neighbor, this approach still may not yield accurate clustering statistics. Unfortunately, there is no way to verify its accuracy without actually obtaining the redshift of every galaxy in SDSS. To get some handle on the impact of our treatment of fiber collisions, however, we also run an MCMC using the unadjusted SDSS measurements. We find that the results have little impact on our general conclusions (see Section 7.1).
Overall, compared to S18, we have succeeded in both achieving tighter constraints on the parameters of our HOD model and in heightening the tension between our model and SDSS. The failure of our five-parameter HOD model to fit the clustering of SDSS for both of volumelimited samples we explore suggests that the assumptions of the standard HOD model may be incorrect, a result in line with several recent works. For instance, halo occupation may depend on secondary properties of the halo, such as age or concentration, in addition to mass (i.e., there may be assembly bias; Zentner et al. 2014Zentner et al. , 2019Vakili & Hahn 2019;Beltz-Mohrmann et al. 2020;Hadzhiyska et al. 2020). Additionally, the number of satellite galaxies in a halo of a given mass may not follow a Poisson distribution (Boylan-Kolchin et al. 2010;Mao et al. 2015;Jiménez et al. 2019). Finally, central galaxies may not always reside at the center of the halo and move with the mean velocity of the halo (i.e., there might be central spatial and/or velocity bias), and satellite galaxies may not trace the spatial and velocity distribution of the dark matter (i.e., there might be satellite spatial and/or velocity bias; Van den Bosch et al. 2005;Watson et al. 2012;Piscionere et al. 2015;Guo et al. 2015a,b;Beltz-Mohrmann et al. 2020). In the case of splashback galaxies, satellites may even exist well outside the radius of their halo (e.g., Wetzel & White 2010;Sinha & Holley-Bockelmann 2012).
In the future, we intend to apply our updated modeling methodology (including our algorithm for choosing observables, treatment of fiber collision errors, and halo mass function corrections) to an expanded HOD model including these additional features. In doing so, we hope to further our understanding of the galaxy-halo connection by exploring its more subtle aspects. An additional possibility that could be responsible for some of the tension that we find is that the cosmological model that we assume in our analysis is incorrect. Fully exploring the cosmological parameter space with our mock-based numerical model is computationally challenging. Nevertheless this research direction is worth pursuing because the high precision that comes from multiple clustering statistics could yield interesting cosmological constraints that are independent from those on larger scales.

A. FIBER COLLISIONS
In this section, we estimate the error in our clustering statistics that results from applying the nearest neighbor correction to SDSS. To do so, we make use of the SDSS plate overlap regions. In an overlap region, the same portion of sky is targeted with two (or more) spectroscopic plates. Approximately 1/3 of the SDSS footprint consists of overlap regions, and thus we have spectra for ∼ 1/3 of galaxies that would otherwise be lost due to fiber collisions. We apply the nearest neighbor correction to galaxies in the overlap regions, effectively pretending some of these galaxies were missed due to fiber collisions. We then examine the resulting impact on our clustering statistics.
For each luminosity threshold, we first construct a baseline fiducial volume-limited sample as described in Section 2. We treat this fiducial sample as an ideal case, pretending that all of the redshifts were spectroscopically obtained. In constructing this fiducial sample we must decide how to treat the galaxies in SDSS which lack spectroscopic redshifts. We choose to assign to these galaxies the redshift of their nearest neighbor. We are thus effectively pretending that the nearest neighbor redshifts are the "correct," spectroscopic redshifts 11 .
Next, we construct modified samples in which we apply the nearest neighbor correction to a fraction of galaxies in the overlap regions. We identify a galaxy as belonging to the overlap regions if both of the following conditions are met: (i) it is within 55 of another galaxy and (ii) spectroscopic redshifts were obtained for both galaxies 12 . To construct our modified samples, we first identify friends-of-friends groups of galaxies, considering only those galaxies in the overlap regions. Galaxies are linked together as part of a group if their angular separation is less than 55 .
We next decide whether each galaxy in the group will either keep its own spectroscopic redshift or receive a nearest neighbor correction. We call a galaxy which is selected at this stage to keep its own redshift a "hit" and one which is not a "non-hit." For each collision group, we maximize the number of galaxies which could have received spectroscopic fibers from only one plate. 11 We repeated this test instead dropping the galaxies without spectroscopic redshifts, and the results did not change. 12 Strictly speaking, these galaxies might not all be in SDSS plate overlap regions. In the parent sample we use, redshifts for some galaxies are obtained from other surveys. This distinction, however, does not matter for purposes of the test we describe.
If there exists more than one possible set of galaxies meeting this condition (e.g., a group only of only two galaxies), then we randomly choose one of these sets and designate the galaxies in this set as hits. The nonhits still may or may not be assigned their spectroscopic redshifts, as we'll soon discuss. We wish to investigate the impact of the nearest neighbor correction as we increase the fraction of galaxies affected by it. We therefore only assign a nearest neighbor correction to a fraction of the non-hits. We consider fractions of 1/4, 1/2, and 1. For fractions of 1/4 and 1/2, we randomly re-designate said fraction of nonhits as "misses." The remaining non-hits are then redesignated as hits. All hits keep their spectroscopic redshifts, while the misses are assigned the redshift of their nearest neighbor (which must be selected from among the hits). For a fraction of 1, all non-hits are designated misses. After assigning a galaxy a new redshift, we recompute this galaxy's k-correction using kcorrect v4 3 (Blanton & Roweis 2007) and absolute magnitude.
For each fraction we consider (1/4, 1/2, 1), we build 100 modified parent samples. Each parent sample is different due to the randomness associated with choosing hits and misses. Finally, we construct one −19 and −21 modified volume-limited sample from each modified parent sample. These samples differ from our fiducial cases because we have modified the redshifts and absolute magnitudes of many of the galaxies. Some galaxies which previously made our sample cuts may be dropped in the new samples and vice versa. We compute clustering statistics on each of these modified samples and compare the results to the fiducial case.
In Figure 11, we show the results of this test for our two volume-limited samples (−19 in blue and −21 in red) for one observable: ξ(s) on a scale of ∼ 0.4 h −1 Mpc. The x-axis is the fraction of galaxies in the total sample to which we reassign a nearest neighbor redshift. For example, when we apply the nearest-neighbor correction to half of the non-hits, this results in a modified −21 sample in which approximately 1.7% of the total galaxies receive the correction. The y-axis is the measurement of ξ(s) on one sample divided by the measurement on the fiducial sample. Each small point is the result for one of the modified samples, while the large points are the mean taken across samples of equal correction fraction. Errors in the mean are plotted but are smaller than the point size. We also show the point (0,1) in which no galaxies receive the correction (i.e., the fiducial case). Figure 11. The estimated error in ξ(s ∼ 0.4 h −1 Mpc) due to using the nearest neighbor correction to handle fiber collisions. The measurement ξ fid is made on a fiducial SDSS sample to which nearest neighbor corrections have been applied. We treat this sample as if redshifts for all galaxies (including the ones which received the nearest neighbor correction) were spectroscopically obtained. The measurement ξNN is made on a modified version of the fiducial sample to which we additionally apply the nearest neighbor correction to galaxies in the SDSS plate overlap regions, reassigning their redshifts. We give ξNN/ξ fid as a function of the fraction of galaxies in the overall sample to which we have reassigned nearest neighbor redshifts. Thus, this fraction for the fiducial sample is zero. The small points represent different realizations of which galaxies receive the correction, while the large points give the mean for a fixed fraction. We fit a line to the means (solid line) and extrapolate (dashed line) to estimate the error for our actual SDSS catalogs (star). Results are shown for −19 and −21 in blue and red, respectively. We perform this same procedure for all of our observables, for both samples. (See text for more details.) We can see in Figure 11 that the error in ξ(s ∼ 0.4 h −1 Mpc) increases linearly as we increase the fraction of galaxies affected by fiber collisions. We estimate this relationship by fitting a line to the means, including the point (0,1) when fitting. This relationship is shown for both samples with the solid lines in Figure 11. Our actual SDSS volume-limited samples, however, contain a higher fraction of galaxies affected by fiber collisions than we can approximate using overlap regions. When building our modified samples, applying the nearest neighbor correction to 100 percent of the non-hits produces −19 and −21 samples in which approximately 2.3 and 3.2 percent of galaxies receive the correction, respectively. In SDSS, however, approximately 5.4 and 7.7 percent of galaxies actually receive the correction for −19 and −21, respectively. To estimate the error for our SDSS samples, we extrapolate the linear fit to the appropriate fraction. We show this extrapolation with a dashed line in Figure 11 and mark the corresponding value of ξ/ξ fid with a star. We call this value c adjust .
Remarkably, the relationship between the error in an observable and the fraction of galaxies affected by fiber collisions is linear for all the observables we consider. We perform the linear fit and extrapolation described above to obtain a value of c adjust for each observable. This value is an estimate of the ratio of the measurement of an observable in our SDSS sample (given the fraction of galaxies which did not receive spectroscopic fibers and thus required the nearest neighbor correction) to the measurement we would make if we had spectroscopic redshifts for all galaxies. Therefore, we can multiply the value of an observable measured on the fiducial SDSS sample by 1/c adjust to obtain an adjusted measurement.
In Figure 12, we show the effect that fiber collisions have on the clustering statistics we use in our analysis. For each observable we take the difference between the measurement on SDSS before adjusting the observables (D NN ) and the measurement after adjusting the observables (D adjust ). We divide this quantity by the cosmic variance uncertainty and plot the result in Figure 12. Results are shown for −19 in the top panel and −21 in the bottom. Each point corresponds to a specific observable, indicated on the x-axis. We indicate the observables we use in our analysis with larger, darker points (see Sections 6 and 7). In general, the nearest neighbor correction works better for the −19 sample than it does for −21. This may be due to the higher fraction of galaxies which received the correction in −21 compared to −19. Although the correction does well for n gal , w p (r p ), and VPF(R), it is inadequate for the small scales of ξ(s) and several scales of n(N ), σ v (N ), mcf(s), and SPF(R).
The results of this test indicate that a failure to adjust the clustering statistics measured on SDSS catalogs that use nearest neighbor corrections can result in significant errors (∼ 1−4 σ obs ) for many of our observables. Therefore, in this paper, we use the adjusted measurements for all observables. We provide both the original nearest neighbor measurements and the adjusted measurements in the machine-readable format. We do not claim that using these adjusted measurements is a complete solution to the issue of fiber collisions. It is simply the best solution we could establish, given the clustering statistics we wish to use in this paper. The assumptions of this test certainly warrant more exploration in future work, but our present goals are not to establish a complete method for dealing with fiber collisions. In any case, in Section 7.1 we perform our analysis on the SDSS datasets with and without the adjustments. While there are some shifts in the posterior results, our general con- Figure 12. Comparison between the original SDSS measurements and the adjusted measurements which account for the error due to using the nearest neighbor correction to handle fiber collisions. We show results for both the -19 (top) and -21 (bottom) volume-limited samples. Each clustering statistic is indicated on the x-axis with a color-coded label. For each observable, the y-axis gives the difference between DNN, the unadjusted measurement on the SDSS catalog (which employs nearest neighbor corrections), and D adjust , the adjusted measurements we use in our analysis. We plot this difference in units of the cosmic variance uncertainty, σ obs . The y-axis thus can be viewed as what happens to the measurement of an observable when we apply nearest neighbor corrections in lieu of having the actual redshifts of all galaxies. The larger points mark the observables we use in our chains on SDSS (Section 7.1).
clusions about the goodness of fit of our five-parameter HOD model do not change.

B. ACCURACY OF MODEL OBSERVABLES
In this work, we expand the methodology of S18 using four additional clustering statistics: the redshift-space correlation function ξ(s), the mean group velocity dispersion σ v (N ), the mark correlation function mcf(s), and counts-in-cells P N (R). For each of these clustering statistics, we must decide on the most accurate way to measure the model observables (see Section 5). The ideal numerical measurement of the model observables is the mean measurement across many high-resolution mock galaxy catalogs constructed from the same HOD and differing only due to cosmic variance. Obtaining such a large suite of high-resolution simulations and using them to explore a large number of HOD points in an MCMC both have a prohibitively high computational cost. We instead have just two high-resolution boxes which we use to estimate the model observables. Still, we have a choice to measure each clustering statistic us-ing either full boxes or mocks carved out from the boxes. Additionally we may use just one box or both boxes. We explore this choice in this section.
While using more volume results in more precisely measured model observables, these measurements may also be more biased (i.e., further from the desired mean than the measurements made on a smaller volume). Because we do not have a large number of high-resolution simulations, we cannot use our high-resolution boxes to test how close a given method of measurement is to the mean measurement across many mocks. The closeness of a given box's measurements to the desired mean measurements is driven by the cosmic variance error of that box. While the measurements themselves are affected by resolution, we do not expect that the cosmic variance of a box is affected by resolution. Therefore, we can use our suite of 100 low-resolution simulation boxes (two of which have the same initial condition seeds and thus the same large-scale density modes as our high-resolution boxes) to determine the most appropriate method for measuring each clustering statistic. We assume that the results of this test would be the same if we had a set of 100 high-resolution boxes with the same seeds as our low-resolution boxes.
To make this choice, we first select a fiducial HOD 13 and populate each of our 100 low-resolution boxes with galaxies. From these boxes we create 400 mock galaxy catalogs and measure each of our clustering statistics on these mocks. We use the term "true mean" to refer to the mean measurement of an observable across these 400 mocks. We next consider different estimations of the model observables using the two boxes for which we have high-resolution counterparts (see Table 2). For all clustering statistics, we consider (1) the mean measurement of four mocks from box 1, (2) the mean measurement of four mocks from box 2, and (3) the mean measurement of eight mocks from both boxes. For n gal , w p (r p ), ξ(s), and P N (R), we also consider (4) the measurement on box 1, (5) the measurement on box 2, and (6) the average of these two measurements. We do not consider box measurements for n(N ), σ v (N ), and mcf(s) because these statistics all depend on galaxy group identification, a process which is affected by survey geometry (Berlind et al. 2006;Campbell et al. 2015). When we measure statistics on a full box, we first distort galaxies along the z-direction using the distant-observer approximation. Once we have each of our estimates of the model observables, we compare them to the "true mean" in order to determine the most appropriate method of measurement.
This test is illustrated for one observable, SPF(R) on a scale of 2 h −1 Mpc, for the −19 sample in Figure 13. The black histogram shows measurements of SPF(R) on 400 mocks, while the vertical black line shows the "true mean." Each measurement x obs in this plot is standardized via the equation x std = x obs − x obs σ obs , where x obs and σ obs are the mean and standard deviation of measurements across 400 mocks, respectively. The estimates using boxes and mocks are shown with the multi-colored dots, as is indicated in the legend. Each time we populate a halo catalog with galaxies, the exact numbers of galaxies in each halo and the placement of these galaxies can vary (see Section 4.2), resulting in a stochastic scatter associated with each estimate of the model observables. In practice, this process is controlled by a random number generator and is seeded with a "population seed." To get a handle on this scatter, we repopulate each of our two low-resolution boxes 100 times, each 13 Here, we show results using the HOD parameters in Table 4. We have repeated this test using different HOD parameters, and the results do not change. The black histogram shows the distribution of measurements on 400 low-resolution mock galaxy catalogs all created with the same HOD, and the black dashed vertical line shows the mean of this distribution. The colored points show estimates of this mean when using fewer boxes/mocks, as is indicated in the legend. The error bar on these points comes from the stochasticity of the process of populating dark matter halos. All measurements are standardized using the mean and standard deviation of the 400 mocks.
time using a different population seed. We then recalculate our estimates of the model 100 times. Lastly, we standardize each of these 100 measurements by subtracting x obs and dividing by σ obs . For each type of estimate, we show the mean and standard deviation of these 100 measurements in Figure 13, with vertical offset added for display purposes. Standardizing allows us to easily compare and interpret the accuracy of model estimates across all of our observables. For example, we can see in Figure 13 that using the 4002 mocks to estimate this scale of SPF(R) results in an estimate that is ∼ 1 σ obs (i.e., the cosmic variance uncertainty) lower than we desire, while using the 4022 mocks produces a much more accurate estimate. Naively, if given the choice of averaging measurements from eight mocks or from four mocks, we would choose to average eight mocks. However, we can see from this plot that the 4002 mocks happen to be outliers for this particular observable. Thus, including them as part of the eight-mock average results in a more biased model observable than just using the four mocks from box 4022. On this standardized scale, we will refer to the difference between an estimate and zero as the "model deviation" and the associated scatter as the "model scatter." Looking at Figure 13, there are several good choices for how to measure this specific scale of SPF(R). To maintain some simplicity in our methodology, however, Figure 14. The standardized error in our estimates of the model, after choosing a measurement method for each clustering statistic. The top and bottom panels show the model deviation for the −19 and −21 samples, respectively. Each clustering statistic is given a different color, as is indicated on the horizontal axis. The error bar associated with each point comes from the stochastic population process and reflects the precision of the model estimate. The more opaque points are the observables we ultimately use in our analysis (see Section 6.3).
we choose one form of measurement per clustering statistic (not per observable). Thus, we must take into account how a given method performs for all scales of a clustering statistic. We treat VPF(R) and SPF(R) as one clustering statistic when making this choice. While we certainly desire precise estimates of the model observables (small model scatter), our priority is to choose estimates that are accurate (small model deviation) across all scales. We make this choice simply by visually comparing the different methods as opposed to establishing any sort of comparison metric. We summarize the ultimate choices we make for how to measure each clustering statistic in Table 3 (see Section 5 for a description).
In Figure 14 for each observable we show the model deviation, given our choice of how to measure said observable. We show results for the −19 sample in the top panel and for the −21 sample in the bottom. The different clustering statistics are assigned different colors as indicated by the label on the x-axis. Each point represents a clustering statistic measured at a different scale, ordered in increasing physical size. The more opaque points represent the observable we use in our analysis. For example, the dark green points in both panels show the different scales of n(N ). For the −19 sample we can see that the smallest scale has a model deviation of about -0.25 σ obs . Thus, we expect our estimate of the model observable to be about 0.25 σ obs lower than it should be on average. The scatter in this estimate is about 0.2 σ obs , meaning that our estimate of this model observable is between 0.05 and 0.45 σ obs lower than it should be for 68% of the HOD points we evaluate in our MCMC 14 . We display the observables we use in our final SDSS analysis (Section 7.1) in bold.
While we could cull observables by establishing thresholds for model deviation and scatter, we do not know what impact these errors have when they are included in a chain. Establishing any sort of cutoff then seems rather arbitrary, especially when considering that the model deviation is generally small ( 0.5 σ obs for almost every observable). We are content to quantify this error and choose methods of measurement that mini-mize its effects. In any case, in Appendix D we show that we are able to recover the correct HOD when running a chain on a mock galaxy catalog using the same observables we use in our SDSS analysis (Section 7.1). This result suggests that the model deviations shown in Figure 14 do not significantly affect the locations of our posterior results.

C. ERRORS IN COVARIANCE MATRIX
While using two high resolution boxes is sufficient for estimating the model observables, the covariance matrix requires many mock galaxy catalogs. Approximate methods, such as jackknifing and bootstrapping, have been shown to produce biased results (Norberg et al. 2009). We follow the approach of S18 in using many low-resolution mock galaxy catalogs to create a covariance matrix. While we have made efforts to improve the covariance matrices we use compared to S18 (see Section 5), our covariance matrices have four identifiable sources of error: 1. The limited number of mocks. 2. The stochastic variation due to population seed. 3. The use of a fixed fiducial HOD. 4. The resolution of the simulations.
We discuss the relative importance of the first three of these points in Appendix C.1 and the last in Appendix C.2. We examine the noise only in the cosmic variance uncertainty σ obs of each observable, ignoring off-diagonal elements of the matrix.

C.1. Comparing relative errors
Assuming Gaussian errors, for a sample of size N and true standard deviation σ, the error in a calculation of the sample standard deviation is σ/ √ 2N . For N = 400 mocks, this gives a fractional error of ∼0.035. We have verified that the distribution of mock measurements for all observables we consider is roughly Gaussian and so this fractional error holds for all values of σ obs .
The elements of our covariance matrix have an additional error due to the stochastic process of populating dark matter halos, which is controlled by the "population seed." Using a different population seed will produce a different matrix. To quantify this error, using our fiducial HOD we remake our 400 mocks 100 different times with 100 different population seeds. This gives us 100 estimates of σ obs . For a particular observable, the fractional error in σ obs is given by where σ obs,i is the i th estimate of the standard deviation of this observable and σ obs is the average of 100 estimates of the standard deviation of this observable. Lastly, there is an error stemming from our decision to not vary the matrix when we consider new HOD parameters in our chain. Ideally, each time we evaluate a new set of HOD parameters, we would reconstruct the covariance matrix from 400 mocks built using those HOD parameters. Because building and measuring clustering statistics on 400 mocks for the ∼500,000 HOD points we evaluate in a chain is computationally infeasible, we instead keep our covariance matrix fixed. To investigate the error we are making in using a fixed matrix, we calculate the fractional error δ(σ obs ) again. This time, however, instead of varying the population seed while keeping HOD fixed, we build 100 different matrices, each with a different HOD. To choose the 100 HODs, we randomly sample points from our OPT chains (see Section 7.1). A broader range of HODs would likely give us a broader error, but we only care about this error insofar as it affects our results. Thus we only consider the range of HODs explored in our analysis.
In Figure 15, we show these three contributions to the fractional error in σ obs . The format of this figure is similar to Figure 14, with clustering statistics assigned different colors and each bar representing a different scale, ordered in increasing physical size. The height of the open bar gives the fractional error in σ obs due to not varying HOD, while the filled bar gives the error due to the stochastic population process. The Poisson noise associated with the number of mocks we use is the same for all observables and is indicated with the horizontal dashed line. The clustering statistics we use in our SDSS analysis (Section 7.1) are indicated with greater opacity in the shaded bars.
The fractional errors in σ obs are all smaller than 10% for both samples. For most observables, the dominant contribution to the noise is the number of mocks we use when constructing the matrix. A few scales of observables (e.g., small scales of w p (r p )) show an exception to this trend, with the dominant error being due to not varying HOD. For many observables (e.g., mcf(s)), there appears to be an "equal" contribution from all sources of noise. It is worth noting that the HOD error effectively contains the stochastic population error because we choose to supply unique population seeds for each HOD 15 . Therefore, in instances where these two errors 15 The way our population process works, even if we supply the same seed, different HODs require the generation of different amounts of random numbers, which affects the random number sequencing. Thus, we simply choose to supply different random numbers. Figure 15. The fractional error in the cosmic variance uncertainty due to various processes. The horizontal line gives the shot noise due to using only 400 mocks to construct the matrix. The height of the hollow vertical bars shows the error from reconstructing the matrix from 100 different HODs, while the height of the filled vertical bars shows the error from varying the "population seed" (i.e., the seed which controls the numbers of centrals/satellites assigned to each halo) 100 times. The clustering statistics are colored and ordered in the same way as in Figure 14. The more opaque bars show the observables we actually use in our analysis. For both samples shot noise is the dominant source of error for most observables, though the error from not varying HOD is dominant for many scales.
appear equal, the stochastic population error may in fact be dominant.

C.2. Resolution effects
An assumption in our modeling procedure is that our covariance matrix is not significantly affected by resolution. In other words, we assume that, if we were to build the same covariance matrix using 400 high resolution mocks, the differences in matrices would be small. We cannot accurately test this assumption without actually having the 100 high-resolution boxes, which would negate the need to test the assumption. However, what we do have is a matching set of two low-resolution and high-resolution boxes. These boxes are identical in their density fields and only differ in resolution. We can therefore use these boxes to test the impact of resolution on our matrix.
To test the effect of resolution, we first divide each box into 27 equal-volume sub-boxes. In the cases of both samples, the volume of each of these sub-boxes is about half the volume of SDSS. With two seeds, this gives us 54 low-resolution and 54 matching high-resolution sub-boxes. Using a fiducial HOD (Table 4), we populate each of these sub-boxes with galaxies and measure clustering statistics. For each observable, we compute a high-resolution standard deviation (σ HD ) and a lowresolution standard deviation (σ LR ) across the 54 subboxes.
As discussed in Appendix B, the process of populating dark matter halos is stochastic and is controlled by a "population seed." This stochasticity introduces noise into any single calculation of the quantity σ LR /σ HD . Therefore, we repeat the population process of both lowresolution and high-resolution sub-boxes with 100 different population seeds. We measure clustering statistics for each of these 100 sets of 54 sub-boxes and get 100 measurements of σ LR and σ HD for each observable. We then average to eliminate the noise due to this stochasticity.
In Figure 16, we show the results of this test for both −19 (top) and −21 (bottom). As in previous figures, different clustering statistics are assigned different colors as indicated by the labels on the x-axis. Each point represents a clustering statistic measured at a different scale, Figure 16. The error in the cosmic variance uncertainty due to resolution. For each observable we calculate the standard deviation across 54 matched high ("HD") and low ("LR") resolution "sub-boxes." This process is repeated 100 times while varying the "population seed" to create a distribution for high and low resolution values of σ. We show here the ratio of means of these distributions. Error bars (barely visible) show the propagation of the errors in the mean of these two distributions. ordered in increasing physical size. For each observable we plot the mean of the 100 values of σ LR divided by the mean of the 100 values of σ HD . We henceforth refer to this ratio as the "resolution error." The error bar (barely visible) associated with each point is found via propagation of the errors in the means of the two distributions. The optimal observables we use in our analysis (Section 7) are indicated with greater opacity.
For the −19 sample, all observables have a resolution error smaller than 6%. For the −21 sample, the resolution errors are larger, but most errors are still less than 10%. The observables most affected by resolution are the small scales of w p (r p ) and ξ(s) in the −21 sample, which show a 10−20% error. The larger resolution error in the −21 sample is due to the fact that halos of mass M min (10 12.72 h −1 M ) in Carmen have fewer particles than halos of mass M min (10 11.54 h −1 M ) in Consuelo. Furthermore, the larger value of σ logM in Carmen leads us to place galaxies in even lower mass halos that are even less resolved. While the resolution error appears in general to be more significant than the errors discussed in Appendix C.1, we emphasize that overall these errors are quite small and likely have very little impact on our analysis.
Given this result, we do not make any attempts to correct for the resolution error. We are in general against adding complexity to our modeling procedure unless it seems necessary and well-motivated. Ultimately our covariance matrix is built from mocks which have a different volume and geometry than the sub-boxes we use to test resolution. We are not confident that, given these differences, the results of this test would be the same in the case of SDSS-like mocks. Additionally, larger cosmic variance uncertainties result in lower χ 2 measurements and in general make it more difficult to rule out incorrect models. If, as this tests suggests, we are using larger cosmic variance uncertainties, then the results of our MCMC will allow a greater space of HOD parameters. We are content to have artificially broad constraints, as opposed to making a correction that might result in artificially tight constraints.
Ultimately, we are concerned with what impact the resolution errors will have on the results of our chains. In Appendix D we show that we are able to recover the correct HOD when running a chain on a mock galaxy catalog using the same observables we use in our SDSS analysis (Section 7.1). This result suggests that the errors in the matrix due to resolution and the processes Figure 17. Similar to Fig. 7, except that these chains are run on mock galaxy catalogs instead of the SDSS. The left and right panels show results for the −19 and −21 mocks, respectively. Both the green contours and the purple contours come from chains run on the same mock galaxy catalog using the same observables as in our OPT SDSS chains. For the "Fid 1" chain (green), the HOD used to make the covariance matrix is not the same HOD used to generate the mock. For "Fid 2" (purple), the matrix and mock are made using the same HOD. The cross-hairs mark the true HOD parameters of the mock galaxy catalog. We are able to recover the true HOD when using our OPT observables, and the difference in covariance matrix has little impact on the results. The HOD parameters for sets Fid 1 and Fid 2 are given in Table 8. discussed in Appendix C.1 do not have a detectable impact on the locations of our posterior results.

D. MOCK CHAINS
After choosing the set of observables for our SDSS chains, we also used the process of importance sampling to conclude that we could reasonably recover the true HOD parameters when running a chain on a mock galaxy catalog (see end of Section 6.3). Still, because all importance-sampled chains contain some noise, we desire a more robust test of the ability of a chain run using our chosen set of observables to recover the true HOD parameters. In this section we explicitly test whether or not we can recover the true HOD parameters when running a chain on a mock galaxy catalog. We also investigate the effect that changing the HOD used to build our fiducial covariance matrix has on our posterior results.
To perform this test, we run an MCMC on a mock galaxy catalog using the same set of observables we use in our SDSS OPT chains. We construct both −19 and Note-For each sample, we run two chains on a mock galaxy catalog generated using the HOD parameters given by Fid 2. In one of these chains, we construct our covariance matrix using Fid 1, and in the other we construct our matrix using Fid 2.
−21 mock galaxy catalogs from the best-fit HOD values of the SDSS OPT chains 16 . We use these HODs (Fid 2 in Table 8) so that our mock has clustering similar to that of SDSS. When running this MCMC, we use the same covariance matrices we used in our SDSS OPT chains (Fid 1 in Table 8). These choices are made so as to best replicate the OPT chains of Section 7.1 but in the case where we know the true HOD parameters.
In Figure 17 we show the results of this test with the green contours (labeled Fid 1 in reference to the matrix). This figure has the exact same general formatting as Figures 7 and 10. The cross-hairs indicate the true HOD of the mock galaxy catalog. We can see that, for both -19 and -21, the true HOD lies inside the joint 68% probability region for both central and satellite parameters. Where exactly the truth should lie is heavily dependent on the cosmic variance of the mock galaxy catalog. Thus, an offset of the middle of the contours from the truth, such as we observe here, is completely within reason. Therefore, we conclude that when using our set of optimal observables, we are able to recover the true HOD parameters.
In Section 7.1, we noted that the HOD used to construct the fiducial covariance matrix was far from the location of our posterior constraints. To investigate what impact this has, we also run (for both −19 and −21 mocks) a second chain in which we use a covariance matrix constructed from the same HOD as the mock galaxy catalog (Fid 2). These results are shown in Figure 17 with the purple contours (labeled Fid 2). When comparing chains run on the same mock but with different matrices (i.e., purple and green contours), there are very minor differences in the posterior results. Critically, the shifts seen for all parameters are not enough to affect the recovery of the true HOD parameters. We considered performing an iterative procedure in which we would reconstruct the covariance matrix from the best-fit HOD parameters from the OPT chain and then rerun the OPT chain using this matrix. In light of these results, we conclude that this iterative procedure is unnecessary and that our choice of fiducial HOD has little impact overall on our posterior results.
Additionally, we check whether our model is able to fit the clustering of the mocks and find that this is indeed the case. In all of these chains, the p-value of the best-fit HOD point is between 0.74 and 0.94. Compare this with the SDSS OPT results, where the p-values of the bestfit point indicate that our model is unable to jointly fit the clustering of SDSS. These results suggest that the tension we find in Section 7 is real and not an artifact of our modeling procedure.