The Stellar-to-Halo Mass Ratios of Passive and Star-Forming Galaxies at z~2-3 from the SMUVS survey

TL;DR

This study investigates the stellar-to-halo mass ratios of passive and star-forming galaxies at redshifts 2-3 using galaxy clustering data, revealing differences in their halo properties and quenching mechanisms.

Contribution

It provides the first detailed analysis of SHMRs for passive and star-forming galaxies at high redshift, incorporating halo modeling and galaxy population division.

Findings

Passive galaxies have higher SHMR normalization than star-forming ones.

Passive galaxies likely reside in halos with higher formation redshifts.

Satellite galaxies occupy more massive halos than centrals at fixed stellar mass.

Abstract

In this work, we use measurements of galaxy stellar mass and two-point angular correlation functions to constrain the stellar-to-halo mass ratios (SHMRs) of passive and \np\ galaxies at $z\sim2-3$, as identified in the \emph{Spitzer} Matching Survey of the UltraVISTA ultra-deep Stripes (SMUVS). We adopt a sophisticated halo modeling approach to statistically divide our two populations into central and satellite galaxies. For central galaxies, we find that the normalization of the SHMR is greater for our passive population. Through the modeling of $\Lambda$ cold dark matter halo mass accretion histories, we show that this can only arise if the conversion of baryons into stars was more efficient at higher redshifts and additionally that passive galaxies can be plausibly explained as residing in halos with the highest formation redshifts (i.e., those with the lowest accretion rates) at a given halo mass. At a fixed stellar mass, satellite galaxies occupy host halos with a greater mass than central galaxies, and we find further that the fraction of passive galaxies that are satellites is higher than for the combined population. This, and our derived satellite quenching timescales, combined with earlier estimates from the literature, support dynamical/environmental mechanisms as the dominant process for satellite quenching at $z\lesssim3$.

Figures

Tables

Table 1. Table 1: Parameter values from MCMC fitting

Parameter	Value	Prior Range	Eqn.
Passive fraction:
$\log \left(M_{h,q}\right)$	${13.03}_{-0.55}^{+0.51}$	$(12.00,14.00)$	6
$f_{\max ,p}$	${0.52}_{-0.24}^{+0.31}$	$(0.00,1.00)$
${\sigma }_q$	${1.02}_{-0.43}^{+0.41}$	$(0.10,2.00)$
Passive Galaxies:
$\log \left({ϵ}_0\right)$	$-{1.59}_{-0.21}^{+0.30}$	$(-2.50,-0.82)$	5
$\log \left(M_c\right)$	${11.83}_{-0.37}^{+0.24}$	$(11.00,12.50)$
$\beta$	${2.72}_{-1.64}^{+1.69}$	$(0.00,6.00)$
$\gamma$	${0.71}_{-0.15}^{+0.15}$	$(0.40,1.10)$
$\log \left(B_{\mathrm{sat}}\right)$	${1.94}_{-0.19}^{+0.28}$	$(1.50,3.00)$	9
${\beta }_{\mathrm{sat}}$	${0.95}_{-0.44}^{+0.54}$	$(0.00,2.00)$
Star-Forming Galaxies:
$\log \left({ϵ}_0\right)$	$-{1.94}_{-0.03}^{+0.03}$	$(-2.50,-0.82)$	5
$\log \left(M_c\right)$	${12.60}_{-0.07}^{+0.07}$	$(11.50,12.80)$
$\beta$	${0.75}_{-0.06}^{+0.06}$	$(0.00,2.50)$
$\gamma$	${0.90}_{-0.15}^{+0.08}$	$(0.40,1.10)$
$\log \left(B_{\mathrm{sat}}\right)$	${1.42}_{-0.04}^{+0.04}$	$(0.50,2.50)$	9
${\beta }_{\mathrm{sat}}$	${1.62}_{-0.21}^{+0.24}$	$(0.50,2.00)$
Fixed Parameters:
${\alpha }_{\mathrm{sat}}$	$1.0$	–	8
${\sigma }_{\log \left(M_{\star }\right)}$	$0.2$ dex	–	–
${\Delta }_z$	$0.035$	–	14

Equations

