Modeling Redshift-Space Clustering with Abundance Matching

We explore the degrees of freedom required to jointly fit projected and redshift-space clustering of galaxies selected in three bins of stellar mass from the Sloan Digital Sky Survey Main Galaxy Sample (SDSS MGS) using a subhalo abundance matching (SHAM) model. We employ emulators for relevant clustering statistics in order to facilitate our analysis, leading to large speed gains with minimal loss of accuracy. We are able to simultaneously fit the projected and redshift-space clustering of the two most massive galaxy samples that we consider with just two free parameters: scatter in stellar mass at fixed SHAM proxy and the dependence of the SHAM proxy on dark matter halo concentration. We find some evidence for models that include velocity bias, but including orphan galaxies improves our fits to the lower mass samples significantly. We also model the clustering signals of specific star formation rate (SSFR) selected samples using conditional abundance matching (CAM). We obtain acceptable fits to projected and redshift-space clustering as a function of SSFR and stellar mass using two CAM variants, although the fits are worse than for stellar mass-selected samples alone. By incorporating non-unity correlations between the CAM proxy and SSFR we are able to resolve previously identified discrepancies between CAM predictions and SDSS observations of the environmental dependence of quenching for isolated central galaxies.


INTRODUCTION
Our theoretical understanding of galaxy formation has advanced significantly with the advent of high resolution N -body simulations that are capable of resolving substructure within dark matter halos. Shortly after the the first such simulations were possible, so called subhalo abundance matching (SHAM) models were devised to take advantage of this newfound ability (Kravtsov et al. 2004;Vale & Ostriker 2004;Conroy et al. 2006). These models place simulated galaxies directly on resolved substructure, and posit that the stellar mass or luminosity of a galaxy is approximately monotonically related to the mass or velocity of the dark matter (sub)halo hosting that galaxy. Despite their simplicity, SHAM mod- * jderose@lbl.gov els are able to reproduce a broad range of galaxy spatial statistics, including projected two-point clustering, conditional luminosity functions, and radial profiles of galaxies within halos (Reddick et al. 2013;Saito et al. 2016;Lehmann et al. 2017). In this work, we demonstrate that SHAM models can also fit redshift-space clustering measurements.
Empirical models such as SHAM now form the basis of many wide-field galaxy survey simulations because of their ability to predict clustering over a broad range of scales and redshifts for a variety of galaxy samples (DeRose et al. 2019;Korytov et al. 2019). The SHAM approach also holds promise as a forward model for cosmological analyses that employ galaxy clustering statistics. Other simulation-based models such as halo occupation distribution (HOD) models (Seljak 2000;Berlind & Weinberg 2002;Bullock et al. 2002) must assume an-alytic models for the number of galaxies in a halo, N , as a function that halo's mass, P (N |M ), and a phase-space distribution of galaxies within halos, P (x, v|M, c), that again depends on the host halo's mass, M , as well as its concentration, c. Such models are able to fit galaxy projected and redshift-space clustering into the one-halo regime (Reid et al. 2014;Zhai et al. 2019;Lange et al. 2021), but usually use five to ten parameters in doing so, making SHAM a potentially more efficient parameterization.
While SHAM has been touted as a highly predictive model for galaxy clustering, questions remain about its ability to fit observables at the statistical precision afforded by modern galaxy surveys. The highest signalto-noise galaxy clustering measurements currently available come from galaxy spectroscopic redshift surveys; to make predictions for the observables from such surveys, models must account for both galaxy positions and line-of-sight velocities. A number of investigations have been conducted into the ability of SHAM to fit redshiftspace galaxy clustering measurements from these surveys. Saito et al. (2016) showed that SHAM combined with a model for stellar mass incompleteness as a function of galaxy color could fit projected clustering measurements from the BOSS CMASS sample. Yamamoto et al. (2015) compared the predictions of a SHAM model using subhalo maximum circular velocity as a proxy for galaxy absolute magnitude to redshift-space clustering statistics measured from the Sloan Digital Sky Survey Main Galaxy Sample (SDSS MGS), and found reasonable agreement, but no quantitative goodness-offit statistics were provided. Guo et al. (2016) performed a more quantitative comparison, investigating the ability of SHAM and HOD models to fit redshift-space distortion (RSD) measurements. They concluded that SHAM models with a single free parameter governing scatter between the SHAM proxy and galaxy absolute magnitude were unable to simultaneously fit projected and redshift-space clustering. Recently Contreras et al. (2021a) found that a SHAM model with parameters governing subhalo artificial disruption and galaxy assembly bias in addition to a parameter governing the scatter between SHAM proxy and stellar mass were able to fit RSD measurements in Illustris TNG Nelson et al. 2018;Marinacci et al. 2018;Naiman et al. 2018;Pillepich et al. 2018), and projected clustering statistics from SDSS. In this work, we perform an extensive study of SHAM's ability to fit RSD measurements from the SDSS MGS, applying quantitative goodness-of-fit metrics and investigating how SHAM extensions such as orphan galaxies and velocity bias affect these fits. The effect of orphan galaxies in SHAM was also explored in Contreras et al. (2021b).
A major drawback of commonly used SHAM approaches are their inability to model galaxy samples that are incomplete in galaxy luminosity or stellar mass. Nearly all samples used for modern cosmology analyses have this property. For example, the main sample used for cosmological constraints in the Baryon Acoustic Oscillation Survey (BOSS) is the CMASS sample (Dawson et al. 2013). CMASS, while approximately complete in stellar mass, has incompleteness that is a function of galaxy star formation rate (SFR) Saito et al. 2016). Galaxy clustering is highly dependent on SFR (Zehavi et al. 2011), so SFR-dependent incompleteness imparts a bias in the CMASS clustering with respect to the clustering of a stellar-mass-complete sample that simple implementations of SHAM cannot model. Upcoming surveys will pose similar or even more severe sample incompleteness problems as many upcoming spectroscopic surveys are designed to target galaxies with strong emission lines (DESI Collaboration et al. 2016;Laureijs et al. 2011;Dore et al. 2019). Thus, methods for incorporating such selections in the SHAM framework must be developed if the approach is to be used to model these samples.
Significant progress has been made in modeling the dependence of galaxy clustering on SFR. Hearin & Watson (2013) introduced an extension to SHAM, dubbed age matching, that correlates galaxy z = 0.1-frame g − r color at fixed r-band luminosity with dark matter subhalo formation time. They showed that age matching reproduces the dependence of projected clustering on galaxy luminosity and color at high accuracy with no free parameters other than the choice of proxies used for the SHAM and age-matching procedures. A similar method was also put forth by Masaki et al. (2013) with comparable results. Generalizations of these first models, now referred to as conditional abundance matching (CAM) models Watson et al. 2015), have been tested against SDSS galaxy-galaxy lensing as well as group and galaxy cluster statistics with encouraging results.
CAM models have also been applied to RSD measurements in the past with little success: Yamamoto et al. (2015) constrained the models presented in Masaki et al. (2013) against SDSS MGS redshift-space clustering finding that none of their CAM models fit the data well when considering chi-squared tests. In the second part of this work we confront CAM with RSD measurements for samples split by stellar mass and SSFR. We focus on two new CAM models, one based on subhalo's accretion rate, and another that uses distance to the near-est host halo above a set mass threshold, examining the goodness-of-fit of these models when used in conjunction with a model for orphan galaxies.
The structure of this paper is as follows. In section 2 we describe the SDSS galaxy samples used in this work and the methodology used to measure RSD statistics including how we account for fiber collisions. In section 3 we describe the simulation employed in this work, which was run at the best fit cosmology from Planck Collaboration et al. (2016). Section 4 introduces the details of our SHAM and CAM models, including how these models are implemented in our simulation. Section 5 presents the results of our SHAM fits to SDSS MGS RSD data and how these fits are impacted by the inclusion of orphan galaxies and velocity bias. In section 6 we fit two CAM models to SDSS MGS RSD measurements as a function of stellar mass and specific star formation rate (SSFR), and use posterior predictive distributions from these fits in order to perform a detailed comparison of these models. In section 7 we summarize our main results and conclude by discussing future directions of investigation.

