The Stellar-to-halo Mass Ratios of Passive and Star-forming Galaxies at z ∼ 2–3 from the SMUVS Survey

The SMUVS program (PI: K. Caputi; Ashby et al. 2018) has collected ultra-deep Spitzer 3.6 and 4.5 μm data within the COSMOS⁶ field (Scoville et al. 2007) covering three of the UltraVISTA ultra-deep stripes (McCracken et al. 2012) with deep optical coverage from the Subaru Telescope (Taniguchi et al. 2007). The UltraVISTA data considered here correspond to the third data release,⁷ which reaches an average depth of K_s = 24.9 ± 0.1 and H = 25.1 ± 0.1 (2 arcsec diameter, 5σ).

A thorough description of the SMUVS multiwavelength source catalog construction and spectral energy distribution (SED) fitting is given in Deshmukh et al. (2018); however, we summarize the main details here.

Sources are extracted from the UltraVISTA HK_s average stack mosaics using SExtractor (Bertin & Arnouts 1996). The positions of these sources are then used as priors to perform iterative point-spread function fitting photometric measurements on the SMUVS 3.6 and 4.5 μm mosaics, using the daophot package (Stetson 1987).

For all of these sources, 2 arcsec diameter circular photometry on 26 broad, intermediate, and narrow bands (U to K_s) is measured (Deshmukh et al. 2018). After cleaning for Galactic stars using a B–J–[3.6] color selection (e.g., Caputi et al. 2011), and masking regions of contaminated light around the brightest sources, the final catalog contains ∼2.9 × 10⁵ UltraVISTA sources with a detection in at least one IRAC band over an area of ∼0.66 square degrees.

The SED fitting is performed with all 28 bands (26 U through K_s as well as Spitzer 3.6 and 4.5 μm) using the template-fitting code LePhare⁸ (Arnouts et al. 1999; Ilbert et al. 2006). We assume Bruzual & Charlot (2003) templates corresponding to a simple stellar population formed with a Chabrier (2003) stellar initial mass function and either solar or sub-solar (1 Z_⊙ or 0.2 Z_⊙) metallicity, and allow for the addition of nebular emission lines. Additionally, we assume exponentially declining star formation histories.

Photometric redshifts and stellar mass estimates are obtained for >99% of our sources. Using ancillary spectroscopic data in COSMOS to assess the quality of the obtained photometric redshifts, we found that the standard deviation, σ_z, of $| {z}_{\mathrm{phot}}-{z}_{\mathrm{spec}}| /(1+{z}_{\mathrm{spec}})$, based on ∼1.4 × 10⁴ galaxies with reliable spectroscopic redshifts in the COSMOS field (see Table 1 in Ilbert et al. 2013 and references therein) is 0.026 (Deshmukh et al. 2018) and is 0.035 for the z ∼ 2–3 population. This statistic is computed excluding outliers, defined as objects for which $| {z}_{\mathrm{phot}}-{z}_{\mathrm{spec}}| /(1+{z}_{\mathrm{spec}})\gt 0.15$, which comprise only ∼5.5% of the whole spectroscopic catalog. These results compare favorably with other photometric surveys in the literature (e.g., Ilbert et al. 2013; Laigle et al. 2016) and highlight the high accuracy of our derived photometric redshifts. For our parent sample of z ∼ 2–3 galaxies, we find σ_z = 0.035.

Throughout this paper, we use the best-fit redshifts and stellar masses computed by LePhare. We also limit ourselves to $\mathrm{log}({M}_{\star }/{M}_{\odot })\gt 9.8$ to ensure a high level of completeness for both our passive and star-forming galaxies. Our combined (i.e., passive + star-forming) galaxy population is 95% complete for $\mathrm{log}({M}_{\star }/{M}_{\odot })\gt 9.6$ at z ∼ 2–3 (Deshmukh et al. 2018). This stellar mass cut leaves us with ∼10³ passive and ∼6.8 × 10³ star-forming galaxies with 2 < z < 3.