$${\delta }^2{P}_{12}={\bar{\eta }}^2[1+w(\theta )]\delta {\Omega }_1\delta {\Omega }_2,$$

$$w(\theta )=\frac{DD(\theta )-2DR(\theta )+RR(\theta )}{RR(\theta )},$$

$$C_{i,j}=\frac{N-1}{N}\sum _{l=1}^N(w_i^l-\overline{w_i})\times (w_j^l-\overline{w_j}),$$

$$\log f_{\mathrm{SHMR}}=\log [M_{\star }(M_{\mathrm{h}})]=\log [\epsilon (M_{\mathrm{h}})\times M_{\mathrm{h}}]$$

$$\epsilon (M_{\mathrm{h}})=2{ϵ}_0{\left[{\left(\frac{M_{\mathrm{h}}}{M_{\mathrm{c}}}\right)}^{-\beta }+{\left(\frac{M_{\mathrm{h}}}{M_{\mathrm{c}}}\right)}^{\gamma }\right]}^{-1}$$

$$f_{\mathrm{p}}(M_{\mathrm{h}})=f_{\max ,\mathrm{p}}\left[\tanh \left({\alpha }_{\mathrm{p}}[\log (M_{\mathrm{h}})-\log (M_{\mathrm{h},\mathrm{p}})]\right)+1\right]/2,$$

$$N_{\mathrm{c}}(M_{\mathrm{h}}|>M_{\star })=\frac{1}{2}\left[1-\mathrm{erf}\left(\frac{\log M_{\star }-\log f_{\mathrm{SHMR}}}{2{\sigma }_{\log M_{\star }}}\right)\right],$$

$$N_{\mathrm{s}}(M_{\mathrm{h}}|>M_{\star })={\left(\frac{M_{\mathrm{h}}}{M_{\mathrm{sat}}}\right)}^{{\alpha }_{\mathrm{sat}}}\exp \left(-\frac{M_{\mathrm{cut}}}{M_{\mathrm{h}}}\right),$$

$$\frac{M_{\mathrm{sat}}}{10^{12}{\mathrm{M}}_{\odot }}=B_{\mathrm{sat}}{\left(\frac{f_{\mathrm{SHMR}}^{-1}}{10^{12}{\mathrm{M}}_{\odot }}\right)}^{{\beta }_{\mathrm{sat}}}.$$

$$\log M_{\mathrm{cut}}=0.76\log M_{\mathrm{sat}}+2.3,$$

$$N(M_{\mathrm{h}}|>M_{\star })=N_{\mathrm{c}}(M_{\mathrm{h}}|>M_{\star })+N_{\mathrm{s}}(M_{\mathrm{h}}|>M_{\star }).$$

$$\Phi (M_{\star ,1},M_{\star ,2})$$

$$w(\theta )=\frac{\int {\left(n_{\mathrm{gal}}(z)\frac{\mathrm{d}V}{\mathrm{d}z}W(z)\right)}^2\frac{\mathrm{d}z}{\mathrm{d}{\chi }_{\mathrm{c}}}\mathrm{d}z\int \xi (r,z)\mathrm{d}u}{{\left(\int n_{\mathrm{gal}}(z)\frac{\mathrm{d}V}{\mathrm{d}z}W(z)\mathrm{d}z\right)}^2},$$

$$W(z)=\frac{1}{2}\left[\mathrm{erf}\left(\frac{z-z_{\mathrm{lo}}}{{\Delta }_z(1+z_{\mathrm{lo}})}\right)-\mathrm{erf}\left(\frac{z-z_{\mathrm{hi}}}{{\Delta }_z(1+z_{\mathrm{hi}})}\right)\right].$$

$${\sigma }_{\mathrm{IC}}^2=\frac{1}{{\Omega }^2}\int \int w_{\mathrm{true}}(\theta )\mathrm{d}{\Omega }_1\mathrm{d}{\Omega }_2,$$

$${\sigma }_{\mathrm{IC}}^2=\frac{\sum w_{\mathrm{true}}(\theta )RR(\theta )}{\sum RR(\theta )}.$$

$${\chi }^2$$