OBSERVATIONAL DATA AND CLUSTERING MEASUREMENTS
In this work we make use of the NYU Value Added Galaxy Catalog (VAGC) (Blanton et al. 2005), which is constructed from SDSS DR7 (Abazajian et al. 2009). We consider three different volume limited samples: 10 9.8 ≤ M * < 10 10.2 , 10 10.2 ≤ M * < 10 10.6 , 10 10.6 ≤ M * < 10 11.2 , where the redshift limits of these samples are given in table 1, and all stellar masses are quoted in units of h −2 M . Although the redshift ranges have been chosen to provide volume complete samples, we have restricted the most massive sample's upper redshift limit to a lower value than otherwise would be possible so that we do not need to account for redshift evolution of the stellar mass function in our SHAM models. We also subdivide these samples into star-forming and quenched galaxies using measurements of their specific star formation rate (SSFR), where we categorize galaxies with SSFR < 10 −11 yr −1 as quenched. In order to avoid complications from the differences between the imaging used for target selection in the north galactic cap (NGC) and south galactic cap, we limit our analyses to the the NGC.
We measure three different clustering statistics for each sample: projected clustering, w p (r p ), and the monopole and quadrupole moments of the redshift-space clustering signal, ξ 0 (s) and ξ 2 (s). The projected correlation function is given by: where π = s·l |l| , r 2 p = s · s − π 2 , and s = s 1 − s 2 , s 1 and s 2 are the redshift-space coordinates of two galaxies, and l = (s 1 + s 2 )/2 (Davis & Peebles 1983;Fisher et al. 1994). We use π max = 40 h −1 Mpc for all of the w p (r p ) measurements presented in this work. We estimate ξ(r p , π) using the Landy-Szalay estimator (Landy & Szalay 1993): where DD, RR, and DR are the number of data-data, random-random and data-random pairs, normalized by the total number of pairs in each radial bin. All of the pair-counting done in this paper makes use of the Corrfunc library (Sinha & Garrison 2017). We use 12 logarithmically spaced bins between r p = 0.13 − 32.6 h −1 Mpc , and 40 linearly spaced bins in π. We assume the SMDPL cosmology (see section 3) to convert redshift to comoving LOS distance, always setting h = 1. The monopole and quadrupole moments of the anisotropic redshift-space clustering signal are given by: where µ = r p /s is the cosine of the angle between the line-of-sight and s. L is the -th Legendre polynomial, with L 0 = 1 and L 2 = (3µ 2 − 1)/2. ξ(s, µ) is also estimated using the Landy-Szalay estimator. We use 12 logarithmically spaced bins between s = 0.13−32.6 h −1 Mpc and 40 linearly spaced bins between 0 ≤ µ < 1.
Due to the finite size of the SDSS fibers used to transmit light from the focal plane to the spectrographs, spectra of pairs of galaxies closer than the 55 arcsecond diameter of a fiber cannot both be observed. This angular size translates to r p = 0.12 h −1 Mpc at the highest redshift edge of the samples considered in this work. For w p (r p ) we simply restrict our analysis to scales larger than this. For ξ (s) we must be more careful, as these statistics can have contributions from r p < 0.12 h −1 Mpc for s > 0.12 h −1 Mpc . In order to remove this dependence, we first assign all galaxies with missing redshifts their "Nearest Neighbor" redshift as implemented in the NYU VAGC. Additionally, we exclude any (s, µ) bins that have any contribution from r p < 0.12 h −1 Mpc before computing multipoles, i.e.
where µ max = r p,max /s. Reid et al. (2014) showed that this approach is excellent at removing the bias imparted to ξ 0/2 (s) by fiber collisions at the expense of sensitivity to clustering at scales smaller than the fiber collision radius. We estimate covariance matrices for all of the measurements presented here using a jackknife procedure. We use N SIDE = 8 HealPix cells (Górski et al. 2005) as jackknife regions. We assign randoms to these HealPix cells, and exclude those cells that have fewer than 50% of the average number density of randoms in them in order to ensure that our jackknife regions are of equal area. This results in 127 equal-area regions at the expense of removing 7253 galaxies that would have otherwise been included in our samples. We can then compute our covariance matrix as: where x i and x j are two elements of our data-vector, x i,k and x j,k are the same elements measured when leaving the k-th jackknife region out, andx i andx j are the means of those elements averaged over all N = 127 jackknife measurements. We apply a Hartlap correction to our inverse covariance matrices (Hartlap et al. 2007) in order ameliorate biases in our analysis due to noise in the covariance matrix. The largest data vector that we consider in this work, the combination of w p andξ 0/2 for all three stellar mass bins, has length N d = 108, so the Hartlap factor leads to a 87% decrease in constraining power in this case.

SIMULATIONS
In this work we make use of the Small Multi-Dark Planck (SMDPL) simulation (Klypin et al. 2016), an Nbody simulation run using L-Gadget2 (Springel 2005) with 3840 3 particles in a (400 h −1 Mpc) 3 volume and a force softening of Plummer = 1.5 h −1 kpc , yielding a mass resolution of 9.63 × 10 7 h −1 Mpc. Using a simulation with a volume of at least (400h −1 Mpc) 3 is important when analyzing SDSS MGS clustering, otherwise sample variance from the simulation becomes a dominant contribution to SHAM parameter constraints (Lehmann et al. 2017). The simulation was initialized using the Zel'dovich approximation at z = 100 and the best fit cosmological parameters from Planck Collaboration et al. (2016). Halo finding was performed using Rockstar (Behroozi et al. 2013a), assuming a virial overdensity definition (Bryan & Norman 1998) and removing unbound particles from the halo mass estimates. Merger trees were generated for these halo catalogs using Consistent Trees, (Behroozi et al. 2013b), and orphan halos were simulated using UniverseMachine (Behroozi et al. 2019a). Only the z = 0 snapshot is used in this work.

MODELS AND MEASUREMENTS FROM SIMULATIONS
Subhalo abundance matching is a technique that assigns galaxy stellar masses or luminosities to resolved dark matter (sub)halos by enforcing the relation: where Φ(M * > x) is the cumulative number density of galaxies more massive than M * , and n(X h > y) is the cumulative number density of halos where some chosen halo property, X h , often times referred to as the SHAM proxy, is greater than y. This implicitly defines a relation M * (X h ) that can be used to assign galaxy stellar masses to halos as a function of X h . We use the stellar mass function measured in Reddick et al. (2013) based on the data described in section 2. Although this stellar mass function was measured over the redshift range 0.026 ≤ z < 0.067, we have confirmed that using this stellar mass function accurately reproduces the number densities of our two most massive samples whose upper redshift limit extends to z = 0.106.
Equation (6) holds in the case that there is zero scatter in the relation between M * and X h , but in any realistic scenario there is scatter induced in this relation, both due to observational uncertainty in the measurement of M * and because of correlations between M * and variables in addition to X h that have been neglected in our model. We account for these sources of scatter by deconvolving a fiducial amount of scatter, σ log M * |X h (abbreviated as σ log M * for the duration of this work), from Φ(x > M * ), using eq. (6) to determine M r (X h ), and then to adding the same amount of scatter back to the assigned values of M * . The deconvolution procedure assumes p(log M * |X h ) to be log-normal, where σ log M * is conventionally quoted in dex. Behroozi et al. (2010) describes this method in greater detail. σ log M * |X h is then left as a free parameter in the SHAM model to be fit to the data under consideration.

SHAM Proxy
In this work we use v α as the proxy for M * i.e. X h = v α in eq. (6). v α is given by: where v max is the maximum circular velocity of the halo, v vir = GMvir Rvir , is the virial velocity of the halo, and α is a free parameter that governs the concentration dependence of the abundance matching proxy. This proxy is evaluated at the epoch at which the halo attains its peak mass, i.e. v max = v Mpeak in eq. (7), and v vir is computed using M peak .
Using quantities evaluated at this epoch mitigates the effects of mergers on v max , which otherwise can cause large temporary spikes in v max that likely do not contain information relevant to long term stellar mass evolution (Behroozi et al. 2014;Chaves-Montero et al. 2016).