Passive galaxies, i.e., those for which the current star formation is negligible relative to their past average, are defined using the criteria of Deshmukh et al. (2018). A galaxy is defined as passive if: (i) it has a color excess E(B − V) ≤ 0.1 (indicating that interstellar dust has little impact on the SED) as computed by LePhare; and (ii) it has a rest-frame color u − r ≥ 1.3, i.e., is "red". The use of the color cut in conjunction with the SED extinction prevents dust-obscured star formation from being the main cause of the red color. This leaves an old stellar population with negligible ongoing star formation as the most plausible explanation of the galaxy's SED.

We compute the galaxy stellar mass functions for our passive and star-forming populations using the 1/V_max technique (Schmidt 1968) as was done by Deshmukh et al. (2018). A V_max correction is computed for each galaxy based on the maximum volume at which it would have a magnitude [4.5 μm] = 26 (or [3.6 μm] = 26 in the case of a non-detection at 4.5 μm). For fainter sources, no V_max correction is applied. A further completeness correction is applied to each galaxy based on its 4.5 μm magnitude (or 3.6 μm in the case of a non-detection at 4.5 μm) and the ratio of SMUVS number counts to the completeness corrected number counts derived from the Spitzer-CANDELS/COSMOS survey (Ashby et al. 2015) at that magnitude.

The errors are derived by summing a Poisson noise term, σ_Poi, an error associated with photometric uncertainties, σ_mc, and a cosmic variance error, σ_cv, in quadrature. The Poisson term, σ_Poi, is derived using the tabulated values of Gehrels (1986). The photometric term, σ_mc, is derived by scattering the photometry input to LePhare within its uncertainties for each galaxy to create 100 realizations of the SMUVS catalog and recomputing the stellar mass function for each; σ_mc is then taken to be the 16–84th percentile of these mock catalog stellar mass functions in each stellar mass bin. The cosmic variance term is derived using the prescription of Moster et al. (2011).

Deshmukh et al. (2018, see their Figure 19) provides a thorough comparison of 2 < z < 6 SMUVS-derived galaxy stellar mass functions to earlier estimates in the literature.

The two-point angular autocorrelation function, w(θ), describes the excess probability, compared to a random (Poisson) distribution, of finding a pair of galaxies at some angular separation θ > 0. It is defined such that

$$\begin{eqnarray}&&{\delta }^{2}{P}_{12}={\bar{\eta }}^{2}\,[1+w(\theta )]\,\delta {{\rm{\Omega }}}_{1}\,\delta {{\rm{\Omega }}}_{2},\end{eqnarray} \tag{ 1 }$$

where $\bar{\eta }$ is the mean surface density of the population per unit solid angle, θ, is the angular separation, and δΩ_i is a solid angle element. Here, the angular correlation function is computed according to the standard Landy & Szalay (1993) estimator

$$\begin{eqnarray}&&w(\theta )=\displaystyle \frac{\mathrm{DD}(\theta )-2\mathrm{DR}(\theta )+\mathrm{RR}(\theta )}{\mathrm{RR}(\theta )},\end{eqnarray} \tag{ 2 }$$

where DD, DR, and RR represent the number of data–data, data–random, and random–random pairs in a bin of angular separation, respectively. The random catalog is constructed to have the same angular selection as the data by generating random points over the survey area and discarding those that fall within masked areas relating to, e.g., foreground stars. Due to the high stellar mass completeness of our samples (≳95%) we do not need to consider any variations in sensitivity/depth across our survey. We use a single random catalog with ∼5 × 10⁵ objects (∼50× greater than our largest galaxy sample) throughout our analyses. The errors on the two-point angular correlation function are estimated from the data using a jackknife method (Norberg et al. 2009). We divide the SMUVS footprint into 60 approximately equal area regions. Removing one region at a time, we compute the covariance matrix as

$$\begin{eqnarray}&&{C}_{i,j}=\displaystyle \frac{N-1}{N}\displaystyle \sum _{l=1}^{N}({w}_{i}^{l}-\bar{{w}_{i}})\times ({w}_{j}^{l}-\bar{{w}_{j}}),\end{eqnarray} \tag{ 3 }$$