$$M(z)=M_0(1+z{)}^{{\beta }_{\mathrm{h}}}{e}^{-{\gamma }_{\mathrm{h}}z},$$

$$\psi (M_{\mathrm{h}},z)={\epsilon }^{\prime }(M_{\mathrm{h}},z)f_{\mathrm{b}}{\dot{M}}_{\mathrm{h}}$$

$$M_{\star }=(1-R)\int \psi (M_{\mathrm{h}},z)\mathrm{d}t,$$

Full text

methods: statistical – galaxies: evolution – galaxies: formation – galaxies: high-redshift – large-scale structure of Universe

††journal: ApJ

1 Introduction

One of the most striking properties of the galaxy population at $z=0$ is its distinct optical color bi-modality (e.g., Strateva et al., 2001; Baldry et al., 2004). Broadly speaking, galaxies are either blue or red, commonly interpreted as star-forming or passively evolving (meaning the current ongoing star formation is significantly lower than the past average). respectively, with a transitional ‘green’ valley (e.g., Schawinski et al., 2014) in between. This bi-modality is observed to be in place until at least $z\sim 1$ (e.g., Bell et al., 2004) and possibly persists at higher redshifts (e.g., Cirasuolo et al., 2007). Although the competing effect of dust obscuration becoming increasingly important at higher redshifts makes the interpretation of the optical color less clear (e.g., Stern et al., 2006).

The relative importance of various physical processes through which this bi-modality emerges is not immediately clear. However, over the last few years, a broad paradigm of two main channels for quenching star formation has emerged. These can be summarized as ‘mass’ (sometimes referred to as ‘intrinsic’) and ‘environment’ quenching and have been plausibly shown to explain the observed Schechter function shape of the local galaxy stellar mass function (Peng et al., 2010). We discuss these two main channels in greater detail in the following two paragraphs.

For mass/intrinsic quenching, some early galaxy formation models produced a distinct red sequence through the incorporation of radio-mode active galactic nucleus (AGN) feedback (e.g. Bower et al., 2006; Croton et al., 2006) where star formation is suppressed in relatively high-mass halos ($\gtrsim 10^{12}$ M*⊙*) in which gas in the hot halo is undergoing quasi-static cooling. However, the bi-modality in the model of Menci et al. (2005) was produced by the mass assembly history of dark matter halos, as their model did not include a prescription for AGN feedback. It should also be noted that whilst physical models can produce a distinct red sequence, they often struggle to precisely reproduce the observed distribution of $g-r$ color locally (e.g., Trayford et al., 2015; Lacey et al., 2016).

It is also accepted that environmental effects acting dynamically on satellite galaxies orbiting within the virial radius of a host halo also play a role in producing a population of passive galaxies. These environmental mechanisms include ram pressure stripping of a galaxy’s hot gas reservoir as it enters that of a larger host halo, strangulation as gas no longer cools onto the satellite galaxy and harassment, tidal interactions between other nearby galaxies (e.g. van den Bosch et al., 2008). Though, as we will discuss later, outflows from supernovae in the satellite may also be important (McGee et al., 2014).

In this work, we use a halo occupation formalism (e.g. Zheng et al., 2005) to investigate the differences between passive and star-forming galaxies at $z=2-3$ as identified in the Spitzer Matching survey of the UltraVISTA ultra-deep Stripes (SMUVS; Ashby et al., 2018; Deshmukh et al., 2018, see also Caputi et al. 2017). This halo formalism provides a statistical framework that describes the expected mean number of galaxies occupying host halos of a given mass. More specifically, we use a halo model similar to that introduced by Leauthaud et al. (2011, see also Tinker et al. 2013), which begins with a parameterization of the stellar-to-halo mass ratio (SHMR) for central galaxies, a function that specifies the mean stellar mass of a central galaxy as a function of halo mass. The free parameters in this model are then calibrated using the observed galaxy stellar mass and two-point angular correlation functions. Although this is a statistical model that lacks a direct physical interpretation, it provides us with a means to investigate the link between stellar and halo mass for passive and star-forming galaxies (divided into central and satellite galaxies). As such, the model’s predictions can then be understood in the context of the potential physical processes involved.