Orphans
It has been shown that subhalos are susceptible to artificial disruption (van den Bosch & Ogiya 2018) even in high resolution N -body simulations such as the one used in this work. Whether a subhalo artificially disrupts or not is strongly dependent on the orbit that the subhalo takes within its host halo, and as such can potentially impart biases on the shape of the clustering signals predicted by SHAM.
The way that we account for this effect in our modeling is two-fold. First, we determine when a subhalo can no longer be tracked by our merger tree algorithm. After this point, we continue to evolve these disrupted subhalos forward in time, modeling their mass and velocity evolution as well as their orbits within their host halo semi-analytically in the manner described in Appendix B2 of Behroozi et al. (2019b). This procedure gives us a catalog of disrupted subhalos, also known as orphan subhalos, in addition to all the subhalos that can still be identified and tracked in the standard manner.
A problem remains that many of the disrupted subhalos that we continue to track are actually disrupted in a physical manner, for example through major mergers. We must decide which of these orphan subhalos may actually still host galaxies. There is a growing literature that attempts to address this problem as a function of properties of the subhalo itself and its orbital parameters (Ogiya et al. 2019;Jiang et al. 2021). In this work, we parameterize the probability that an orphan subhalo physically disrupts, i.e. that it is not available to be populated with a galaxy, as a function of the maximum velocity of the subhalo at its peak mass, v Mpeak : where Θ is the Heavyside step-function, and where T disr,low , and T disr,high are the asymptotes of T disr (v Mpeak ) at low and high v Mpeak respectively, σ disr determines the gradient of T disr (v Mpeak ) as a function of v Mpeak and v disr,mean sets the v Mpeak at which the slope of T disr (v Mpeak ) is greatest. We also investigated a number of other forms for T disr , including dependence on host v Mpeak and a joint dependence on host and subhalo v Mpeak , none of which improved our ability to fit all of the stellar mass bins in our data simultaneously, as discussed in section 5.2. It is also possible that baryonic effects disrupt subhalos as a function central galaxy morphology (e.g. the presence of a disk), star formation rate or super-massive black hole activity. In this work we have assumed that all subhalos that are resolved at z = 0 in our simulation can potentially host galaxies, and thus baryonic effects cannot disrupt resolved subhalos. This seems a reasonable assumption given the stellar mass ranges considered here, especially as there is only a very small population of halos that are significantly stripped but still resolved in our simulations.

Velocity Bias
As an additional extension of our model we consider the possibility that the velocity distributions of central and satellite galaxies differ from the velocity distributions of the (sub)halos in our simulations. We introduce two parameters to govern these potential deviations, α c and α s , which determine central and satellite velocity bias respectively.
We include central velocity bias by assigning central galaxies LOS velocities, v c , with respect to the LOS host halo center-of-mass velocity, v h , drawing from the dis-tribution: where σ c = α c σ h / √ 3, σ h is the three-dimensional host halo velocity dispersion as measured by Rockstar, and the factor of √ 3 comes from the conversion between three-dimensional velocity dispersion and LOS velocity.
Satellite galaxies, i.e. galaxies that have been assigned to subhalos, are treated separately. In the absence of satellite velocity bias, satellite galaxies are assigned the velocities of the subhalos that host them. Our implementation of satellite velocity bias simply re-scales the satellites velocities in the host-halo frame of reference: where v s is the satellite LOS velocity, and v sub is the LOS velocity of the subhalo that the satellite galaxy is assigned to.

Conditional Abundance Matching
In order to assign galaxy SSFRs to simulated galaxies, we adopt a conditional abundance matching (CAM) framework. CAM posits that at fixed stellar mass, SSFR is monotonically related to a second halo property, Y halo , i.e.
where F (SSFR|M * ) is the cumulative distribution of galaxy SSFR at fixed M * . No functional form is assumed for P (< SSFR|M * ) or P (< Y halo |M * ), rather, they are measured directly from the data and the simulations respectively. As a generalization of CAM, we allow for an imperfect correlation between Y halo and SSFR. This is accomplished by enforcing where This ensures that Y halo is a noisy version of Y halo , where the Pearson correlation coefficient between Rank(Y halo ) and Rank( Y halo ) is set to r. We use a "bin-free" version of CAM as implemented in halotools , where P (< Y halo |M * ) and P (< SSFR|M * ) are determined in sliding windows around the M * of each simulated halo and galaxy in the data respectively. This allows us to avoid the discreteness effects imparted by wide bins in M * evident in earlier implementations of CAM.
In a similar way to eq. (6), this defines a relation P (SSFR|Y halo , M * ). In this work we consider two different quantities for Y halo . The first, which we refer to as ∆v max , is a measure of recent accretion of matter onto subhalos defined as where a Mpeak is the scale factor at which the halo attains its peak mass, and a dyn is the scale factor one dynamical time before a, where a dynamical time is given by t dyn = 4 3 πGρ vir −1/2 . ∆v max is also used as a proxy for SFR in UniverseMachine (Behroozi et al. 2019a), albeit in a significantly more intricate manner.
The second proxy that we make use of is R h . R h is defined as the distance between a given galaxy and the closest host halo with mass greater than a specified threshold, M cut , where M cut is left as a free parameter that is fit to data. This proxy is motivated by claims that quenching has a strong dependence on proximity to massive halos, and that there is a particular mass scale for this quenching (Peng et al. 2010;Behroozi et al. 2013a;Zu & Mandelbaum 2016).

Measurements in Simulations
We take a slightly different approach to RSD measurements in our simulation than what we use on the SDSS data. In particular, the periodic nature of our simulations allows measurements to be made without the use of random points. This is important in order to improve the accuracy of the emulators we build in this work, which make use of small scale clustering measurements. Otherwise, a very large number of random points would be required to remove the contribution of their finite sampling from our measurements. Because of this difference between the measurements in our simulations, and those made in the data, the angular and redshift window functions of the data are not accounted for in our model predictions. Nevertheless, these are largely accounted for in our jackknife covariance matrix measured from the SDSS data.
For each stellar mass bin in table 1, we select an analogous sample in our simulations by cutting on the stellar masses that have been assigned via our SHAM model. We produce redshift-space coordinates for each galaxy from real-space positions by modifying the LOS coordinate as follows: where x z,rsd and x z,cos are the redshift-space and comoving line-of-sight coordinates and v z is the line-of-sight velocity in km s −1 . Periodic boundary conditions are applied after this transformation. We can then use the natural estimator: for ξ(r p , π) and ξ(s, µ), where DD and RR are again data-data and random-random pair counts normalized by the number of pairs in each bin. The RR term can now be expressed analytically, given our periodic boundary conditions. w p (r p ),ξ 0 andξ 2 are then calculated in the same way as described in section 2. We measure w p (r p ),ξ 0 andξ 2 three times, once using each of the x, y, and z axes of our simulation as the line-of-sight, using the average over the three lines of sight to build our emulators. For the purposes of estimating covariance matrices for our simulation measurements, which is necessary due to the non-negligible contribution of sample variance from our simulations to the total error in our analysis, we again use a jackknife procedure. Here we make use of 125 equal volume sub-boxes as our jackknife regions rather than the regions used for measurements on SDSS. Our covariance matrix is again estimated from our jackknife measurements using eq. (5). We compute this jackknife covariance matrix at the joint best fit parameter values for w p ,ξ 0 andξ 2 for the ∆v max CAM model discussed in section 6.1.
We must account for one additional source of uncertainty in our analysis: the stochasticity imparted on our measurements due to our use of Monte Carlo draws from random variables in our SHAM and CAM implementations. We do this by repopulating our simulations 10 times at each point in SHAM and CAM model space and recomputing w p ,ξ 0 andξ 2 for each re-population. The measurements presented in the following sections are the mean of these 10 re-populations. We neglect any residual contribution of this stochasticity in our error budget.

SHAM RESULTS
Instead of directly populating our simulations at each point in parameter space in order to make clustering predictions, we construct surrogate models for our SHAM and CAM models and use these to predict observables as a function of our model parameters. Appendix A describes the surrogate modeling framework we use to perform all of our parameter estimation in more detail.
All posterior parameter distributions and evidences in this work are derived using the nested sampling algorithm dynesty (Speagle 2019). We use a convergence criterion of ∆ log Z < 0.01 for all analyses, where Z is the evidence. All analyses assume a Gaussian likelihood, where the covariance used is the combination of the jackknife covariance matrix estimated from the data and the jackknife covariance matrix from the simulations used as an approximation for the error on our surrogate models. We assume flat priors on all parameters, whose edges are listed in table 2.