where N is the total number of regions, $\bar{w}$ is the mean correlation function (${\sum }_{l}{w}^{l}/N$), and w^l is the estimate of w with the lth region removed (i.e., computed from 59 of the 60 regions). We compute w(θ) using eight evenly spaced logarithmic bins of angular separation in the range $-3\lt {\mathrm{log}}_{10}(\theta /\deg )\lt -1$.

From the HOD (Equation (11)), we can calculate the stellar mass function Φ(M_⋆,1, M_⋆,2), i.e., the abundance of galaxies with M_⋆,1 < M_⋆ < M_⋆,2, using

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Phi }}({M}_{\star ,1},{M}_{\star ,2}) & = & {\displaystyle \int }_{0}^{\infty }[N({M}_{{\rm{h}}}| \gt {M}_{\star ,1})-N({M}_{{\rm{h}}}| \gt {M}_{\star ,2})]\\ & & \times \displaystyle \frac{{dn}}{{{dM}}_{{\rm{h}}}}\,{{dM}}_{{\rm{h}}},\end{array}\end{eqnarray} \tag{ 12 }$$

where dn/dM_h is the halo mass function, from which we can construct the total stellar mass function. This is then convolved with a Gaussian with $\sigma ={\sigma }_{\mathrm{log}{M}_{\star }}$ to account for scatter before comparing to our observational data.

As a further constraint on our model, we use the measured two-point angular correlation functions of passive and star-forming galaxies selected by their stellar mass. Using the angular correlation function, w(θ), rather than its spatial analog, ξ(r), is a necessary consequence of our photometric redshifts, which are not accurate enough to measure the galaxy distribution in three dimensions. Therefore, a two-point spatial correlation function, ξ(r), computed from an HOD, needs to be projected along the line of sight into two dimensions for comparison with our data.

First, we compute ξ(r) from our HOD assuming Navarro–Frenk–White (Navarro et al. 1997) halo density profile, with the concentration relation of Bullock et al. (2001).⁹ The halo mass function used is the parameterization of Tinker et al. (2008), with the high-redshift correction of Behroozi et al. (2013), as well as the large-scale dark matter halo bias parameterization of Tinker et al. (2010a). These Tinker et al. relations adopt the definition of a halo as a spherical overdensity of 200 relative to the mean cosmic density at the epoch of interest (e.g., Lacey & Cole 1994; Tinker et al. 2008).

The matter power spectrum is computed according to Eisenstein & Hu (1999) with the nonlinear correction of Smith et al. (2003) applied. Additionally, we implement the two-halo exclusion model of Tinker et al. (2005), which improves on that presented in Zheng (2004).

For projecting our two-point spatial correlation functions, we use the Limber (1953) equation,

$$\begin{eqnarray}&&w(\theta )=\displaystyle \frac{\int {\left({n}_{\mathrm{gal}}(z)\tfrac{{dV}}{{dz}}W(z)\right)}^{2}\tfrac{{dz}}{d{\chi }_{{\rm{c}}}}\,{dz}\int \xi (r,z){du}}{{\left(\int {n}_{\mathrm{gal}}(z)\tfrac{{dV}}{{dz}}W(z){dz}\right)}^{2}},\end{eqnarray} \tag{ 13 }$$

where n_gal is the number density of galaxies predicted by the HOD, dV/dz is the comoving volume element, dz/dχ_c = H₀E(z)/c where E(z) = [Ω_m(1 + z)³ + Ω_Λ]^1/2, and χ_c corresponds to the comoving radial distance to redshift z. The comoving line-of-sight separation, u, is defined by $r\,={[{u}^{2}+{\chi }_{{\rm{c}}}^{2}{\varpi }^{2}]}^{1/2}$ where ${\varpi }^{2}/2=[1-\cos (\theta )]$.

The W(z) term in Equation (13) was introduced by Cowley et al. (2018) and relates to the redshift window that is being probed by the survey. If the redshifts were known precisely, then (ignoring further complications such as redshift space distortions) this would be a top-hat function equal to unity between the limits of the redshift range probed and zero elsewhere. However, with photometric redshifts, this is not the case, and the top-hat window should be convolved with an error kernel that is generally unknown a priori.

In order to mitigate this, here we assume a Gaussian error kernel and approximate a (1 + z) evolution in the error kernel dispersion, Δ_z, such that