We restrict our analysis to $z=2-3$ as the SMUVS photometry is optimized for $z\gtrsim 2$ galaxies (Deshmukh et al., 2018). At higher redshifts ($z\gtrsim 3$), the number of passive galaxies is insufficient to robustly measure the angular two-point correlation function for this population, which is a key constraint on our model. We thus extend similar analyses that have been performed at lower redshifts using a combination of galaxy-galaxy lensing, galaxy two-point correlation functions, stellar mass functions and halo model analyses (e.g. Mandelbaum et al., 2006; van Uitert et al., 2011; Tinker et al., 2013; Rodríguez-Puebla et al., 2015). These earlier works generally find that local ($z\sim 0$) differences in the SHMR between passive and star-forming galaxies are generally small, with blue/star-forming galaxies having an SHMR a factor of $\sim 2$ higher than red/passive galaxies (e.g. Rodríguez-Puebla et al., 2015). However, this difference diminishes toward higher redshift and actually inverts such that at a fixed stellar mass, star-forming galaxies are in more massive halos by $z\sim 1$ (e.g. Tinker et al., 2013), particularly for galaxies with high stellar masses ($\gtrsim 10^{10.5}$ M*⊙*).

Throughout this work, we assume a flat $\Lambda$ cold dark matter (CDM) cosmological model with ($\Omega_{\rm m}$, $\Omega_{\rm b}$, $\Omega_{\Lambda}$, $h$, $\sigma_{8}$, $n_{\rm s}$) = ($0.3$, $0.04$, $0.7$, $0.8$, $0.95$). This paper is structured as follows: in Section 2, we summarize the main details of the SMUVS survey, describe the data used and the determination of our observational constraints, the galaxy stellar mass, and angular two-point correlation functions. In Section 3, we describe our halo occupation model, and we present our main results in Section 4. We summarize and conclude in Section 5. Some further information regarding our model fitting is given Appendix A.

2 The Data

2.1 The SMUVS survey