Baseline SHAM
We now study the behavior of our fiducial SHAM model including as parameters only the scatter in stellar mass at fixed SHAM proxy, σ log M * , and the parameter controlling the concentration dependence of the SHAM proxy, α.
The effect that these parameters have on projected and redshift-space clustering can be seen in fig. 2, where we show the effect of varying σ log M * and α (among other parameters to be discussed later) by 2σ around their best fit values for each stellar mass bin. It is clear that σ log M * has the greatest impact on the most massive bin, affecting w p andξ 0/2 at similar levels. This is expected, as the effect that variations in σ log M * have on clustering is highly dependent on the slope of the stellar mass function, and the bias-halo-mass relation b(M ). This is because increasing σ log M * preferentially brings galaxies hosted by lower mass halos into a given stellar mass selection. This effect is stronger the steeper the slope of the stellar mass function. As galaxies hosted by halos of lower mass scatter into the selection, the large scale bias of that sample is reduced because of the positive slope of b(M ). Thus, changes in σ log M * have the greatest effect where the slope of b(M ) and stellar mass function are simultaneously large, which occurs at high stellar masses.
Variations in α become quite important for the two less massive bins, where there is more scatter in v Mpeak − v vir relation (Lehmann et al. 2017). Increasing α ranks halos with larger concentrations higher at fixed M peak , thus preferentially selecting satellite galaxies, which generally form earlier and are thus more concentrated on average than host halos of the same M peak . As such, increasing α boosts the satellite fraction f sat of our samples, preferentially increasing both the one halo term of our clustering signals, which scales as the number of satellite galaxies squared, and the large scale bias of our samples, as satellites preferentially reside in more mas- sive and highly biased halos. Increasing α also preferentially selects more concentrated central halos at fixed M peak , again boosting the large scale bias due to the secondary dependence of halo bias on concentration, which is a positively sloped relation at this halo mass scale Mao et al. 2018). We will revisit the effects of these parameters on our clustering statistics when discussing model extensions in sections 5.2 and 5.3, but we note here that additional SHAM parameters such as artificial subhalo disruption and velocity bias severely hamper our ability to constrain α. Table 3 lists the reduced chi-squared values of the best fit models for a number of different data vector and galaxy sample combinations. For the two most massive galaxy samples (10.2 ≤ log 10 M * < 10.6 and 10.6 ≤ log 10 M * < 11.2, combined referred to as "Top two") our model obtains a good fit to w p ,ξ 0 andξ 2 simultaneously. All three clustering measurements, w p andξ 0/2 , show similar reduced chi-squared values, suggesting that there is not significant tension in the model when trying to fit all three statistics simultaneously. Nevertheless, we show in sections 5.2 and 5.3 that some model extensions can improve these fits.
We also see in table 3 that the goodness of fit for all clustering statistics is significantly worse for the 9.8 ≤ log 10 M * < 10.2 sample than for the more massive samples. This is particularly true for the redshift-space clustering statistics, and can also be seen in the relatively larger values for reduced chi-squared seen in the joint fit to all three samples. In the following sections, we will discuss further extensions to our fiducial SHAM model that improve our ability to fit this least massive galaxy sample, although no extension does a good job of simultaneously fitting all three samples. Figure 3 shows constraints on the fiducial SHAM model parameters when fit to statistics measured from the two most massive galaxy samples. We see that all three statistics prefer consistent values for the fiducial SHAM parameters. The parameter constraints from each statistic have very similar degeneracies, with σ log 10 M * and α showing little correlation. ξ 0 and ξ 2 are significantly more constraining than w p , mostly due to the relatively larger signal to noise of the multipole measurements.
We also show the best fit model to all clustering statistics in each stellar mass bin in fig. 1. This figure conveys the same impression as the goodness-of-fit statistics in table 3, namely that our model is a good fit to the two most massive galaxy samples and that there are more significant residuals in the least massive sample.
Our ability to fit the more massive samples runs contrary to the results presented in Guo et al. (2016), who used the same simulation to show that a SHAM model accounting for only scatter in absolute magnitude at fixed SHAM proxy cannot simultaneously fit projected and redshift-space clustering in the SDSS MGS. This discrepancy can be explained by considering the constraints that we obtain on α, which governs the behavior of the mass proxy used in our abundance matching models. In fig. 1 we see that our data rule out α = 0, which corresponds to M peak abundance matching, at high significance. M peak is approximately the same as the M acc abundance matching model presented in Guo et al. (2016), and so we see that when using a similar SHAM proxy we obtain similar results. Our constraints rule out α = 1, corresponding to v Mpeak at 1.4σ confidence. v Mpeak abundance matching is similar, but not identical, to the v acc and v peak abundance matching models that Guo et al. (2016) also rule out, although Guo et al. (2016) rules them out at much higher significance. This discrepancy may be a result of slight differences between v Mpeak , v acc and v peak , and may also be a result of different fitting procedures and galaxy samples. Guo et al. (2016) also use a different method to correct for fiber collisions than the one used here, allowing them to use smaller scale clustering measurements, which could also account for some of the differences between the findings presented here and those presented in their work.
Nevertheless, we see that allowing for a continuous degree of freedom governing the concentration dependence of our abundance matching proxy is important for obtaining good fits to all clustering statistics simultaneously. We have tested that this still holds when selecting galaxies using absolute magnitude in appendix B, as was done in Guo et al. (2016). Although we obtain improved fits with respect to those presented in Guo et al. (2016), we shall show in sections 5.2 and 5.3 that extensions to the baseline SHAM model presented here are still preferred by the data in some cases.

Orphan galaxies
We now examine physically motivated extensions to our fiducial SHAM model to see how they impact our parameter constraints and goodness-of-fit statistics. We consider two different model extensions: orphan galaxies and velocity bias.
First, we consider a SHAM extension that corrects for potential artificial subhalo disruption in SMDPL. We will often refer to this extension as an "orphan" model. Artificial subhalo disruption has received attention in the recent literature, with van den Bosch & Ogiya (2018) pointing out that artificial disruption is commonplace in cosmological N -body simulations at virtually all resolutions, even for subhalos that are resolved with hundreds of thousands of particles.
The model that we use to correct for this potential simulation artifact is presented in section 4.2. Orphan subhalos are tracked after disruption until their current maximum circular velocity, v now , falls below a fraction, T disr , of their v Mpeak value. The four parameters of the orphan model govern the v Mpeak dependence of T disr , although only one of them, the parameter controlling the low v Mpeak asymptote of T disr , T disr,low , has a large ef-fect on the SHAM model predictions. As we decrease T disr,low , we bring more subhalos into our sample, thus increasing f sat , again boosting the one and two-halo terms of our clustering signals. Changes in T disr,low almost exlusively change f sat , and thus give the dependence of our clustering signals slightly different scale and environmental dependence than changes in α and σ log M * . This can be seen in fig. 2, where T disr,low , is largely degenerate with α and σ log 10 M * on large scales, but on small scales has a different scale dependence especially inξ 0,2 , where the effect of f sat has an additional effect in boosting the amplitude of the finger-of-god, allowing for constraints on subhalo disruption that are not entirely degenerate with α and σ log M * .
In order to perform model comparisons, we make use of Bayes factors: i.e. the ratio of the Bayesian evidence for the data given an extended model, P (d|M ext ) to that obtained using our fiducial model, P (d|M fid ). If the ratio is larger than unity, then the extended model is preferred over the fiducial model, and if it is less than unity then the data prefer the fiducial model. In table 4, we list the Bayes factors for each of our extended models computed when fitting the clustering measurements of each stellar mass bin, the two most massive mass bins, and all mass bins simultaneously. For the most massive bin we find R < 1, but for the two less massive bins the SDSS data do prefer the presence of orphan galaxies. It is not surprising that the less massive bins considered here require orphan galaxies while the most massive bin does not. The effect of artificial subhalo disruption has been shown to have a halo mass dependence at fixed simulation mass resolution, with less massive and more poorly resolved subhalos disrupting sooner than their counterparts that are resolved with more particles (Ogiya et al. 2019).
The larger impact of orphans on less massive samples is also observable in the effect that including orphan galaxies has on our σ log 10 M * and α constraints, shown in fig. 4. The constraints on these parameters for the most massive galaxy sample are significantly less affected by including orphans in our model than the constraints from the less massive bins. For the least massive sample both α and σ log 10 M * become almost entirely unconstrained once orphans are included, and for the second most massive bin, α becomes largely unconstrained.
The changes in the constraints on α for the less massive bins are particularly interesting. The fiducial model favors α consistent with 1 for these bins, while including   Figure 1. Comparison of the best fit projected correlation function, wp(rp) , monopole ξ0(s) and quadrupole ξ2(s) models (lines) models the SDSS measurements (points) in three bins of stellar mass, as listed in the bottom row. Rows alternate between clustering measurements, and fractional residuals of our model from the measured data. Gray bars in the fractional residual panels represent the 1σ errors. Error bars on the points include both sample variance from the data and from our simulation. Shaded regions around the solid lines represent the 1 − σ posteriors of our model. The blue lines are the best fit model including artificial subhalo disruption and velocity bias as discussed in sections 5.2 and 5.3, while the orange line is our fiducial model, discussed in section 5. Velocity bias is important in order to fit the most massive sample, and orphan galaxies are preferred by the two lower mass samples. The solid points with error bars are the best fit model predictions, while the colored lines are the model predictions varying each parameter one by one. We see that the most important parameter for each stellar mass bin is different, with σ log M * , α, and T disr,low causing the largest variations in the most, second most, and least massive stellar mass bins respectively. Satellite velocity bias is important forξ 0/2 for all masses, having a small residual impact on wp because we do not project over infinitely long distances along the LOS. α, and T disr,low have similar impacts on large scales, but differ for r < 2h −1 Mpc , allowing some ability to constrain α in the presence of subhalo disruption systematics.  Figure 3. Constraints on the parameters of the fiducial SHAM model, σ log M * , log-normal scatter in M * at fixed vα, and α, the parameter that controls the concentration dependence of the SHAM proxy, fit to the two most massive galaxy samples. Various contours are constraints from different clustering measurements, all of which are statistically consistent. The combined result rules out M peak (α = 0) SHAM at many sigma, and v M peak SHAM (α = 1)at 1.4 sigma.
orphan galaxies removes nearly all constraining power on α suggesting that previous results favoring v Mpeak abundance matching over M peak were largely a result of those models' lack of orphan galaxies, similar to the claims made in Campbell et al. (2018). Finally, we examine the v Mpeak dependence of our orphan model for these three different stellar mass selections in fig. 5. For the most massive bin, we again see the preference for no orphan galaxies in the fact that T disr ≥ 1 for all v Mpeak , meaning that the data do not require the inclusion of any orphans for any values of v Mpeak populated by our model. The best fit in the second bin shows some preference for orphan galaxies, but the no-orphan scenario is not ruled out. On the other hand, for the least massive sample, we see that the mean constraint on T disr is less than one for all values of v Mpeak , indicating that the SDSS data prefer orphan galaxies.
We also see that there is some tension between the orphan model constraints from this least massive sample and the most massive sample, with the least massive sample requiring orphan galaxies in v Mpeak ranges where the most massive sample rules out orphans. The bottom panel of fig. 5 shows the v Mpeak distributions for the three mass bins. We see that all three have significant overlap in v Mpeak and so the different preferences for T disr at fixed v Mpeak lead to an inability to simultaneously model the least and most massive bins. This Constraints on σ log 10 M * for all three stellar mass bins individually, the two most massive bins analyzed simultaneously ("Top two") and all three bins analyzed simultaneously ("All"). Points are the 1 dimensional posterior means, and error bars depict 95% confidence intervals. Different colored points represent different model choices, including our fiducial two parameter SHAM model (blue), an extension that allows for orphan galaxies (orange) and an extension allowing for both orphans and velocity bias (green). The constraints on σ log 10 M * are not significantly affected by varying these modeling choices, except for the least massive bin where inclusion of orphans significantly degrades the ability to constrain σ log 10 M * . (Left) Same as right panel, but showing constraints on α. Allowing for orphan galaxies significantly degrades the constraints in α for the two less massive galaxy samples. This indicates that previous preferences for v Mpeak abundance matching over M peak abundance matching in the literature may have been driven by artificial subhalo disruption.
tension is also borne out by the very small Bayes factors for this model extension when considering all mass bins. This finding suggests that there must be a secondary variable that controls artificial subhalo disruption that our model does not account for. We have tried modifying our disruption model to account for more complex dependencies, including dependence on host halo v Mpeak and concentration, none of which alleviate this tension. We leave more extensive investigation of these additional dependencies to future work.