$$\begin{eqnarray}&&W(z)=\displaystyle \frac{1}{2}\left[\mathrm{erf}\left(\displaystyle \frac{z-{z}_{\mathrm{lo}}}{{{\rm{\Delta }}}_{z}(1+{z}_{\mathrm{lo}})}\right)-\mathrm{erf}\left(\displaystyle \frac{z-{z}_{\mathrm{hi}}}{{{\rm{\Delta }}}_{z}(1+{z}_{\mathrm{hi}})}\right)\right].\end{eqnarray} \tag{ 14 }$$

Here, z_lo and z_hi represent the lower and upper redshift limits of the photometric redshift bin, respectively. While the integral in Equation (13) is, in principle, over all redshift, it only has significant contributions from redshifts where W(z) is appreciably nonzero.

Once projected according to Equation (13), we account for the integral constraint (e.g., Groth & Peebles 1977). This is the correction required as the measured angular correlation function will integrate to zero over the whole field, by construction of the estimator for w(θ) used. This results in the measured angular correlation function being underestimated by an average amount

$$\begin{eqnarray}&&{\sigma }_{\mathrm{IC}}^{2}=\displaystyle \frac{1}{{{\rm{\Omega }}}^{2}}\int \int {w}_{\mathrm{true}}(\theta )d{{\rm{\Omega }}}_{1}d{{\rm{\Omega }}}_{2},\end{eqnarray} \tag{ 15 }$$

where w_true is the true angular correlation function and the angular integrations are performed over a field of area Ω.

Here, we evaluate Equation (15) according to the numerical method proposed by Roche & Eales (1999),

$$\begin{eqnarray}&&{\sigma }_{\mathrm{IC}}^{2}=\displaystyle \frac{\sum {w}_{\mathrm{true}}(\theta )\mathrm{RR}(\theta )}{\sum \mathrm{RR}(\theta )}.\end{eqnarray} \tag{ 16 }$$

We take w_true to be the angular correlation function predicted by our HOD model and subtract ${\sigma }_{\mathrm{IC}}^{2}$ from it before comparing it to our observed correlation functions, rather than subtract a value from our observed correlation functions. We do this for each evaluation of w(θ) in our fitting procedure, which is described below. These corrections are typically ∼10⁻².

Our model has a total of 15 free parameters. Three of these (f_max,p, α_p, M_h,p) relate to f_p(M_h) (Equation (6)). A further six for each population (passive and star-forming) come from the four f_SHMR parameters (Equation (4)) (₀, M_c, β, γ) and two for satellite galaxies (B_sat, β_sat). We fix α_sat = 1 (e.g., Kravtsov et al. 2004), ${\sigma }_{\mathrm{log}{M}_{\star }}=0.2$ as discussed above and Δ_z = 0.035 (based on a comparison to our spectroscopic redshift sample) for both populations. In the future, we expect to investigate more complex models for Δ_z (e.g., a dependence on stellar mass); however, we find that our data are unable to robustly constrain these, so we have opted for the simplest model and have fixed Δ_z to be constant here.

In our fitting, we calibrate our model parameters such that we minimize the χ², which we define as

$$\begin{eqnarray}\begin{array}{rcl}{\chi }^{2} & = & \displaystyle \sum _{{\rm{p}},\mathrm{sf}}\left[\displaystyle \sum _{l}\displaystyle \frac{\mathrm{log}{{\rm{\Phi }}}_{l}^{{\rm{o}}}-\mathrm{log}{{\rm{\Phi }}}_{l}^{{\rm{m}}}}{\mathrm{log}{\sigma }_{\ {\rm{o}},l}^{2}}\right.\\ & & +\ \left.\displaystyle \sum _{k}\displaystyle \sum _{i,j}[{w}_{k}^{{\rm{o}}}{({\theta }_{i})-{w}_{k}^{{\rm{m}}}({\theta }_{i})]({C}_{k}^{-1})}_{i,j}[{w}_{k}^{{\rm{o}}}({\theta }_{j})-{w}_{k}^{{\rm{m}}}({\theta }_{j})]\right]\\ & & +\ \displaystyle \sum _{l}\displaystyle \frac{\mathrm{log}({{\rm{\Phi }}}_{l}^{{\rm{o}},\mathrm{sf}}/{{\rm{\Phi }}}_{l}^{{\rm{o}},{\rm{p}}})-\mathrm{log}({{\rm{\Phi }}}_{l}^{{\rm{m}},\mathrm{sf}}/{{\rm{\Phi }}}_{l}^{{\rm{m}},{\rm{p}}})}{\mathrm{log}{\sigma }_{\ {\rm{o}},\mathrm{sf},l}^{2}+\mathrm{log}{\sigma }_{{\rm{o}},{\rm{p}},l}^{2}}.\end{array}\end{eqnarray} \tag{ 17 }$$