The SMUVS program (PI: K. Caputi; Ashby et al. 2018) has collected ultra-deep Spitzer $3.6$ and $4.5$ $\mu$m data within the COSMOS111http://cosmos.astro.caltech.edu 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,222http://www.eso.org/sci/observing/phase3/data_releases/uvista_dr3.pdf which reaches an average depth of $K_{\rm s}=24.9\pm 0.1$ and $H=25.1\pm 0.1$ ($2$ arcsec diameter, $5\sigma$).

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_{\rm 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 (PSF) fitting photometric measurements on the SMUVS $3.6$ and $4.5$ $\mu$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_{\rm 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 $\sim 2.9\times 10^{5}$ UltraVISTA sources with a detection in at least one IRAC band over an area of $\sim 0.66$ square degrees.

The SED fitting is performed with all $28$ bands ($26$ $U$ through $K_{\rm s}$ as well as Spitzer $3.6$ and $4.5$ $\mu$m) using the template-fitting code LePhare333http://www.cfht.hawaii.edu/~arnouts/LEPHARE/lephare.html (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, $\sigma_{z}$, of $|z_{\rm phot}-z_{\rm spec}|/(1+z_{\rm spec})$, based on $\sim 1.4\times 10^{4}$ 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\sim 2-3$ population. This statistic is computed excluding outliers, defined as objects for which $|z_{\rm phot}-z_{\rm spec}|/(1+z_{\rm spec})>0.15$, which comprise only $\sim 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\sim 2-3$ galaxies, we find $\sigma_{z}=0.035$.

Throughout this paper, we use the best-fit redshifts and stellar masses computed by LePhare. We also limit ourselves to $\log(M_{\star}/\mathrm{M}_{\odot})>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 $\log(M_{\star}/\mathrm{M}_{\odot})>9.6$ at $z\sim 2-3$ (Deshmukh et al., 2018). This stellar mass cut leaves us with $\sim 10^{3}$ passive and $\sim 6.8\times 10^{3}$ star-forming galaxies with $2<z<3$.

2.2 Definition of passive galaxies

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)\leq 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\geq 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.

2.3 The galaxy stellar mass function

We compute the galaxy stellar mass functions for our passive and star-forming populations using the $1/V_{\mathrm{max}}$ technique (Schmidt, 1968) as was done by Deshmukh et al. (2018). A $V_{\mathrm{max}}$ correction is computed for each galaxy based on the maximum volume at which it would have a magnitude $[4.5$ $\mu\mathrm{m}]=26$ (or $[3.6$ $\mu\mathrm{m}]=26$ in the case of a non-detection at 4.5 $\mu$m). For fainter sources, no $V_{\mathrm{max}}$ correction is applied. A further completeness correction is applied to each galaxy based on its $4.5$ $\mu$m magnitude (or $3.6$ $\mu$m in the case of a non-detection at $4.5$ $\mu$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, $\sigma_{\mathrm{Poi}}$, an error associated with photometric uncertainties, $\sigma_{\mathrm{mc}}$, and a cosmic variance error, $\sigma_{\mathrm{cv}}$, in quadrature. The Poisson term, $\sigma_{\rm Poi}$, is derived using the tabulated values of Gehrels (1986). The photometric term, $\sigma_{\rm 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; $\sigma_{\rm mc}$ is then taken to be the $16-84^{\rm th}$ 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 Fig. 19) provides a thorough comparison of $2<z<6$ SMUVS-derived galaxy stellar mass functions to earlier estimates in the literature.

2.4 The two-point angular correlation function

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

[TABLE]

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

[TABLE]

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 ($\gtrsim 95$ %) we do not need to consider any variations in sensitivity/depth across our survey. We use a single random catalog with $\sim 5\times 10^{5}$ objects ($\sim 50\times$ 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

[TABLE]

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 $l^{\rm th}$ region removed (i.e. computed from 59 of the 60 regions). We compute $w(\theta)$ using eight evenly spaced logarithmic bins of angular separation in the range $-3<\log_{10}(\theta/\mathrm{deg})<-1$.

3 The Theoretical Model

The model used here is based on that presented by Tinker et al. (2013), which in turn extends on that first presented by Leauthaud et al. (2011). We begin by defining the mean stellar-to-halo mass relation, $f_{\rm SHMR}$, for central galaxies. For this we adopt

[TABLE]

where the SHMR, $\varepsilon(M_{\rm h})$, is described by

[TABLE]

That is, a double power law with a low-mass slope, $\beta$, and a high-mass slope, $-\gamma$, normalized to $\epsilon_{0}$ at halo mass $M_{\rm c}$ (e.g., Yang et al., 2003; Moster et al., 2010). More complex forms for this relation have been suggested in the literature (e.g., Behroozi et al., 2013); however, we find that a double power law is sufficient for our purposes. We do not place any constraint on $f_{\rm SHMR}$ other than the requirement that it is monotonic with halo mass such that it can be inverted, which in practice mainly limits the possible values of $\gamma$. We allow for an intrinsic Gaussian scatter from this relation, $\sigma_{\log M_{\star\rm,int}}$, of $0.16$ dex, which we assume is constant with respect to halo mass (e.g., More et al., 2011). We also assume a measurement-based Gaussian scatter of $\sigma_{\log M_{\star\rm,meas}}=0.04\times(1+z)$ as suggested by Ilbert et al. (2013), which is $0.13$ dex based on a median redshift of $z=2.3$ for our $\log(M_{\star}/M_{\odot})>9.8$ galaxies. We thus use a total scatter of $\sigma_{\log M_{\star}}=\left[\sigma_{\log M_{\star\rm,int}}^{2}+\sigma_{\log M_{\star\rm,meas}}^{2}\right]^{1/2}\approx 0.2$ dex. For simplicity, we apply this fixed value to both passive and star-forming galaxies.

At each halo mass, we assume that a fraction, $f_{\rm p}$, of our central galaxies are passive. For this we assume a smooth hyperbolic tan function

[TABLE]

that transitions between having zero and a fraction of $f_{\rm max,p}$ passive central galaxies, with a fraction of $f_{\rm max,p}/2$ galaxies quenched at $M_{\rm h,p}$ and the width of the transition being controlled by $\alpha_{\rm p}$. The choice of this functional form is somewhat ad hoc and exponential functional forms have also been used elsewhere in the literature (e.g. Tinker & Wetzel, 2010), though we find that implementing these here have a minimal impact on our results.

The halo occupation distribution (HOD) function of central galaxies is described by

[TABLE]

(Leauthaud et al., 2011). For passive central galaxies, this is multiplied by $f_{\rm p}(M_{\rm h})$ and by $1-f_{\rm p}(M_{\rm h})$ for star-forming centrals.

The occupation of satellite galaxies is modeled as a power law with an exponential cut-off,

[TABLE]

(Leauthaud et al., 2011) where $\alpha_{\rm sat}$ is the power-law slope the characteristic satellite halo mass, $M_{\rm sat}$, is described by

[TABLE]

Here, $f_{\rm SHMR}^{-1}$ is the inverse of Equation 4 and $B_{\rm sat}$ and $\beta_{\rm sat}$ describe the normalization and slope of the power law, respectively. The parameter $M_{\rm cut}$ in Equation 8 is often modeled in a similar fashion. However, for simplicity, here we fix this parameter according to

[TABLE]

following Conroy et al. (2006). Our combined HOD is the sum of the central and satellite distributions, i.e.,

[TABLE]

3.1 The stellar mass function

From the HOD (Equation 11), we can calculate the stellar mass function $\Phi(M_{\star,1},M_{\star,2})$, i.e., the abundance of galaxies with $M_{\star,1}<M_{\star}<M_{\star,2}$, using

[TABLE]

where $\mathrm{d}n/\mathrm{d}M_{\rm 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_{\log M_{\star}}$ to account for scatter before comparing to our observational data.

3.2 The two-point angular correlation function and the effect of photometric redshift errors

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(\theta)$, rather than its spatial analog, $\xi(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, $\xi(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 $\xi(r)$ from our HOD assuming Navarro-Frenk-White (NFW; Navarro et al., 1997) halo density profile, with the concentration relation of Bullock et al. (2001).444Martinez-Manso et al. (2015) found that the predictions of halo models, such as the one used here, are fairly insensitive to the adopted concentration relation (see their Appendix D). 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,

[TABLE]

where $n_{\rm gal}$ is the number density of galaxies predicted by the HOD, $\mathrm{d}V/\mathrm{d}z$ is the comoving volume element, $\mathrm{d}z/\mathrm{d}\chi_{\rm c}=H_{0}E(z)/c$ where $E(z)=[\Omega_{\mathrm{m}}(1+z)^{3}+\Omega_{\Lambda}]^{1/2}$, and $\chi_{\rm 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, $\Delta_{z}$, such that

[TABLE]

Here, $z_{\rm lo}$ and $z_{\rm hi}$ represent the lower and upper redshift limits of the photometric redshift bin, respectively. Whilst 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(\theta)$ used. This results in the measured angular correlation function being underestimated by an average amount

[TABLE]

where $w_{\rm true}$ is the true angular correlation function and the angular integrations are performed over a field of area $\Omega$.

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

[TABLE]

We take $w_{\rm true}$ to be the angular correlation function predicted by our HOD model and subtract $\sigma_{\rm 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(\theta)$ in our fitting procedure, which is described below. These corrections are typically $\sim 10^{-2}$.

3.3 Fitting

Our model has a total of $15$ free parameters. Three of these ($f_{\rm max,p}$, $\alpha_{\rm p}$, $M_{\rm h,p}$) relate to $f_{\rm p}(M_{\rm h})$ (Equation 6). A further six for each population (passive and star-forming) come from the four $f_{\rm SHMR}$ parameters (Equation 4) ($\epsilon_{0}$, $M_{\rm c}$, $\beta$, $\gamma$) and two for satellite galaxies ($B_{\rm sat}$, $\beta_{\rm sat}$). We fix $\alpha_{\rm sat}=1$ (e.g., Kravtsov et al., 2004), $\sigma_{\log M_{\star}}=0.2$ as discussed above and $\Delta_{z}=0.035$ (based on a comparison to our spectroscopic redshift sample) for both populations. In future, we expect to investigate more complex models for $\Delta_{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 $\Delta_{z}$ to be constant here.

In our fitting, we calibrate our model parameters such that we minimize the $\chi^{2}$, which we define as

[TABLE]

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 $\log(M_{\star}/\mathrm{M}_{\odot})>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 $\chi^{2}$. We thus have $59$ constraining data points, which results in $59-15=44$ degrees of freedom ($N_{\rm 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^{4}$ iterations after a burn-in phase of $2\times 10^{3}$ and test for convergence using the $\hat{R}$ statistic (Gelman & Rubin, 1992). Our one-dimensional and two-dimensional marginalised posterior distributions are shown in Appendix A, 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-84^{\rm th}$ percentiles. These are summarized in Table 1 in Appendix A. Our best-fit model achieves a reduced $\chi^{2}$ of $\chi^{2}/N_{\rm dof}\sim 1.9$, indicating a reasonable fit to the data.

4 Results

In this section, we present our main results. The performance of the best-fit model with respect to the calibration data (the observed stellar mass functions, their ratio, and the angular two-point correlation functions) is shown in Section 4.1. In Section 4.2 we investigate the resulting SHMRs for central galaxies (see Equation 5) and perform some further modeling of these with $\Lambda$CDM halo mass accretion histories. Satellite galaxies are explored in Section 4.3.

4.1 The Best-fit Model

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.

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). Whilst 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) $\sim 10^{-2}$ 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 ($\gtrsim 2\times 10^{-2}$ deg) for star-forming galaxies with $\log(M_{\star}/\mathrm{M}_{\odot})\lesssim 10.2$. Here, the discrepancy may be related to our decision to fix $\Delta_{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 $\Delta_{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.

In general, these discrepancies are fairly minor and the results shown here, combined with the reduced chi-squared value of $\chi^{2}/N_{\rm dof}\sim 1.9$ mentioned earlier, indicate that our model can reproduce the calibration data to a good degree of accuracy.

4.2 The Stellar-to-halo Mass Ratios

Figure 3 shows the best-fit SHMR, $\varepsilon(M_{\rm 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 $\sim 0.5$ dex) than the passive one (though at a similar stellar mass, $\sim 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_{\rm h}\lesssim 10^{13}$ M⊙*. We now investigate this difference in normalization further using a simple model comprising a realistic ensemble of halo histories.

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

[TABLE]

where $M_{0}$ is the halo mass at $z=0$. We generate a realistic ensemble of $\beta_{\rm h}$ and $\gamma_{\rm 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, $\psi(M_{\rm h},z)$, in each halo according to

[TABLE]

(e.g. Tacchella et al., 2018) where $f_{\rm b}=\Omega_{\rm b}/\Omega_{\rm m}$ is the baryon fraction and $\varepsilon^{\prime}(M_{\rm 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

[TABLE]

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_{\log M_{\star,\rm 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 $\varepsilon^{\prime}(M_{\rm h},z)$ and find that the parameter values ($\epsilon_{0}^{\prime}$, $M_{\rm c}^{\prime}$, $\beta^{\prime}$, $\gamma^{\prime}$) = ($0.15$, $2.1\times 10^{12}$ $h^{-1}$ 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 $\varepsilon^{\prime}(M_{\rm h},z)$ is likely more complex than investigated here, but we have illustrated how such differences can occur.

For each halo history, we now calculate a halo formation redshift, $z_{\rm 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 $\varepsilon^{\prime}(M_{\rm h},z)$ model in Figure 5, though we note that this Figure would be substantively similar if the constant $\varepsilon^{\prime}(M_{\rm h},z)$ model were considered. Interestingly, we find that 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.555We note that Feldmann et al. (2016) found a similar result studying passive galaxies in the FIRE simulations (see also Feldmann & Mayer 2015). 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.

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_{\star}\gtrsim 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.

4.3 Satellite Fractions and Quenching

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.

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.

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 passive666In reality, this fraction will be nonzero but less than it is at $z=2.3$ because galaxies were accreted at higher redshifts when the fraction of central galaxies that are passive will have been lower. (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 $40^{\rm th}$ percentile of the satellite age distribution, whereas in v1 it would be the $50^{\rm th}$ 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, $\propto(1+z)^{-3/2}$ [for a halo in virial equilibrium the dynamical time scales as $\rho^{-1/2}$ and the density, $\rho$, scales as $(1+z)^{3}$]. 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, $\eta$, 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_{\star}=10^{10.5}$ M*⊙* galaxies at $z\sim 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_{\star}=10^{10.5}$ M*⊙* we find $t_{\rm quench}=0.77^{+0.33}_{-0.22}$ Gyr for v1 and $1.03^{+0.60}_{-0.27}$ Gyr for v2.

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.

5 Summary

We have used a sophisticated halo occupation model to investigate the stellar-to-halo mass ratios (SHMRs) of passive and star-forming galaxies, as identified in the SMUVS survey, at $z\sim 2-3$. The free parameters in the model are calibrated through fitting them to the observed galaxy stellar mass and angular two-point correlation functions for each population. The model provides a meaningful way to statistically assign an average halo mass to an observed stellar mass and to divide our galaxy populations into central and satellite galaxies.

We find that the central galaxy SHMRs are different for passive and star-forming galaxies. Interestingly, the normalization of this relation around its peak is higher for passive galaxies than that of the star-forming population. Using a simple model based on $\Lambda$CDM halo mass accretion histories, we show that this can arise if the efficiency with which baryons are converted into stars at a given halo mass evolves with redshift, such that it is more efficient at higher redshifts. This mass accretion model also shows that passive central galaxies can be plausibly explained as simply a selection of the host halos with the highest formation redshifts, i.e., those with the lowest mass accretion rates. The halo formation redshift is very dependent on stellar mass for passive galaxies and this may help to explain the very different evolution of high- and intermediate-stellar mass selected passive galaxies found by Deshmukh et al. (2018).

At a given stellar mass, satellite galaxies occupy host halos of a greater mass than their central counterparts. Additionally, the satellite fraction of passive galaxies is greater than that of the combined population (passive and star-forming) indicating that environmental processes are expected to play a role in galaxy quenching.

Assuming that passive (i.e., quenched) satellites will be those that have been satellites for the longest and using the distribution of satellite (sub-halo) ages from an $N$-body simulation we derive satellite quenching timescales. Through comparing our estimates, and others from the literature at lower redshifts, to the results from very simple dynamical- and outflow-based models, we find that our data seem to favor a dynamically driven quenching model, e.g., ram pressure stripping. However, it is possible to draw different conclusions depending on the observational data considered, as they appear to support a range of possible scenarios. It appears that further studies of the satellite quenching timescale at higher redshifts will allow for a more definitive answer to this particular question.

Acknowledgements

The authors would like to thank the anonymous referee for a comprehensive and constructive report that allowed us to improve the overall quality of the manuscript. This work is based in part on observations carried out with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. It also uses data products from observations conducted with ESO Telescopes at the Paranal Observatory under ESO program ID 179.A-2005, data products produced by TERAPIX and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium, observations carried out by NASA/ESA Hubble Space Telescope obtained and archived at the Space Telescope Science Institute; and the Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. This research has made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA.

This work has made use of the following open-source python packages: numpy (van der Walt et al., 2011), scipy (Oliphant, 2007), matplotlib (Hunter, 2007) and ipython (Pérez & Granger, 2007). The color-vision-impaired-friendly color schemes used throughout can be found at https://personal.sron.nl/~pault/.

W.I.C., S.D., and K.I.C. acknowledge funding from the European Research Council through the award of the Consolidator grant ID 681627-BUILDUP. The Cosmic Dawn Center is funded by the Danish National Research Foundation.

Appendix A MCMC Fitting

In this Appendix, we give some brief details regarding the results of our MCMC fitting procedure. We use the Gelman & Rubin (1992) $\hat{R}$ statistic for each free parameter to assess whether our MCMC chains have converged. These are shown in the top right panel of Figure 9. All parameters have achieved values of $\hat{R}<1.2$, indicating that they have converged.

The marginalized $1$D and $2$D posterior distributions for our free parameters are shown in Figure 9. Whilst there are some degeneracies between parameters, this is not often the case, and the majority of the $1$D posterior distributions form well-defined peaks. Our best-fit parameter values and $16-84^{\rm th}$ percentile ranges from our MCMC chains are given in Table 1.