Velocity Bias
We have also explored an extension to our fiducial SHAM model that allows galaxies to have different velocity distributions from the subhalo populations that they are assigned to. Such an extension is commonly referred to as velocity bias. We note that the constraints on velocity bias presented in this work are not directly comparable to constraints obtained from HOD models that place galaxies on particles or model satellite veloc- Constraints with T disr ≥ 1 imply no need for orphan galaxies. (Bottom) v Mpeak distribution preferred by the best fit models for each sample. The significant overlap in these distributions for the least and most massive samples, along with the different inferred T disr for these two samples suggests that additional complexity is required in our orphan model in order to account for subhalo disruption consistently across the whole stellar mass range considered in this work.
ity distributions assuming isotropic Jeans equilibrium in an NFW profile, as the unbiased velocity distributions used in this work come from subhalo populations. Here we consider two additional parameters. The first, α c , allows central galaxies to have non-zero velocity with respect to the center-of-mass of the halo they are hosted by. The second parameter re-scales the velocity dispersions of satellites with respect to their host halo by a multiplicative factor, α s . The implementations of these models are described in section 4.3.
In table 4 we see R > 1 for model extensions that include only velocity bias other than for combinations that include the middle mass bin. Figure 6 shows constraints on the velocity bias parameters for each stellar mass selection. The blue points are for model extensions including only velocity bias, and orange points show constraints when including orphan galaxies as well as velocity bias. The inclusion of orphan galaxies does not appreciably change our velocity bias constraints. It is apparent that the preference for the velocity bias model in the most and least massive galaxy samples is driven by a deviation of α s from unity at a significance of 3.1σ and 2.6σ for the least and most massive bins respectively. When analyzing the two most massive samples simultaneously with one set of velocity bias parameters, this preference for non-unity satellite velocity bias is significantly decreased. This may indicate that the detec- tion of velocity bias in the most massive galaxy sample is due to noise, or it may simply indicate that the 10.2 ≤ log 10 M * < 10.6 bin is driving the fit back towards α s = 1. When all three bins are analyzed together, α s is consistent with unity, but we caution that this result may not be robust due to the poor overall fit to the combination of all three galaxy samples. The central velocity bias parameter, α c is consistent with 0 at 95% confidence for all samples when analysed individually as shown by the blue points in in the righthand panel of fig. 6. The constraints when analyzing all samples simultaneously prefer non-zero central velocity bias at greater than two sigma but once orphan galaxies are included we find α c consistent with zero for all samples.
As can be seen in fig. 4, including velocity bias has little effect on parameter constraints for σ log 10 M * and α. The slight exception to this is for the σ log 10 M * from the top two most massive samples, where including velocity bias changes the constraint by approximately 1σ.

CONDITIONAL ABUNDANCE MATCHING RESULTS
Here we investigate whether conditional abundance matching (CAM), is capable of modeling clustering as a function of stellar mass and specific star formation rate (SSFR). Previous sections explored SHAM extensions that may be required to model stellar mass complete galaxy selections, but such samples are rarely gathered in large spectroscopic surveys where redshift-space clustering can be measured with high precision. If CAM were capable of modeling selections in stellar mass and SSFR, then it would hold potential for modeling the clustering signals of luminous red galaxy and emission line galaxy samples observed in upcoming surveys such as DESI, 4MOST, Euclid and PFS, whose selections will implicitly depend on both stellar mass and SSFR in addition to other variables.
In the following subsections we investigate two distinct conditional abundance matching models: one that ties SSFR to the rate of change in potential well depth of subhalos as traced by their maximum circular velocities, and one that matches SSFR with distance to the nearest halo above a specified mass threshold.
The following analyses are limited to individual stellar mass selections, as we cannot estimate a covariance matrix that is simultaneously invertible for all stellar mass and SSFR bins. As such, we cannot present CAM results using the combined constraining power of all stellar mass bins, as we did for the SHAM models. Instead, the following analyses are always based on models fit simultaneously to the quenched, star forming and combined samples for each individual stellar mass bin. As we shall see, the CAM models for individual mass bins are not entirely consistent, and so analysis of all bins in combination would have limited utility even if we could obtain a reliable covariance matrix for such an analysis.