Here, the indices "p" and "sf" relate to our passive and star-forming populations, l to bins of stellar mass, i and j to bins of angular separation and k to stellar mass selections. The superscripts "o" and "m" refer to observational data and model prediction, respectively. The first line relates to our stellar mass functions, the second to our two-point angular correlation functions, and the third to the ratio of our stellar mass functions.

We  obtain good results only using the two-point correlation function for $\mathrm{log}({M}_{\star }/{M}_{\odot })\gt 9.8$ and  >10.6 populations as constraints in our calibration. We have verified that including additional mass bins results in negligible changes to the resulting best-fit model, so we opt for a simpler model in that it has fewer observational constraints. This choice also means that the correlation (two-point) and stellar mass function (one-point) data contribute roughly equal amounts to the overall χ². We thus have 59 constraining data points, which results in 59 − 15 = 44 degrees of freedom (N_dof).

We use the affine-invariant python implementation for Markov chain Monte Carlo (MCMC), emcee (Foreman-Mackey et al. 2013) in order to fit our model. We employ 60 walkers for 10⁴ iterations after a burn-in phase of 2 × 10³ and test for convergence using the $\hat{R}$ statistic (Gelman & Rubin 1992). Our one-dimensional (1D) and two-dimensional (2D) marginalized posterior distributions are shown in the Appendix, along with our convergence tests. We assume flat priors for each free parameter. The best-fit values are taken to be the median of the posterior distribution, and the uncertainties the 16–84th percentiles. These are summarized in Table 1 in the Appendix. Our best-fit model achieves a reduced χ² of χ²/N_dof ∼ 1.9, indicating a reasonable fit to the data.

In Figure 1, we show the observed galaxy stellar mass functions for our passive and star-forming populations, as well as the ratio between the two, as with the results from our best-fit model shown for comparison. The shape of the high-mass end of the stellar mass functions are substantively similar and the main difference between them is that the slope of the function for passive galaxies at masses smaller than the "knee" is much flatter. Both stellar mass functions can be accurately reproduced by the model.

Zoom In Zoom Out Reset image size

In Figure 2, we present the two-point angular correlation functions for passive and star-forming galaxies, with results from the best-fit model also shown. The large-scale clustering (i.e., two-halo term) of the two populations is similar. However, there is more power on small scales for the passive population, which suggests a larger fraction of these galaxies are satellites (we will explore this further in Section 4.3). While passive galaxies are recognized as being more clustered in general at lower redshift (e.g., Zehavi et al. 2011), this appears to be less pronounced on larger scales and with increasing redshift (e.g., Hartley et al. 2010). We also note that Tinker et al. (2013) find similar differences between their passive and star-forming correlation functions for 0.2 < z < 1.0. Again, we can see that the model can reasonably reproduce the observed correlation functions, though there appear to be some mild systematic discrepancies at (i) ∼10⁻² deg for passive galaxies, which could be related to complex nonlinear biases (this scale approximately coincides with the transition from the one- to two-halo term) not properly accounted for in our model; and (ii) at large scales (≳2 × 10⁻² deg) for star-forming galaxies with $\mathrm{log}({M}_{\star }/{M}_{\odot })\lesssim 10.2$. Here, the discrepancy may be related to our decision to fix Δ_z, which can control the amplitude of the correlation function independent of the halos the galaxies occupy (e.g., Cowley et al. 2018). We have investigated more complex models in which Δ_z was a free parameter for both passive and star-forming samples and had a stellar mass dependence. However, as mentioned earlier, we find that our current data are insufficient to robustly constrain these models, although this may help alleviate this mild discrepancy.