∆v max CAM
The first CAM model that we explore is one that ties SSFR to ∆v max (see eq. (14)). This model is motivated by the success of the UniverseMachine in fitting the clustering and number densities of stellar mass and SSFR selected galaxies as a function of time.
UniverseMachine also ties galaxy SSFR to ∆v max albeit in a model with many more free parameters in order to model star formation over a much broader range in time. Here we focus on a simplified version of the UniverseMachine model, where we allow for a single free parameter: the linear correlation coefficient between SSFR and ∆v max at fixed galaxy stellar mass. The implementation of this model is discussed in section 4.4.
The solid lines in fig. 7 show the fits of the ∆v max CAM model to w p , andξ 0/2 for quenched (red), starforming (blue), and combined (black) samples for each stellar mass bin. Similar to fig. 1, the fits shown in each stellar mass bin are performed independently.
The correlations that allow the ∆v max CAM model to fit these signals are similar to those at play in the baseline SHAM model. For ∆v max CAM with r ∆vmax = 1, galaxies with lower values of ∆v max at fixed M * are monotonically assigned galaxies with lower SSFRs. Subhalos preferentially have lower values of ∆v max at fixed halo mass, and so subhalos are preferentially assigned lower values of SSFR, thus boosting f sat for more quenched samples. Additionally, ∆v max CAM assigns low SSFR values to the oldest host halos that are no longer accreting mass, thus boosting the central assembly bias signal in a similar way to increasing α in SHAM (Zentner et al. 2014). As r ∆vmax is decreased, the correlation between SSFR and f sat , assembly bias, and, by proxy, clustering amplitude becomes flatter. Thus, r ∆vmax governs the ratio of clustering amplitudes of the various SSFR selections at fixed M * .
By eye, we see that the fits to the quenched and star forming samples are reasonable, but slightly worse in general than the fits to the combined samples. This is shown more quantitatively by the reduced chi-squared values shown in table 5. The reduced chi-squared values from the w p +ξ 0/2 row of table 3 are included for reference and listed as "Combined SHAM" indicating that the SHAM model is fit to the "combined" sample of galaxies binned in stellar mass but not SSFR. Fitting to the quenched and star-forming samples simultaneously with the combined sample degrades the fit to the combined sample for all stellar mass bins, indicating that adjustments from the best fit model from section 5 are required in order to simultaneously fit the quenched and star-forming galaxy samples. The most striking example of this is a large shift in the orphan model parameters for the least massive sample, where the ∆v max CAM model prefers essentially all artificially disrupted subhalos to be kept as orphans in order to fit the very large one halo term for the quenched sample in this stellar mass bin. Evidently, ∆v max CAM on its own cannot boost the one-halo clustering signal of this least massive bin enough on it's own. Thus, more orphan galaxies must be included in order to further boost f sat , and these orphans must be preferentially quenched. This happens naturally in ∆v max CAM, as orphans will be preferentially stripped and thus assigned lower SSFRs. In general, the constraints on the SHAM parameters are also broadened when fit simultaneously with r ∆vmax , again indicating a slight tension in the CAM model when fit to the data presented here. Figure 8 shows constraints on r ∆vmax in blue for the three stellar mass bins considered here. Notably, there is little evolution of this parameter as a function of stellar mass, with the least massive bin yielding constraints of . Comparison of the best fit ∆vmax (solid) and R h (dashed) CAM models to wp andξ 0/2 measured for quenched (red), star-forming (blue) and combined (black) samples for each stellar mass bin. Shaded regions around the solid lines are 1σ posterior model predictions. 1σ posterior distributions are not included around the R h models for clarity as they are comparable to the ∆vmax distributions in size. Overall the fits for both models are comparable, although the ∆vmax model generally outperforms the R h model for star-forming samples. Both models struggle to fit the quenched galaxy measurements for the two less massive bins on scales r ≤ 1 h −1 Mpc .  Figure 8. Posteriors of the linear correlation coefficient, rSSF R, between SSFR and the two CAM proxies, ∆vmax (blue), and R h (yellow) for each stellar mass bin. Except in the least massive bin, ∆vmax is more correlated with SSFR than R h , again indicating that ∆vmax is a better proxy for SSFR. We also see that rSSF R varies much more strongly with stellar mass for the R h model than the ∆vmax model. r ∆vmax = 0.58 ± 0.01, while the two most massive samples prefer r ∆vmax = 0.52±0.02 and r ∆vmax = 0.54±0.13 respectively.

R h CAM
The second CAM model we consider is one that ties SSFR to R h , the distance to the nearest host halo with mass above a threshold M cut where M cut is left as a free parameter. Here, the physics that creates the differences between high and low R h samples is quite different than in the ∆v max CAM case. Subhalos that are close to massive host halos and thus have low values for R h are assigned lower values of SSFR. This again acts to boost f sat at fixed M * for quenched galaxy samples, but here there is no direct correlation between halo age and quenching, thus this model will not impart central galaxy assembly bias signals in the way that the ∆v max CAM model does. On the other hand, this model will introduce large correlations between quenching of nearby galaxies, also known as galactic conformity (Weinmann et al. 2006), which can have a large effect on one and two halo clustering amplitudes . Like the ∆v max CAM model, the R h CAM model also allows the correlation coefficient between R h and SSFR, r R h , to be free, again controlling the ratio of clustering amplitudes for quenched and star-forming samples.
The dashed lines in fig. 7 show the best fit R h model to each stellar mass bin. Overall the fits are comparable to the ∆v max model, although with noticeably worse performance for star-forming samples. The reduced chisquared values for these fits are displayed in table 5, where we indeed see that the R h model is a worse fit to the star-forming samples. Table 6 shows Bayes factors comparing the ∆v max and R h model. The ∆v max model is preferred by the data for the two most massive bins, although the significance of this preference is by far the greatest in the second most massive galaxy sample. The R h model is preferred in the least massive sample, but neither model is a particularly good fit to the data for this bin. It should be noted that these Bayes factors depend on the prior that we have chosen for M cut , and a significantly smaller prior on this parameter would result in smaller Bayes factors and less of a preference for the ∆v max model. Joint constraints on the R h CAM model parameters are shown in fig. 9, where a strong correlation between r R h and M cut can be seen for all stellar mass bins. Unlike with the ∆v max model, there is a strong trend in the correlation between R h and SSFR as a function of stellar mass. There is also an apparent trend in M cut with stellar mass. More massive samples tend to be quenched by their proximity to more massive halos, as indicated by the larger preferred M cut values as a function of stellar mass. The 1-dimensional posteriors on r R h are compared with those obtained from the ∆v max CAM model in fig. 8. Table 6. CAM Bayes Factors 9.8 ≤ log 10 M * < 10.2 10.2 ≤ log 10 M * < 10.6 10.6 ≤ log 10 M * < 11.2 5.8 × 10 −6 ± 1.7 × 10 −6 9.9 × 10 7 ± 3.1 × 10 7 21 ± 5.2

Model predictions for finer color selections
Given the comparable ability of both CAM models to fit the clustering statistics considered thus far, we now look more closely at the clustering of samples selected more finely in SSFR-M * space. To do so, we divide the SDSS galaxies into four equal quartiles in SSFR, independent of M * , yielding four bins in SSFR: with Y i ∈ {−99, −11.84, −10.98, −10.072, −6.39}. It should be noted that the four bins in SSFR do not contain equal numbers of galaxies for each stellar mass bin because the SSFR bin edges are determined using all galaxies with 9.8 ≤ log 10 M * < 11.2. We use the same M * and redshift bins as those discussed in section 2. Errors on the data are again estimated via jackknife, and the simulation errors are estimated using the jackknife procedure discussed in section 3. Instead of fitting our CAM models to these new finer SSFR selections, we make predictions for these selections given the MAP CAM parameters presented in the previous section using the joint fit to w p ,ξ 0 andξ 2 for each stellar mass bin. These predictions are shown in fig. 10, where the points represent the SDSS measurements, and the solid lines and contours are the best fit and 1σ posterior of the ∆v max CAM model, while dashed lines are the MAP R h model.
Given the small differences in the performance of these models when fit to the quenched and star-forming samples, the marked difference in the predictions for this finer binning in SSFR is surprising. For the least massive sample, the ∆v max model significantly out-performs the R h model, with the R h model predicting a very small difference in clustering amplitude for the lower SSFR samples, and a very large difference in clustering amplitude for the higher SSFR samples. These predictions average out to be roughly correct when performing the binary splits in SSFR, but these finer splits make it clear that the reasonable fit to the binary split samples is merely coincidental. For the two most massive bins, the predictions for the ∆v max and R h model are more comparable, with the ∆v max model outperforming the R h model for star-forming samples, while the R h model is a better fit to the third lowest SSFR bin.
In the two less massive bins, the large difference in clustering amplitude for the two lower SSFR bins is in-teresting in its own right, suggesting that there is a diversity in the quenched galaxy population that is explained to a large extent by diversity in ∆v max . This can be compared to the two most star forming samples, whose clustering amplitudes are much more comparable, suggesting less diversity of environments when galaxies are on the star-forming main sequence than for quenched galaxies. This further suggests that that the data may prefer a model with separate correlation coefficients between SSFR and ∆v max for quenched and star-forming galaxies, with ∆v max being more correlated with SSFR in the quenched regime than the star-forming regime.

Environmental Dependence of Galaxy Quenching
Because galaxy auto-correlation functions up-weight regions of high galaxy density, they are largely sensitive to the clustering of satellite galaxies, which preferentially live in these regions. For this reason, the clustering measurements presented thus far are not particularly sensitive to the processes governing central galaxy quenching, and thus not sensitive to one of the main pitfalls of the original formulation of CAM: its inclusion of very strong environmental dependence on quenching of low-mass isolated galaxies. Tinker et al. (2017) showed that for low mass galaxies the quenched fraction, f q = n quenched n all , predicted by age matching depends very strongly on local galaxy overdensity, while in SDSS this dependence is much weaker. Here we perform our own comparison of similar measurements to the updated CAM models presented in this work to see if the reported tension persists. We adopt slightly different definitions than those presented in Tinker et al. (2017) because we use SSFR rather than D 4000 , and because we have not run a group finder on our simulations in order to identify isolated galaxies. We define the quenched fraction as the ratio of the number of galaxies with log 10 SSFR ≤ −11 to the total number of galaxies. We measure galaxy density, ρ, as the number of galaxies in an annulus in redshift-space defined by 0.5 h −1 Mpc <= r p < 4 h −1 Mpc and δcz < 1000 km/s, where we only count galaxies that are less massive than the galaxy in question, and more massive massive than 0.3M * , where M * is the mass of the galaxy around which we are counting satellites.
The top row of fig. 11 shows a comparison of f q as a function of galaxy density between our best fit CAM models from section 4.4 for all galaxies in each stellar  Figure 10. Predictions of the CAM models that provide the best fit to the combination of wp,ξ0 andξ2 measurements using a finer SSFR selection. Lower SSFR selections are represented by redder points and lines, while higher SSFR selections are represented by bluer points and lines. The highest SSFR sample for the most massive bin is not shown as it is extremely noisy. We see that the ∆vmax model outperforms the R h model for most samples, especially for the most star-forming and most quenched samples in the two less massive bins. mass bin. We have limited the SDSS sample's redshift range to z < 0.064 to ensure volume completeness for all samples simultaneously, and we perform the CAM procedure using the color distributions from this sample. We see that agreement between the best fit ∆v max model and the data is very good for all stellar masses, while the agreement between the best fit R h model and the data is significantly worse, with the R h model predicting a different shape of f q as a function of galaxy density, especially for low densities.
The bottom row of fig. 11 shows similar measurements, but this time only considering primary galaxies, defined as those galaxies that do not have a more massive galaxy within a region of r p < 0.5 h −1 Mpc and ∆cz < 1000km s −1 . This is the same definition used in Behroozi et al. (2019b). We see that the trends in f q with ρ become less steep than their counterparts measured for all galaxies for both the SDSS data and the CAM models. Agreement between SDSS and the ∆v max model is still very good for all samples, while the R h model shows a similar shape mismatch as seen when considering all galaxies. Thus the tension reported in Tinker et al. (2017) is alleviated for the ∆v max model. In order to demonstrate the aspect of the updated CAM model that alleviates this tension, we also show ∆v max model predictions for a model with r ∆vmax = 1 as dashed lines. This r ∆vmax = 1 model is closer to what was used in Tinker et al. (2017). We see that there is a much stronger trend in f q as a function of galaxy density when assuming this perfect correlation, which is very similar to that reported in Tinker et al. (2017) and in more sig-nificant tension with the SDSS measurements presented here than the best fit ∆v max model.
Recently O'Donnell et al. (2021) reported that measurements of number density profiles around isolated Milky-way mass galaxies in SDSS suggest that, if anything, star-formation rates are anti-correlated with dark matter accretion rates, seemingly contradicting the findings presented here. There are a number of differences between the measurements presented here and in O'Donnell et al. (2021), specifically the isolation criteria used, which are more conservative in O'Donnell et al. (2021), and the use of mean number profiles as a function of scale, rather than the counts-in-cells-like measurements that we use here. O'Donnell et al. (2021) also use accretion rates rather than ∆v max . Due to these differences it is unclear whether O'Donnell et al. (2021) contradicts our results, and further investigation confronting the statistics presented in their work with the models presented here is warranted.
Another apparent tension between CAM and SDSS data was discovered in the joint analysis of galaxygalaxy lensing and projected clustering presented in Zu & Mandelbaum (2016). They showed that at fixed stellar mass, galaxy quenching mostly depends on halo mass, with galaxies hosted by more massive halos exhibiting stronger quenching, whereas the age-matching model discussed in their work exhibits the opposite trend, i.e. at fixed stellar mass, central galaxy quenching is inversely related to halo mass. More massive galaxies accumulate their mass at late times, and thus are younger and assigned bluer galaxies. We have investigated this tension with the best fit ∆v max model, and found that unlike in age matching with a perfect correlation between z starve and SSFR, our best fit ∆v max model shows no signs of this inverted quenching relationship.

CONCLUSION AND FUTURE WORK
In this work we have explored the ability of abundance matching models to fit redshift-space clustering data measured from the Sloan Digital Sky Survey Main Galaxy Sample. The parameter constraints in this work make use of polynomial chaos expansions as surrogates for abundance matching models, facilitating many aspects of our analysis. Our main findings are as follows: 1. SHAM models including free parameters for scatter and the concentration dependence of the SHAM proxy can fit redshift-space clustering measurements.
2. Orphan galaxies improve our SHAM fits, reducing the preference for v Mpeak abundance matching over M peak .
3. We obtain > 2σ detections of satellite velocity bias, emphasizing the need for this degree of freedom to be included when fitting RSD statistics using SHAM.
4. Using ∆v max as a proxy for SSFR in CAM models provides a good fit to the clustering data presented in this paper.
In more detail, we fit a fiducial SHAM model with free parameters governing the scatter in stellar mass at fixed SHAM proxy, σ log 10 M * , and the concentration dependence of the SHAM proxy, α, to SDSS RSD measurements in three bins of stellar mass. We find that this model is able to reproduce these measurements with high fidelity except for in the least massive galaxy sample.
Section 5.2 and section 5.3 explore extensions to this fiducial SHAM model, allowing for orphan galaxies and velocity bias. Orphan galaxies are slightly preferred by the data for the two less massive galaxy samples, and satellite velocity bias is preferred by the least and most massive galaxy samples, as quantified by Bayes factors listed in table 4. Due to a lack of a consistent trend in preference or exclusion of these models, we suggest that they be included in future analyses using SHAM. Furthermore, including orphan galaxies degrades the ability to constrain α, suggesting that previous SHAM constraints that prefer v Mpeak over M peak may have been driven by artificial subhalo disruption.
While orphan galaxies allow SHAM to fit the least massive samples investigated in this work, the orphan model constraints from the three different samples are in tension. The least massive sample prefers a large number of orphans, boosting halo occupations of galaxies by a factor of ∼ 1.5 across a broad halo mass range with respect to a no-orphan model, while the most massive sample prefers no orphans. This suggests that there are additional variables governing artificial subhalo disruption that are not included in our model, but are required in order to achieve self consistency over the full stellar mass range considered in this work.
Having investigated stellar mass complete galaxy samples, we turn to samples selected in stellar mass and SSFR, as many samples used for cosmology in ongoing and upcoming surveys will be. We present two new conditional abundance matching models, one that uses a proxy for subhalo's matter accretion rate, ∆v max , and another that uses distance to massive halos, R h , as proxies for SSFR at fixed stellar mass. We have allowed for non-unity correlation between these proxies and SSFR in the form of a linear correlation coefficient. We find that both models are able to fit RSD measurements split 10 −1 10 0 10.6 < log M * < 11.2 Figure 11. Predictions of quenched fraction as a function of galaxy density compared to measurements in SDSS for the ∆vmax (blue) and R h (orange) CAM models fit to wp andξ 0/2Ṫ he top row shows measurements of the quenched fraction of all galaxies, while the bottom row considers only primary galaxies as defined in section 6.4. The different columns are different stellar mass bins, as labeled in the figure. The ∆vmax model is in good agreement with the SDSS data for all magnitudes and densities, whereas the R h model shows a significantly stronger correlation between galaxy density and fq than the data both for all galaxies, and the primary only measurements. The dashed blue line shows a ∆vmax model with r∆v max = 1, similar to that used in Tinker et al. (2017), demonstrating that the inclusion of non-unity correlation between ∆vmax and SSFR allows the data presented in this figure to be fit significantly better than the perfect correlation case.
into quenched and star-forming samples reasonably well, but the ∆v max model generally performs better than the R h model. While promising, neither CAM model is able to fit the quenched and star-forming samples as well as our SHAM model is able to fit the stellar mass complete samples.
In order to perform further comparison between the CAM models, we have made predictions for additional statistics at the models' maximum a posteriori parameter values when constrained against w p ,ξ 0 andξ 2 . First, we examined the dependence of w p andξ 0/2 in four bins of SSFR for each stellar mass bin, finding that the ∆v max model predictions match the SDSS measurements sig-nificantly better than the R h model, particularly for the higher SSFR samples.
We have also examined the a posteriori predictions of our CAM models for quenched fraction as a function of galaxy overdensity. Here, the ∆v max model prediction is nearly perfect and is again significantly better than the R h model prediction, especially at low mass and low density when considering both primary and non-primary galaxies. If we force the correlation between SSFR and ∆v max to be perfect, then we recover the tensions between CAM and SDSS data reported in previous work, such as Tinker et al. (2017), suggesting that allowing for non-unity correlation between SSFR and ∆v max is what allows us to fit this data well.
This work suggests that SHAM is a promising alternative to other simulation-based models for galaxy redshift-space clustering, particularly for higher stellar masses where only two free parameters are required to fit the SDSS data considered in this work. Modeling orphan galaxies is important for the less massive samples, but significant additional investigation is required in order to obtain an orphan model that is flexible enough to model the whole stellar mass range considered here. Including non-unity correlations in CAM alleviates many of the tensions with SDSS data that have previously been reported. CAM provides a reasonable fit to the the RSD measurements from stellar mass and SSFR selected galaxy samples, but further work is necessary in order to improve these fits to the level required for robust cosmological inference. The development of SHAM and CAM models lay the bedrock for the construction of realistic mock galaxy catalogs. With more concerted effort, SHAM and CAM may provide a robust forward model, usable for constraining the growth of structure in our universe from a broad range of galaxy clustering statistics. the Partnership for Advanced Supercomputing in Europe (PRACE, www.prace-ri.eu) for funding the Mul-tiDark simulation project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, www.lrz.de).
The Bolshoi simulations have been performed within the Bolshoi project of the University of California High-Performance AstroComputing Center (UC-HiPACC) and were run at the NASA Ames Research Center. We are grateful to Stuart Marshall and the rest of the SLAC computing team for extensive support of this work.