Zoom In Zoom Out Reset image size

In general, these discrepancies are fairly minor and the results shown here, combined with the reduced chi-squared value of χ²/N_dof ∼ 1.9 mentioned earlier, indicate that our model can reproduce the calibration data to a good degree of accuracy.

Figure 3 shows the best-fit SHMR, ε(M_h) (Equation (5)), for passive and star-forming central galaxies and the fraction of central galaxies that are passive as a function of host halo mass (Equation (6)). The high-mass slopes of the ratios are similar for both populations, as would be expected from the similarity in the high-mass regime of the stellar mass functions. Interestingly, the star-forming population relation peaks at a larger halo mass (by ∼0.5 dex) than the passive one (though at a similar stellar mass, ∼10^10.7 M_⊙, roughly coinciding with the knee of the observed stellar mass functions) and the passive relation has a higher normalization for M_h ≲ 10¹³ M_⊙. We now investigate this difference in normalization further using a simple model comprising a realistic ensemble of halo histories.

Zoom In Zoom Out Reset image size

McBride et al. (2009) found that the mass assembly histories of dark matter halos in the Millennium N-body simulations (Springel et al. 2005) could be well described with the two parameter functional form

$$\begin{eqnarray}&&M{(z)={M}_{0}(1+z)}^{{\beta }_{{\rm{h}}}}\,{e}^{-{\gamma }_{{\rm{h}}}z},\end{eqnarray} \tag{ 18 }$$

where M₀ is the halo mass at z = 0. We generate a realistic ensemble of β_h and γ_h values according to Appendix A of McBride et al. and weight the halos such that our halo mass function at z = 2.3 is reproduced. From this, we generate the star formation histories, ψ(M_h, z), in each halo according to

$$\begin{eqnarray}&&\psi ({M}_{{\rm{h}}},z)={\varepsilon }^{{\prime} }({M}_{{\rm{h}}},z)\,{f}_{{\rm{b}}}\,{\dot{M}}_{{\rm{h}}}\end{eqnarray} \tag{ 19 }$$

(e.g., Tacchella et al. 2018) where f_b = Ω_b/Ω_m is the baryon fraction and ε'(M_h, z) has the same form as Equation (5), though has a different meaning than in Equation (4) so is denoted with the "prime" symbol here for clarity, as are its parameters. Additionally, we compute the stellar mass of each galaxy, assuming the instantaneous recycling approximation, through

$$\begin{eqnarray}&&{M}_{\star }=(1-R)\int \,\psi ({M}_{{\rm{h}}},z){dt},\end{eqnarray} \tag{ 20 }$$

where R is the returned fraction, the fraction of initial stellar mass that is returned to the interstellar medium by mass loss from dying stars. We use a value of R = 0.41, which corresponds to our choice of a Chabrier (2003) IMF. For computing stellar mass functions, we add a Gaussian scatter with ${\sigma }_{\mathrm{log}{M}_{\star ,\mathrm{meas}}}=0.13$ dex onto our intrinsic stellar masses, as we did in Section 3. In bins of stellar mass, we then label galaxies with the lowest intrinsic specific star formation rates as passive such that the difference between the passive and star-forming central galaxy stellar mass functions is reproduced. We assume no redshift evolution in ε'(M_h, z) and find that the parameter values (${\epsilon }_{0}^{{\prime} }$, ${M}_{{\rm{c}}}^{{\prime} }$, β', γ') = (0.15, 2.1 × 10¹² h⁻¹ M_⊙, 0.8, 0.8) give reasonable results, i.e., can reproduce our observed stellar mass function.