A. SURROGATE MODELING
In order to facilitate our analyses, we build surrogates for the SHAM and CAM models considered in this work. In particular, we have a set of measurements, Y = {y 1 , ..., y N } ∈ R, at points x = {x 1 , ..., x N } ∈ R D over a domain D X , and some mapping f : X → Y . The mapping f is the SHAM/CAM model predictions for our clustering statistics when applied to SMDPL, and is slow to evaluate. We wish to construct a fast and accurate surrogate modelf for f . In order to do so we construct a Polynomial Chaos Expansion for f (Xiu 2010). If the points in X follow a joint probability density function p(x), and f is a finite variance model, i.e.
and f is drawn from a Hilbert space of finite variance functions, H, endowed with the inner product then we can construct a PCE of f such that where α = {α 1 , ..., α d } is a d-dimensional index and Ψ α is an element of the multivariate orthonormal polynomial basis of H. In order to evaluate the expansion in an efficient manner it must be truncated in order to yield a finite sum, yielding an approximation to f : where A ⊂ N. In this work, we choose a truncation rule such that A = {α |α| < p}, i.e. we keep all polynomials with total order less than p, and we optimize p to compromise between accuracy and computational efficiency. For finite p, the PCE coefficients can be determined by optimizing a specified loss function. In our case, we take our loss function to be the mean squared error over all X, in which case the coefficients m α are determined by solving the normal equation.
We build independent PCEs for each scale, statistic, galaxy sample and model. In order to determine coefficients for each PCE we take X to be 1000 points in SHAM parameter space drawn from a Latin Hypercube. We have found that 1000 points is sufficient to achieve an accuracy consistent with the errors on the measurements in our simulations for all the models considered in this work, as discussed below.
We populate SMDPL at each point x i ∈ X and measure w p (r p ),ξ 0 (s) andξ 2 (s) as outlined in section 4.5. This produces 1000 training points y i ∈ Y for each scale, statistic, galaxy sample and model. Given this data, we determine unique sets of coefficients m β α for each PCE, where β now denumerates the scale, statistic, sample and model under consideration. In order to determine the generalization error of our PCEs to points outside of X we perform a cross-validation procedure, iteratively leaving out each point x i ∈ X, refitting m β α , and calculating the residual of the PCE at the removed point, y i −f (x i ). Figure 12 shows the residual distribution of the PCEs for all scales and statistics averaged over the three different stellar mass bins considered in this work. We only show the residuals for PCEs trained on the model with the greatest number of parameters, i.e. the R h CAM model with velocity bias and subhalo disruption discussed in section 6.2, as we have found the residuals for this model to be the largest of all models considered in this work. Residuals are quoted as fractions of the error on each measurement, where the errors are determined via the jackknifing procedure outlined in section 4.5. The mean error at each point, shown as the white point in each violin, is consistent with the error on the measurement for all but the two smallest radial bins in w p (r p ). We compare 2nd and 3rd order truncation schemes, and find that the second order scheme is sufficient to produce residuals consistent with the error on our training set, and so we use this scheme for all analyses in this paper.

B. ABSOLUTE MAGNITUDE AND COLOR SELECTED SAMPLES
In addition to the stellar mass and SSFR selected samples discussed in the rest of this work, we have also analyzed samples selected by z = 0.1-frame r-band absolute magnitude and g − r color. The M r selections that we use are listed in table 7, and we divide red and blue samples according to g − r > 0.21 − 0.03 M r .
We make use of the absolute magnitude function described in Wechsler et al. (2021) for our SHAM models, but otherwise the models are implemented identically to the descriptions in section 4. Figure 13 shows the fiducial and extended SHAM model fits to each absolute magnitude bin. We see similar trends in these fits to those discussed in section 5, with the two fainter  Figure 12. Test of accuracy of the emulator compared to sample variance of SMDPL. Each panel shows the distribution of emulator residual divided by estimated sample variance in SMDPL for each radial bin in wp,ξ0 andξ2 in the top, middle and bottom panels respectively. Distributions are constructed by performing leave-one-out tests where, each point in our Latin Hypercube is left out one at a time as described in section 4. The mean residual (white point in the violin) is smaller than sample variance (log 10 χ 2 ≤ 0), for all but the two smallest rp bins for wp showing that our emulator is performing at the level of sample variance in our simulations or better. The blue and green distributions represent the residuals for a 2nd and 3rd order PCE respectively. We see that the 2nd order expansion is sufficient to model this data and actually outperforms the 3rd order expansion.  bins better fit by the extended model, while the brighter bin is fit equally well by the fiducial model and the extended model, although Bayes factors do not prefer the extended models in any of the individual M r bins. The goodness of fit for the M r selected samples is overall slightly worse than for the stellar mass selected samples but qualitatively the fits are still quite reasonable. Fits of the ∆v max and R h CAM models to galaxy samples subdivided by M r and g − r color are shown in Figure 14. Again the results are comparable to those found for the M * and SSFR selected samples presented in section 6, with the ∆v max model slightly outperforming the R h model. The main differences that are apparent are both models' ability to fit the faintest red sample for the M r selections, while the CAM models struggled to fit  fig. 1, but fitting to samples selected by Mr instead of M * . Here, the models perform comparably to the stellar mass selected case although with slightly worse reduced chi-squared values for the fainter two samples. None of the model extensions are preferred by Bayes factors, but there is a significant improvement in chi-squared for the extended model in the second Mr bin. the analogous sample for the M * and SSFR selections. Unlike for the M * and SSFR selections, the model fits presented here cannot fit the smallest scales for the two fainter blue samples. We also find that the correlation coefficients between g − r and our CAM proxies are significantly larger than in the SSFR case, suggesting that these CAM models describe the g − r split data better than the SSFR selections.  fig. 7, but for Mr and g − r selections, where red lines represent red galaxy selections as defined in eq. (23), and blue lines are blue galaxy selections. Overall the fits for both models are comparable, although the ∆vmax model generally out-performs the R h model for both blue and red galaxy samples. Unlike for the stellar mass and SSFR selections, the ∆vmax model performs comparably across all Mr bins, although with some systematic deviations for r ≤ 1h −1 Mpc for blue galaxy samples.