The SHMRs for passive and star-forming galaxies from this model are shown in the top panel of Figure 4. They are essentially identical, implying that for the difference in normalization we found in Figure 3 to arise ${\varepsilon }^{{\prime} }({M}_{{\rm{h}}},z)$ must evolve with redshift. Indeed, we illustrate this further by repeating the above exercise but now allowing ${\epsilon }_{0}^{{\prime} }$ to vary as a function of redshift. For simplicity we choose ${\epsilon }_{0}^{{\prime} }(z)\,={\epsilon }_{{\rm{N}}}^{{\prime} }+{\epsilon }_{z}^{{\prime} }\times (1+z)$ and values of (${\epsilon }_{{\rm{N}}}^{{\prime} }$, ${\epsilon }_{z}^{{\prime} })=(0.04,0.03)$ such that the conversion of baryons into stars becomes more efficient at higher redshifts (and our observed stellar mass functions are still reasonably reproduced), all other parameters are kept the same. The results from this model are shown in the bottom panel of Figure 4. As expected, we find that the normalization of the SHMR for passive galaxies is elevated relative to that of the star-forming galaxies. We note that we have not precisely reproduced the differences found in Figure 3, as the evolution of ε'(M_h, z) is likely more complex than investigated here, but we have illustrated how such differences can occur.

Zoom In Zoom Out Reset image size

For each halo history, we now calculate a halo formation redshift, z_f, which we define as the redshift at which that halo had assembled half of its z = 2.3 mass. We show these for passive and star-forming galaxies from the evolving ε'(M_h, z) model in Figure 5, though we note that this figure would be substantively similar if the constant ε'(M_h, z) model were considered. Interestingly, selecting passive galaxies by choosing those with the lowest specific star formation rates results in us selecting halos with the highest formation redshifts, i.e., those that are assembling mass at the lowest rates¹⁰ . This helps further explain the result shown in Figure 4. Passive galaxies tend to reside in "older" halos, thus gaining a larger proportion of their final stellar mass at higher redshifts when the conversion from baryons into stars was more efficient. This results in a higher SHMR, such as that found in Figure 3.

Zoom In Zoom Out Reset image size

Additionally, on further inspection of Figure 5, we see that the halo formation redshift for star-forming galaxies is not strongly dependent on stellar mass but for passive galaxies, there is a much stronger trend that flattens for M_⋆ ≳ 10^10.75 M_⊙. This could explain the difference in the evolution of the number density of intermediate- and high-mass passive galaxies found by Deshmukh et al. (2018), as these stellar mass ranges appear to select halos with very different rates of growth.

In this section, we address the abundance of satellites in our galaxy populations and discuss potential physical quenching mechanisms for them.

As a result of our halo formalism, we can compute the contribution to the stellar mass function from a given halo mass by appropriately changing the integration limits in Equation (12). We show the contribution to the stellar mass function from a selection of halo masses and divided into passive/star-forming and central/satellite in Figure 6. Here we can see that at a fixed stellar mass the satellite mass function is dominated by greater mass halos than the central one, and this is especially pronounced for passive galaxies. The stellar mass of central galaxies forms a relatively tight relationship with the mass of their host halos (e.g., Equation (5)) but this becomes offset once they are accreted into a more massive host halo and become a satellite. Their stellar mass will not necessarily change much following the halo merger but their host halo mass will suddenly have increased, so this result is unsurprising. Additionally, we can see that the contributions from each halo mass range form a fairly well-defined peaked function for central galaxies, but this is not the case for satellites, indicating that a relationship between the stellar mass of a satellite galaxy and the mass of its host halo must have a much greater scatter than for centrals. This can be understood as the merger event and subsequent environmental effects scattering satellite galaxies from the initial (relatively tight) SHMR they had as central galaxies.

Zoom In Zoom Out Reset image size

We show the satellite fractions for our passive population and the complete mass selected population (passive + star-forming) in Figure 7. We see that our results for the complete population compare favorably with our earlier analysis in Cowley et al. (2018), which used a less-sophisticated halo occupation model than we have adopted here. We can also see that the satellite fractions for passive galaxies are generally greater (as has been found at lower redshifts, e.g., Knobel et al. 2013), though the difference is fairly small and there are significant uncertainties on the satellite fraction for passive galaxies. This is unsurprising, as this quantity is constrained mainly by the small-scale clustering data for passive galaxies that have significant uncertainties. However, a greater satellite fraction indicates that a passive galaxy is more likely to be a satellite than one selected from the combined population, and implies that  the environmental effects satellite galaxies are subject to can play a role in quenching ongoing star formation. We explore this idea in greater detail below.

Zoom In Zoom Out Reset image size

4.3.1. The Astrophysics of Satellite Quenching

In this Section, we infer satellite quenching timescales and compare these with other estimates from the literature at lower redshift, to briefly discuss the possible physical mechanisms for satellite quenching, i.e., dynamical/orbit- or outflow-based, which can generally be recognized as "environmental" and "intrinsic/mass-related", respectively. In doing so, we are broadly revisiting the discussion of McGee et al. (2014), though arrive at a different conclusion.

We derive a quenching timescale for our satellite population based on the method of Tinker & Wetzel (2010) and Tinker et al. (2013). At each stellar mass, we compute the fraction of satellite galaxies that are passive, and then compare this to the age distribution of satellite halos (i.e., the time since they were last a distinct halo) at this redshift taken from the recent P-Millennium simulation (Baugh et al. 2018). Using the ansatz that the quenched (i.e., passive) galaxies are the ones that have been satellites for the longest period of time, we can then derive a quenching timescale. For example, if 50% of satellite galaxies are passive and 50% of satellite halos have been a satellite for 1 Gyr or more, then our quenching timescale would be 1 Gyr. However, some satellites may have been passive prior to accretion onto their current host. To account for this, we follow Tinker et al. (2013) and assume two extremes: either that no galaxies were passive prior to being accreted (v1) or that the fraction of central galaxies that are passive at z = 2.3 was already passive¹¹ (v2). For example, if 50% of satellite galaxies are passive but 10% of central galaxies are passive, then the quenching timescale we are concerned with in v2 will be the 40th percentile of the satellite age distribution, whereas in v1 it would be the 50th percentile.

In the case that quenching depends only on the orbit of the satellite within its host halo then the quenching timescale will evolve as the ratio of the inverse densities, ∝(1 + z)^−3/2 [for a halo in virial equilibrium the dynamical timescales as ρ^−1/2 and the density, ρ, scales as (1 + z)³]. However, if supernovae-driven outflows are primarily responsible then the quenching time will scale more strongly with the star formation rate. McGee et al. used the fitted prescriptions for the star-forming "main sequence" of Peng et al. (2010) and Whitaker et al. (2012) to investigate the potential evolution of the quenching timescale under this scenario, assuming that the outflow rate is invariant with redshift for a given star formation rate and stellar mass. They also modified the prescription of Whitaker et al. to account for different mass-loading factors, η, where a larger value means more cold gas mass is ejected from a galaxy per unit star formation, thus inhibiting subsequent star formation. We show the scaling of these outflow models [normalized to the quenching timescale of Wetzel et al. (2013) for M_⋆ = 10^10.5 M_⊙ galaxies at z ∼ 0] in Figure 8. We also include in this figure a broad range of other estimates of the quenching timescale of 10^10.5 M_⊙ satellite galaxies, including our own. At M_⋆ = 10^10.5 M_⊙ we find ${t}_{\mathrm{quench}}={0.77}_{-0.22}^{+0.33}$ Gyr for v1 and ${1.03}_{-0.27}^{+0.60}$ Gyr for v2.

Zoom In Zoom Out Reset image size

The various observational data seem to cover a plethora of possible scenarios, most likely related to the use of different techniques and methods to derive a quenching timescale. This highlights the need for a homogeneous analysis to be performed over a broad range of redshifts and means that we are not able to draw strong conclusions here. McGee et al. considered only the data of Wetzel et al. (2013), McGee et al. (2011), Mok et al. (2014), and Muzzin et al. (2014) and so arrived at the conclusion that outflow-based models were in better agreement with the data. In contrast, our estimates, taken in conjunction with those of Wetzel et al. (2013), seem to strongly favor a dynamical mechanism for satellite quenching, e.g., ram pressure stripping. This is the same conclusion that Tinker & Wetzel (2010) arrived at and it appears from inspection of Figure 8 that more high-redshift data, such as that presented in this work, could help further discriminate between these two potential processes.