Open Access
Issue
A&A
Volume 664, August 2022
Article Number A61
Number of page(s) 28
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202243136
Published online 05 August 2022

© M. Shuntov et al. 2022

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe-to-Open model. Subscribe to A&A to support open access publication.

1. Introduction

Within the current paradigm of structure and galaxy formation, galaxies form and evolve within dark matter halos (White & Rees 1978). Their properties are inextricably connected in what is known as the galaxy-halo connection (see e.g., Wechsler & Tinker 2018, for a review). One facet of the galaxy-halo connection is the relationship between the mass of the dark matter halo and the stellar mass of the galaxy it hosts, referred to as the stellar-to-halo mass relation (SHMR). The SHMR expresses the efficiency of the stellar mass assembly of a galaxy integrated over the halo’s lifetime and it is, at first-order, a function of the halo mass and is shaped by the physical mechanisms of galaxy formation (e.g., Somerville & Davé 2015, for a review).

Galaxy formation is an inefficient process: the ratio of stellar mass to halo mass, M*/Mh, is quite low (Shankar et al. 2006; Mandelbaum et al. 2006; Zheng et al. 2007; Conroy & Wechsler 2009; Behroozi et al. 2010). This quantity is a strong function of halo mass and rises to a peak at a characteristic peak halo mass, suggesting that at most only 20% of all the available baryons in the halo have turned into stars. At lower and higher halo masses, the M*/Mh ratio decreases rapidly, which is seen as a signature of different feedback processes that suppress star formation and act at different halo mass scales: stellar feedback in low-mass halos and active galactic nuclei (AGN) in high-mass halos (e.g., Silk & Mamon 2012, for a review).

Central and satellite galaxies contribute to the total stellar mass content of halos. As a consequence, the total SHMR can be decomposed in the contributions from both. In lower mass halos, the central galaxy makes up most of the stellar mass content, but its growth is regulated by stellar feedback mechanisms such as supernovae (SNe), stellar winds, radiation pressure, and photoheating. All of these mechanisms are important in explaining the observed amount of stellar mass in lower mass galaxies; otherwise, the galaxy masses end up overpredicted (Hopkins et al. 2012). On the other hand, in cluster-scale halos, the satellite galaxies dominate the stellar mass budget, mostly due to their high number (Leauthaud et al. 2012; Coupon et al. 2015). Within large halos, the stellar mass assembly in satellites is described by interplay between the hierarchical merger tree, where smaller halos accrete into larger ones and in situ star formation facilitated by cold gas inflows and the various environmental quenching mechanisms that act to suppress further growth (Peng et al. 2010; Gabor & Davé 2015). These include the so called “hot-halo”, where the infalling gas is heated by virial shock heating (Birnboim & Dekel 2003), as well as mechanisms such as strangulation (Larson et al. 1980; Balogh et al. 2000), ram-pressure stripping (Gunn & Gott 1972), and harassment (Moore et al. 1996). In such massive halos, AGN feedback keeps the gas hot through the so-called “radio” mode and is necessary for explaining the break in the local galaxy luminosity/stellar mass function at the bright/massive end (Croton et al. 2006; Bower et al. 2006). All of these factors shape the total (central + satellite) SHMR, and its evolution with redshift can indicate the relative importance of these processes in shaping the star-formation efficiency as a function of the halo mass.

Some of the most predictive models that give direct insight into the physical processes that shape the galaxy-halo connection come from hydrodynamical simulations (e.g., Vogelsberger et al. 2014; Dubois et al. 2014; Schaye et al. 2015; Pillepich et al. 2018a). However, the known physics implemented in the simulations imposes a strong prior on the galaxy evolution and cannot provide information about physical processes that had not been previously expected to contribute to the galaxy-halo connection (Behroozi & Silk 2018). Additionally, an important caveat of simulations is that they cannot simulate the full physics at all scales, so they rely on various parametrizations below the resolution scale, namely, subgrid physics. Subgrid physics varies from one simulation to another, and so do its conclusions on the galaxy-halo connection. Even though hydrodynamical simulations are the most predictive, they are computationally expensive, which makes it difficult for them to be constrained against observations using techniques such as Markov chain Monte Carlo (MCMC).

Empirical models offer the most flexibility (Wechsler & Tinker 2018). Many methods have been described in the literature to compute the SHMR from observational data sets or numerical simulations. One of the most widely employed techniques is abundance matching (AM), where the abundance (i.e., number density in a comoving volume) of galaxies above a given mass is matched to the abundance of dark matter (sub)halos which then gives the halo mass (Marinoni & Hudson 2002; Kravtsov et al. 2004; Vale & Ostriker 2004; Tasitsiomi et al. 2004). Another method based on a statistical description of the way galaxies populate halos is the halo occupation distribution (HOD) model. The HOD can model one- and two-point statistical observables such as the galaxy stellar mass function (GSMF) and galaxy clustering as measured by the correlation function (2PCF; e.g., Peacock & Smith 2000; Seljak 2000; Scoccimarro et al. 2001; Berlind & Weinberg 2002). Within the HOD formalism, we can also model observables that directly probe dark matter halos, such as galaxy-galaxy lensing (Leauthaud et al. 2011). In galaxy-galaxy lensing, one measures the subtle coherent distortions of the shapes of background galaxies induced by the foreground matter distribution to probe the latter. However, it relies on accurate measurements of galaxy ellipticities, which becomes increasingly difficult at z >  1 and poses a limitation with regard to its implementation in probing a redshift evolution over a vast interval.

The phenomenological approaches of AM and HOD are based on statistical measures (e.g., GSMF and 2PCF) that indirectly probe the dark matter halos and are agnostic with regard to the physical processes that shape their relation with the galaxies. An additional drawback of these phenomenological models is that they rely on accurate knowledge of the halo mass function, which has to be calibrated using numerical simulations; in addition, they are sensitive to various definitions of the halo profile, radius, mass, concentration, and bias. However, this type of modeling is not constrained by strong priors from the physics of galaxy evolution – it is almost entirely constrained by observations and, as such, they can reveal signatures of new physical processes that shape galaxy properties.

This paper aims to constrain the redshift evolution of the SHMR up to z ∼ 5 by applying an HOD-based analysis consistently on a homogeneous data set: the COSMOS2020 photometric catalog. COSMOS2020 (Weaver et al. 2022) is the latest iteration of the photometric catalog of the COSMOS (Scoville et al. 2007) survey that includes the latest data-releases of deep imaging, covering wavelengths from the ultraviolet to the near-infrared. The deep multi-band photometry allows for the estimation of accurate photometric redshifts and stellar masses, along with a selection of complete samples up to high redshift. By adopting an HOD-based model to jointly fit for galaxy abundance and clustering, our analysis is aimed at constraining the satellite contribution to the total stellar mass budget in halos across a vast redshift range. This allows us to infer a coherent picture of how the stellar mass assembles as a function of the halo mass throughout cosmic history.

The novel aspect of our work is comprised by the use of a single dataset to perform all the measurements and probe the SHMR to z ∼ 5 for both central and satellite galaxies. Most of the investigations in the literature have relied on observables from heterogeneous data sets to constrain their models (e.g., Behroozi et al. 2013, 2019; Moster et al. 2018). Different data sets can have different selection functions and methods of estimating galaxies’ physical parameters that have the capacity to propagate various systematic biases, which may muddle the interpretations (Behroozi et al. 2010). Therefore, our work is free from such “inter-observational” systematics.

Our work builds up on the literature in several ways. Legrand et al. (2019) is the only work that has measured the SHMR using a single data set (COSMOS2015 of Laigle et al. 2016) up to z ∼ 5 using sub-halo abundance matching. A shortcoming of this approach is the satellites are treated as centrals in their own sub-halo so that they only predict the SHMR of central galaxies. Our approach allows for centrals and satellites to be decoupled in order to compute the contribution of both to the total mass content of halos. Previous works that have measured both central and satellite SHMR are limited only to z <  1 (e.g., Leauthaud et al. 2012; Coupon et al. 2015), or a single z-bin measurement at 2 <  z <  3 as in Cowley et al. (2019). Therefore, this paper presents the only measurement of the SHMR for both centrals and satellites up to z ∼ 5 using a homogeneous dataset: COSMOS2020.

The organization of this paper is as follows. In Sect. 2, we describe the COSMOS dataset applied here and the mass-selected samples in the ten redshift bins comprised by our analysis. In Sect. 3, we present the methods we employ to perform our measurements of galaxy abundance and clustering. In Sect. 4, we lay out the HOD-based modeling of our observables with a parametrization of the SHMR as a starting point. In Sect. 5, we present our measurements of the observables and the results on the redshift evolution of the HOD and SHMR. In Sect. 6, we discuss the physical mechanisms that may regulate the growth of central and satellite galaxies in dark matter halos. We also compare our results with hydrodynamical simulations and discuss the possible origins of the discrepancies that we find. Finally, we summarize our finding in Sect. 7

Throughout this paper we adopt a standard ΛCDM cosmology with H0 = 70 km s−1 Mpc−1, Ωm, 0 = 0.3, where Ωb, 0 = 0.04, ΩΛ, 0 = 0.7, σ8 = 0.82 and ns = 0.97. Galaxy stellar masses, when derived from spectral energy distribution (SED) fitting, scale as the square of the luminosity distance (i.e., ), therefore, as h−2; dark matter halo masses, usually derived from dynamics in numerical simulations, scale as h−1. The h scaling factors are retained implicitly for all relevant measurements, unless explicitly noted otherwise (see Croton 2013, for an overview of h and best practices). When making comparisons to the literature, we rescale all the measurements to the cosmology adopted for this paper. All magnitudes are expressed in the AB system (Oke & Gunn 1983). Stellar masses are obtained assuming Chabrier (2003) initial mass functions (IMF) and when comparing to the literature, stellar masses are rescaled to match the IMF adopted in this paper.

2. Data

2.1. COSMOS2020 catalog

This work makes use of the COSMOS2020 catalog (Weaver et al. 2022). This deep multi-wavelength near-infrared selected catalog uses deep observations over the 2 deg2 COSMOS field in 35 photometric bands from ultraviolet (UV) to near-infrared (NIR). This unique combination of depth, area, and wavelength coverage allows an accurate estimation of photometry, photometric redshifts, and stellar masses for around a million sources up to z ∼ 10.

Briefly, COSMOS2020 comprises two photometric catalogs extracted on a “chi-squared” combination (Szalay et al. 1999) of deep near-IR images in izYJHKs (AB mag 3σ depth in 2″ apertures of 27.6,  27.2,  25.3,  25.2,  24.9,  25.3). One catalog uses the traditional approach of measuring fluxes in fixed apertures using SExtractor (Bertin & Arnouts 1996), while the second uses a profile-fitting technique using The Farmer (Weaver et al., in prep.) built around The Tractor. Photometric redshifts (photo-zs), stellar masses, and other physical parameters are estimated for these two photometric catalogs using two spectral energy distribution (SED) fitting codes LePhare (Arnouts et al. 2002; Ilbert et al. 2006) and EAZY (Brammer et al. 2008).

In this analysis, we use the CLASSIC catalog with photometric redshifts estimated using LePhare. This is principally because photometric measurements carried out with The Farmer, which is based on a model fitting technique, can fail to converge in certain cases, especially in crowded regions or near bright sources. This spatially variable completeness is problematic for measurements of the angular correlation function; we have measured relative correlation function differences of up to about 30% on scales of 1′ between the two catalogs. The CLASSIC catalog, on the other hand, contains photometric measurements for almost all sources within the survey masks. There is a caveat, however, since aperture photometry is not very reliable in crowder regions and around bright sources neither. So whereas FARMER is a pure catalog, since it photometers all the reliable sources, CLASSIC is a more complete catalog, photometrying almost all sources.

Compared to COSMOS2015, we know that COSMOS2020 reaches similar photometric redshift precision almost one magnitude fainter – this is shown in Fig. 17 of Weaver et al. (2022), where the 1σ uncertainty of the photo-zs is plotted as a function of redshift and magnitude bin. The normalized median absolute deviation (σNMAD1) at i <  22.5 is below 0.01 (1 + z) and stays below 0.05 (1 + z) to 25 <  i <  27. The outlier fraction2 is below 1% and 20% for the corresponding magnitude bins. The bias3 ranges from −0.003 to −0.014 in the bright and faint magnitude bins, respectively.

Accurate and complete stellar mass estimates over a wide redshift range are necessary for obtaining an accurate measurement of the SHMR. In COSMOS2020, this is enabled by the inclusion of the deep near-IR data from the UltraVISTA survey (McCracken et al. 2012), DR4 in YJHKs, and mid-IR data from the Cosmic Dawn Survey (Moneti et al. 2022), as well as Spitzer/IRAC observations in channels 1–2 (3.6 μm, 4.5 μm). Stellar masses are estimated with LePhare using SED templates produced from stellar population synthesis models by Bruzual & Charlot (2003) and initial mass functions by Chabrier (2003). The SEDs are fixed at z = zphot then and fitted to the multi-wavelength photometry (for more details, see Weaver et al. 2022 and Laigle et al. 2016). The point estimate of the stellar mass is the median of the resulting PDF marginalized over all other parameters, with the 16th and 84th percentiles of the PDF giving the 1σ confidence interval. The improved depth (e.g., Ks = 25.3, [3.6]=26.4 at 3σ) translates to higher stellar mass completeness compared to the previous versions of the catalog. This enables a selection of samples based on stellar mass that are complete down to log M*/M ∼ 8.2 at z ∼ 0.3 and log M*/M ∼ 9.3 at z ∼ 4.

Throughout this work, we used 2″ aperture magnitudes and apply aperture-to-total and Milky Way extinction corrections using the Schlafly & Finkbeiner (2011) dust map. We applied masks to remove sources near bright stars and in regions near artefacts. This leaves an effective area of 1.27 deg2 that corresponds to the footprint of the UltraVISTA survey4.

Finally, we use two star-galaxy classifications to remove uncorrelated stars from the catalog. One uses morphological information from HST/ACS and Subaru/HSC images, where half-light radii and magnitudes classify as stars all point-like sources at i <  23 and i <  21.5 in ACS and HSC images, respectively. This criterion is also satisfied by point-like AGN sources. The second, which is an SED-based criterion, classifies as stars those sources with the χ2 of the best-fit stellar template lower than the χ2 of the best-fit galaxy template (for more details see Weaver et al. 2022). We performed tests by measuring the correlation function of sources classified as stars, while further removing point-like AGN sources based on their χ2. The correlation function of our stellar sample is zero, indicating a clean separation.

2.2. Sample selection

Measuring galaxy clustering and abundance requires complete stellar mass-selected and volume-limited samples. To select stellar mass-complete samples we use the stellar mass completeness limit, which is computed following the method prescribed by Pozzetti et al. (2010) and is described in catalog paper of Weaver et al. (2022). To ensure complete samples, the channel 1 limiting magnitude [3.6]lim = 26 is computed with the help of the deeper CANDELS-COSMOS catalog (Nayyeri et al. 2017) which is used for completeness check. All samples throughout this work are selected to be brighter than [3.6]lim = 26 across the full redshift range, despite the fact that a Ks based selection is also suitable at low redshifts.

To probe the cosmic evolution of the observables we bin the samples in ten redshift bins from z = 0 to z = 5.5 with varying widths. The widths were chosen to ensure roughly the same number of galaxies in each bin. The redshift bins are listed in Table 1. Additionally, we require that each galaxy has its lower and upper 1σ values (zlow, zup) within ±0.5 of the z-bin limits. This criterion removes any galaxies with highly uncertain redshifts that can introduce errors. To ensure that the samples remain as complete as possible, we don’t impose any other selection criteria, for example, based on S/N, number of bands in which a source was photometered, χ2 of the SED fit etc.

Table 1.

Sample selection in redshift and stellar mass thresholds.

One of the ingredients of the model of the galaxy correlation function is the redshift distribution N(z). We used the z-likelihood from LePhare to build N(z). Formally, for each galaxy, there is a likelihood of the observed photometry (denoted by the vector of fluxes o) given the redshift ℒ(o|z). N(z) is then constructed by simply stacking the individual ℒ(o|z) in each z-bin:

(1)

where the sum is over the number of objects of the redshift and stellar mass threshold-selected sample Nsample.

Constructing the N(z) in this way directly accounts for the uncertainty in the photometric redshift when selecting galaxies in a bin, namely, the fact that galaxies can still have their true z outside the bin limits. However, it has been shown by Ilbert et al. (2021) that such a procedure can lead to biases in the mean redshift that can be inferred from the N(z). They have quantified the notion that these biases can reach up to ∼0.01 (1 + z). By using the model of w(θ) (described in Sect. 4.6) we tested the effects of bias in the mean redshift, as well as different estimates of N(z). Our conclusion is that biases on an order of magnitude of ∼0.01 (1 + z) result in relative differences in the value of w(θ) of less than about 3%, which is considerably smaller than the typical relative error of the measurement (about 10%). On the other hand, the shape of N(z) (notably, the width of its wings) can lead to significant differences in w(θ) that can bias the inferred SHMR parameters. The reason for this is the mix of physical scales when considering larger volumes: the angular correlation of galaxies selected in a wider radial interval is inevitably lower, since they can be far apart in the radial direction but close in angular separation. For example, considering N(z) to be a Gaussian distribution centered at the z-bin mean and with width half of the z-bin width, can lead to a relative difference of about 20% at z ∼ 0.4 and more than 50% at z ∼ 2.5 (see Appendix A). Ilbert et al. (2021) have shown that the ℒi(o|z) as output from LePhare can lead to biased and miss-calibrated N(z) as evidenced by comparison with the true redshift histogram and the Probability Integral Transform (PIT) statistic (see Fig. 4 in Ilbert et al. 2021). The authors show that an N(z) that is better representative of the true distribution can be obtained by using a posterior distribution, such as:

(2)

where Pr(z|m0) is the so called ‘photo-z prior’ (Brodwin et al. 2006). This prior can be constructed for every z-bin by summing the likelihoods per magnitude bins such as:

(3)

where Θ(m0, i|m0) is equal to 1 if the object’s magnitude m0, i is within the magnitude bin centred at m0 and zero otherwise. The outcome of this procedure and the effects on w(θ) are presented in more detail in Appendix A. We adopt the N(z) obtained using Eq. (2) for our analyses. The N(z) is constructed for each considered sample including the mass threshold-selected samples. The resulting distributions for all z-bins for galaxies above the mass-completeness limit are shown in Fig. 2. We note that there are some narrow (of a width of ∼0.01) dips at several z-values (e.g., z ∼ 1.3,  2.9,  4.0). These come from the individual likelihoods being close to zero at these exact values. The template fitting outputs minimal likelihoods at exactly these three z-points, and as they are narrow (∼0.01), they do not affect the selection nor the modeling.

To probe the clustering strength as a function of galaxy mass, we further select samples in stellar mass thresholds. In each z-bin we define samples selected in several stellar mass thresholds starting from the mass completeness limit; these are indicated by the horizontal solid lines in Fig. 1 and listed in Table 1. We chose the thresholds rather arbitrarily to ensure a good signal-to-noise (S/N) for the clustering measurement of each mass threshold-selected sample within a z-bin.

thumbnail Fig. 1.

Sample selection in the stellar mass-redshift plane. The solid grid lines show the mass threshold for each sample. The solid violet curve indicates the stellar mass completeness limit. We note that the histogram does not correspond to the final selection since further selection criteria (e.g., based on PDF(z) width) are applied.

thumbnail Fig. 2.

Redshift distribution of the ten galaxy samples used for the clustering and abundance measurements. The redshift distribution is obtained by stacking the posterior photo-z distributions for all the sources in a given bin as described in Sect. 2.2.

3. Measurements

3.1. Galaxy clustering

We measured the two-point angular correlation function w(θ) using the Landy & Szalay (1993) estimator5

(4)

where DD is the number of data-data pairs in a given angular separation bin [θ, θ + δθ], RR is the number of pairs in the random catalog in the same bin, and DR is the number of pairs between the data and the random catalog. The data and random pairs are normalized by the total number of galaxies and random objects. We constructed the random catalog with the same survey geometry mask as the data catalog and Nrandom ∼ 3 × 106, which is more than 50 times the number of galaxies in each considered bin.

The covariance matrix is computed using the jackknife method by subdividing the full area in Npatch = 22 patches and recomputing the correlation function removing one patch at a time. The covariance matrix is then estimated as:

(5)

where is the mean correlation function and wi is the correlation function with the ith patch removed. The final covariance matrix thus includes uncertainties due to sample (cosmic) variance, dominating at larger scales, and due to Poisson statistics from counting objects in bins dominating mostly at small scales. However, due to the limited survey size, this method still cannot accurately estimate the uncertainties due to cosmic variance. To capture the cosmic variance effects in the covariance matrix one can use computationally expensive simulations, which is out of the scope of this work. This, for example, is done in Leauthaud et al. (2011), where they show that cosmic variance has an impact on the covariance matrix on large scales (i.e., in the two-halo regime), namely, as the cosmic variance increases the correlation in the data at large scales.

Due to the limited size of the survey (1.27 deg2), the clustering measurements suffer from the effects of the integral constraint (IC, Groth & Peebles 1977). This leads to an underestimation of w(θ) at large scales comparable to angular size of the survey by a constant factor wIC, such that the true correlation function is:

(6)

We incorporate this correction into our model, which is described in detail later in this work.

3.2. Galaxy abundance

We measured stellar mass functions across ten redshift bins using the 1/Vmax technique (Schmidt 1968). This estimator weighs each galaxy by the maximum volume in which it would be observed given the redshift range of the sample and magnitude limit of the survey. The comoving volume Vmax is computed between zmin and zmax, where zmin is the lower redshift limit of the z-bin and zmax = min(zbin, up, zlim), then zbin, up is the upper redshift limit of z-bin, and zlim is the maximum redshift up to which a galaxy of a given magnitude can be observed given the magnitude limit of the survey. For this purpose, we used IRAC channel-1 magnitudes and a limit of [3.6]=26. We computed the SMFs in the mass range starting from the mass completeness limit of each z-bin, with a bin size of Δlog M = 0.25.

Uncertainties in the SMF include contributions from Poisson noise (σPois), cosmic variance (σcv), and SED fitting uncertainties (σfit). We computed the uncertainties due to cosmic variance following Steinhardt et al. (2021). The starting point is the cosmic variance “cookbook” code of Moster et al. (2011), that computes the stellar-mass dependent cosmic variance and performs well up to intermediate redshifts and masses but becomes increasingly underestimated at high redshift and mass. Steinhardt et al. (2021) extend the recipes to the early universe (z >  3) by using linear perturbation theory.

Additionally, due to photometric errors and degeneracies in the SED fits, there are uncertainties in the M* measurements that propagate to the SMF. To estimate SED fitting uncertainties on the SMF we use the PDF(M*). We assign to each galaxy a weight that corresponds to its probability to be found in the given mass bin:

(7)

We then used this weight to compute the combined Poisson and SED fitting uncertainties in the following way:

(8)

The final uncertainty is then . The measurements for the clustering and SMF are presented in Fig. 4 and discussed in Sect. 5.2

4. Theory and modeling

The observables of galaxy clustering and abundance can be predicted within the framework of the halo model. The halo model of the large scale structures postulates that all matter in the universe is contained in virialized dark matter halos. Using halos as the basic unit, it provides a method to statistically describe the distribution of dark matter and analytically compute its clustering. Combined with a model that describes how galaxies populate halos – including both central and satellite galaxies, one can model various statistics of the galaxy distribution. One of the principal advantages of this modeling is the capability to infer a range of properties about the satellite galaxies, which inferred in such a large redshift span is one of the novelties of our work. The power in the clustering of all galaxies on small scales and the number densities at the low mass end of the GSMF are both shaped by satellites. Therefore, by fitting the observed correlation function and GSMF with models that consider the satellites, one can infer the parameters that govern the statistical distribution of satellites within dark matter halos. We note that, the observations don’t distinguish centrals from satellites.

4.1. Stellar-to-halo mass relation

The stellar mass assembled by a galaxy depends most strongly on the mass assembly of the host halo. Consequently, a strong relationship between them is evident – the stellar-to-halo mass relation (SHMR). Our goal is to constrain the SHMR. Leauthaud et al. (2011), hereafter L11, laid the theoretical framework to model galaxy clustering and SMFs based on the HOD formalism with as a starting point a functional form of the SHMR. This function must capture the different growth rates of galaxies as a function of the halo mass that is shaped by various feedback processes that operate at different mass scales. Behroozi et al. (2010) presented a functional form of the SHMR, which we have adopted for our analysis:

(9)

This relation is parameterized by a characteristic halo mass and stellar mass scales given by the parameters M1 and M*, 0, respectively. Here, M1 controls the normalization of Mh as a function of M* and M*, 0 controls the position along the M*-axis. Together, two parameters govern the transition mass scale between the low-mass and high-mass regime of the SHMR. The low-mass regime (M* ≲ 1010.5M) is described by a power-law regulated by the parameter β. The high-mass regime follows a sub-exponential law regulated by the parameter δ, and the transition regime is shaped by the parameter γ.

The ratio between the stellar mass and the halo mass obtained with Eq. (9), (i.e., M*/Mh) can be considered as the efficiency of the galaxy formation process that encapsulates all the processes that lead to the conversion of baryons to stars (from gas cooling and star formation to stellar and AGN feedback); this can be referred to as the baryon conversion efficiency or star formation efficiency (SFE). Since we can consider that baryonic matter content of halos is equal to the universal baryonic fraction fb = Ωbm ≈ 0.16, the M*/Mh ratio will inform us of the fraction of the baryons available in the dark matter (DM) halo that have converted into stars. At a given redshift, the M*/Mh ratio gives the baryon conversion efficiency integrated over the lifetime of the halo and therefore includes the combination of all the different physical processes that regulate star formation throughout the halo life (e.g., gas accretion, mergers, feedback). The shape of the Mh/M* ratio is a strong function of halo mass, which indicates that various feedback mechanisms operate on different halo mass scales to regulate star-formation. This is the principal quantity of interest in this work and we extensively discuss it in Sect. 5.

As L11, we assume that this SHMR concerns only central galaxies occupying the central regions of the dark matter halos. These halos can contain smaller sub-halos that orbit the potential well and can host satellite galaxies. Satellites can undergo different stellar growth than centrals, since various distinct processes affect satellites in the dark matter halos that can regulate their growth (such as stripping and harassment). This will reflect in the SHMR, therefore a different relation for satellites is a more accurate assumption. This is one of the advantages of this methodology to infer the SHMR, as compared to the commonly implemented abundance matching (AM) technique. In AM, one treats satellites as centrals in their own sub-halo, which then assumes that centrals and satellites follow the same SHMR. The formalism that we employ is able to constrain the contribution of both centrals and satellites in the total stellar mass content in halos of given mass (Sect. 4.4).

Due to the effects of the various galaxy formation processes, there is a scatter of galaxy stellar mass that exists (as well as other galaxy properties) at a fixed halo mass. This stochastic nature of the SHMR can be modeled with a conditional function that describes the probability of observing a central galaxy with M* at a given Mh, which can be chosen to be a log-normal distribution. The dispersion of the log-normal distribution σLog M* describes the scatter in stellar mass at fixed halo mass and is a free parameter that can be fitted with the data. Following previous works (e.g., Moster et al. 2010; Leauthaud et al. 2012; Tinker et al. 2013), we consider σLog M* to be independent of the halo mass, and we leave it as a free parameter to be fitted in each z-bin. Hydrodynamical simulations, however, show that at z = 0, σLog M* generally decreases with Mh going from ∼0.32 at Mh ∼ 1011 to ∼0.15 at Mh ∼ 1012 and staying constant to higher masses (Pillepich et al. 2018a). We checked that varying σLog M* over this range changes the clustering correlation function insignificantly (and within the measurements errorbars) and given the fact that this parameter is mostly constrained by the GSMF, taking σLog M* independent of the halo mass is a safe assumption for our purposes.

4.2. Central occupation distribution

The HOD describes the statistical occupation of galaxies in dark matter halos. It assumes a probability distribution of the number of galaxies residing in halos conditioned on some criteria, usually on the mass P(N|Mh). Typically, centrals are assumed to follow a Bernoulli distribution, while the number of satellites follows a Poisson distribution (see e.g., Zheng et al. 2005, and references therein). Under these assumptions, the HOD is described by the average number of galaxies with stellar masses higher than some threshold in halos of a given mass ,

(10)

is a monotonic function increasing from 0 to 1. fSHMR(Mh), whose inverse function is defined by Eq. (9) and gives the stellar mass at the halo mass, while σLog M* is the scatter in stellar mass at fixed halo mass. We note that all the other parameters regulating the central HOD parametrize the functional form of the SHMR. An alternative approach employed by many studies (e.g., Zheng et al. 2007; Zehavi et al. 2011; Coupon et al. 2012; McCracken et al. 2015; Ishikawa et al. 2020) specifies the central HOD assuming the SHMR to be a simple power law, with a parameter quantifying the minimum halo mass to host a galaxy Mmin. The downside of this model is the difficulty in the interpretation of Mmin as the halo mass at the stellar mass threshold, that is, , especially at high masses where the deviation from a power-law of the SHMR is clear (see L11 for more details).

4.3. Satellite occupation distribution

The occupation of halos by satellites can be modeled by a power-law at high halo masses with an exponential cut-off at low masses, given by:

(11)

where αsat is the power-law slope, Msat is the halo mass scale for the satellites defining the amplitude of the power law and Mcut is the cutoff scale. HOD studies have shown that the satellite mass scale is proportional to at the threshold stellar mass (e.g. Zheng et al. 2007; Zehavi et al. 2011). This allows us to parametrize Msat and Mcut as power laws by introducing four more parameters:

(12)

The HOD is fully specified by the average occupation number of galaxies in halos, as given by Eqs. (10) and (11). Finally, the total number of galaxies including centrals and satellites is simply ⟨Ntot⟩=⟨Ncent⟩+⟨Nsat⟩. We can also compute the average number of galaxies in a mass bin of by simply taking the difference

(13)

The model has a total of 11 parameters. The SHMR for the centrals has five parameters (M1, M*, 0, β, δ, γ) with one additional parameter that describes the scatter in stellar mass at a fixed halo mass of σLog M*. The occupation distribution for satellites is modeled with five parameters (αsat, Bsat, βsat, Bcut, βcut). We did not parametrize the redshift evolution of these parameters; instead we inferred them for the redshift bins from 0.2 <  z <  5.5 and then looked for their evolution with a value of z determined a posteriori.

4.4. Total stellar content in halos

From the model of the conditional mass function, it is possible to compute the total stellar mass contained in halos of a given mass by performing an integration over the stellar mass. Since we do not have a model of Φs(M*|Mh), we can use the occupation distributions of centrals and satellites because they are also integrals of the conditional mass function. Therefore, the contribution of centrals and satellites to the total stellar mass content in halos can be computed as:

(14)

This equation computes the contribution of galaxies (centrals and satellites) in a stellar mass bin to the total stellar mass content in halos of Mh. The total SHMR then shows the overall efficiency of the galaxy formation process in halos, integrated over the halo’s history, that is a combination of the in situ conversion of gas to stars and ex-situ from merging with satellites.

4.5. Model of the galaxy stellar mass function

From the defined occupation distribution of halos (Eqs. (10) and (11)), we can obtain the number density of galaxies in a given mass bin, , by integrating over the halo mass function (HMF) dn/dMh:

(15)

This allows us to also compute the GSMF of centrals and satellites by using their respective occupation distributions (shown in Eq. (15)). The literature abounds with prescriptions of the HMF obtained under various assumptions and methods, and for our work we applied the HMF of Despali et al. (2016). The adoption of different HMFs has an effect on the modeled GSMF and inevitably on the inferred model parameters. The HMF also depends on the choice of halo mass definition, and since the models are computed using the HMF, the final results will also depend on these definitions. For example, a different halo mass definition would result in a systematic shift in halo masses. To compute the HMF, we used the COLOSSUS code (Diemer 2018) and we used the virial overdensity (Bryan & Norman 1998) halo mass definition for the results we present in this paper. The model for the SMF shows a high sensitivity to the parameters describing the central SHMR, and coupled with the high signal-to-noise of the measurements has the most constraining power.

4.6. Model of the two-point angular correlation function

The model of the two-pt angular correlation function follows closely the usual prescriptions (see e.g., Cooray & Sheth 2002). For completeness, we detail the principal equations in Appendix B. For the computation of the two-pt angular correlation function, we rely on LSST Dark Energy Science Collaboration’s Core Cosmology Library (CCL)6. It offers a library of routines to calculate a range of cosmological observables and is still under active development. The validation of the software along with a range of benchmark tests are presented in Chisari et al. (2019). The main ingredients that enter the modeling along with the prescriptions and assumptions we adopt here are given in Table 2.

Table 2.

Adopted ingredients in the halo model.

Due to the relatively small volume probed by the COSMOS survey, the integral constraint affects w(θ) at large scales. We adjust the model to take this into account. The correction factor due to the IC can be estimated from the double integration of the true correlation function over the survey area:

(16)

This integration can be carried out using the random-random pairs from the random catalog following Roche & Eales (1999)

(17)

where wtrue(θ) is HOD-predicted model. Finally, the model that we fit against the data is simply w(θ)=wtrue(θ)−wIC.

5. Results and analysis

5.1. Fitting procedure

We fit the models of the w(θ) and the SMF to our measurements using a Markov chain Monte Carlo (MCMC) approach, minimizing χ2 as:

(18)

where w represents the measurement vector containing w at θ for all mass thresholds, while and are the models for a given set of parameter values. The first line of Eq. (18) corresponds to the clustering likelihood and the second line to GSMF likelihood. We use the affine-invariant ensemble sampler implemented in the emcee code (Foreman-Mackey et al. 2013). We used 200 walkers for our 11 parameters and relied on the auto-correlation time τ to assess the convergence of the chain. To consider the chains converged, we require that the auto-correlation time is at least 60 times the length of the chain and that the change in τ is less than 5%. We discarded the first 2 × max(τ) points of the chain as the burn-in phase and thin the resulting chain by 0.5 × min(τ). We imposed flat priors on all parameters; for the mass parameters, the flat priors are on the log quantities.

For the best-fit parameters values, we take the medians of the resulting posterior distribution, with the 16th and 84th percentiles giving the lower and upper uncertainty estimates. The best-fit parameters and their uncertainties for all the 10 z-bin are listed in the Appendix F. The posterior distributions for the 11 parameters in the redshift bins are shown in Appendix E.

5.2. Measurements and best-fit models

The measurements of the w(θ) and GSMFs are shown in Fig. 3, where we isolate the measurement in 0.5 <  z <  0.8 and compare it with clustering measurements from the literature, as well as in Fig. 4, where we show the measurements in all the other z-bins. In the upper panel, we show the clustering measurements (open circles with errorbars) for different mass-threshold samples and in the bottom panel the SMF measurement. The solid lines show the best-fit models with the color code corresponding to the mass-threshold measurement (in the top panel). Table F.1 also shows the reduced chi-squared value for the-best fit parameters. Their values range from 2.5–6 for most bins except for 0.8 <  z <  1.1, where . Given the number of data points, ranging from 75 at 0.8 <  z <  1.1 to 30 at 4.5 <  z <  5.5, these values of indicate a reasonably good fit. One possible explanation can be the greater complexity of the data, which might not be completely captured by the fits. For example, we are simultaneously fitting for several mass-selected clustering measurements. Due to uncertainties in the stellar masses, there some super-covariance may exist between all the measurements, however, this type of modeling is beyond the scope of this work.

thumbnail Fig. 3.

Best-fit models for clustering and abundance compared with measurements at 0.5 <  z <  0.8. Top: clustering measurements for the 6 mass threshold-selected samples (empty circles with errorbars) along with the best-fit models in solid lines in corresponding colors. Bottom: measurements of the stellar mass function together with the best fit model. The dashed lines show the SMF in the same redshift bin obtained by Davidzon et al. (2017) for comparison.

thumbnail Fig. 4.

Best-fit models of clustering and abundance plotted over the measurements for all z-bins apart from 0.5 <  z <  0.8, which is shown in Fig. 3.

Description of the clustering measurements. The clustering measurements exhibit the usual behavior, with w(θ) following a power-law at small scales that breaks at intermediate scales (∼1′). The origin of this break comes from the fact that the power at small scales is dominated by galaxies residing in the same halo (1-halo term), which drops off quickly at intermediate scales; at larger scales, the power mainly comes from large-scale clustering of halos (two-halo term) and the transition between these two regime creates the characteristic shape (Zehavi et al. 2004). The relative contribution from these two terms becomes more apparent at high masses and high redshifts, with the one-halo term dominating the power at small scales with a steep slope and the two-halo term dominating the large scales with a shallower slope at the transition. The clustering amplitude increases with increasing mass threshold – a familiar behavior based on the fact that massive galaxies trace high-density and more clustered regions.

At scales larger than 0.1 deg, there is a sharp drop in power due to the effects of the integral constraint. The best-fit models shown in solid lines generally agree well with the measurements. It should be noted that the measurements show an excess of power at scales of ≳0.02 deg (∼0.5 Mpc at z ∼ 7) and the fits are systematically below the data points. This may be due to the effects of the non-linear halo bias effect (scale-dependent halo bias). While in this work we use the scale-independent halo bias of Tinker et al. (2010), some studies suggest that the halo bias is scale dependent in the quasi-linear regime at scales of about 1 Mpc (Jose et al. 2017). This will add power in the correlation at scales of ∼0.04 deg. Furthermore, at 0.5 <  z <  1.5, a contribution may also come from the known overabundance of rich structures in the COSMOS field at these redshifts, as discussed by McCracken et al. (2007, 2015) and Meneux et al. (2009). This excess of power will decrease the SHMR, indicating an even lower efficiency of converting baryons to stars. Clustering is particularly sensitive to the satellite content within the halo, therefore the parameters regulating the satellite HOD will be constrained by clustering.

Literature comparison. In Fig. 3, we show COSMOS clustering measurements at 0.5 <  z <  0.8 from McCracken et al. (2015) and Leauthaud et al. (2012), and in the HSC-SSP Wide survey from Ishikawa et al. (2020). McCracken et al. (2015) used the 1.5 deg2 COSMOS footprint of UltraVISTA DR1 (McCracken et al. 2012) to measure clustering for mass-threshold selected samples. We also show the correlation function for galaxies with log M*/M >  9.4 and log M*/M >  11.0 in black wedges and triangles. Qualitatively, measurements from the literature are in agreement with our work, although the 2PCF for the log M*/M >  11.0 sample has a slightly higher amplitude, especially at small scales. Leauthaud et al. (2012) used the 1.64 deg2 of COSMOS, as imaged by HST/ACS in F814W (Koekemoer et al. 2007), to measure the 2PCF for mass-threshold samples at 0.48 <  z <  0.74. The measurement for log M*/M >  11.1 is shown as gray hexagons in Fig. 3, which are consistent with our measurements and those of McCracken et al. (2015). Finally, in green diamonds, we show the measurements from Ishikawa et al. (2020) in 145 deg2 in the HSC-SSP Wide for a sample of log M*/M >  10.1 galaxies at 0.55 <  z <  0.80. The amplitude of Ishikawa et al. (2020) 2PCF corresponds to what we measure in this work between log M*/M >  10.3 and log M*/M >  11.0. This could come from incompleteness in their samples, or/and uncertain stellar masses that were estimated with optical (grizy) bands only.

Redshift evolution of the clustering. With respect to redshift, to show a possible evolution of the clustering amplitude, in Fig. 5 we recompute the correlation function for galaxies selected above the same mass threshold in all z-bins: M* >  1010M. Although the clustering amplitude of dark matter decreases with increasing redshift, the evolution of the clustering amplitude for galaxy samples selected at the same mass-threshold depends on the galaxy formation model. The clustering of galaxies depends on how galaxies occupy DM halos, which can change with redshift. N-body simulations combined with semi-analytical models of galaxy formation indicate that the clustering amplitude of similarly selected galaxies first decreases from z = 0 to z = 1.5, remains constant up to z = 2.5, and then increases again at higher redshifts (Kauffmann et al. 1999). Qualitatively, this behavior can be observed in our measurements in Fig. 5: the correlation amplitude is the highest in the lowest redshift bins, reaches the lowest amplitude for intermediate z-bins of about z ∼ 1.5 and then increases again at z >  2.0.

thumbnail Fig. 5.

Correlation of galaxies with M* >  1010M (left panel) and GSMF (right panel) for all ten redshift bins. The green dashed lines in the right panel correspond to the GSMF of Davidzon et al. (2017).

Redshift evolution of the SMF. The SMF measurements in Fig. 5 also show the usual evolution with redshift (see e.g., Ilbert et al. 2013; Davidzon et al. 2017): the normalization decreases and the knee at M* ∼ 1011M becomes less and less prominent with increasing redshift; the slope of the low-mass end remains constant up to z = 2 but steepens at higher redshifts where the SMFs resemble more a power-law (e.g., in the z >  4.5 bin); the redshift evolution is strongly dependent on mass: the low-mass end evolves more rapidly than the high-mass end. The SMFs, having the most constraining power over the model parameters (due to the small measurement errors and sensitivity of the model), show an excellent fit of the models to the measurements. The dashed lines in Fig. 4 show the SMFs measured by Davidzon et al. (2017) using the previous version of the catalog, COSMOS2015. Overall, they are in agreement with our measurements over the whole redshift range.

5.3. Evolution of the mean halo occupation with redshift

The mean halo occupations, as defined by Eqs. (10) and (11), are shown in Fig. 6 for 0.2 <  z <  0.5. We show the mean number of galaxies in four mass bins: log M*/M = {[9.0,  9.5], [9.5,  10.0], [10.0,  10.5], [10.5,  11.0]}, as a function of halo mass for all galaxies in the thick solid lines, and for satellites and centrals in dotted and dashed lines, respectively. It is immediately evident that the mean halo occupation shifts toward high halo masses for more massive galaxies, as it requires more massive halos to host more massive galaxies. Furthermore, the central occupation peaks at some characteristic mass. Halos that have this characteristic mass can be considered as most likely to host a central galaxy in a given stellar mass bin.

thumbnail Fig. 6.

Mean number of galaxies with stellar masses in a given mass bin as a function of the mass of the halo that they occupy. We show the mean halo occupation function for galaxies in 4 stellar mass bins (color coded accordingly) at 0.2 <  z <  0.5. The thick solid lines show the total ⟨Ntot⟩ and the dashed and dotted lines show the centrals and satellites.

As the halo mass increases, the number of satellites starts to increase sharply. The mean occupation for low-mass galaxies shows that there can be halos of intermediate mass that do not host any low-mass galaxies. For example, halos of Mh ∼ 1012M have a very low probability of hosting of 109.5 <  M*/M <  1010 galaxies. The central and satellite decompositions (dashed and dotted lines) show that this is because galaxies in this mass bin cannot be centrals in Mh ∼ 1012M halos and can only be satellites in even more massive halos. We also note that as their stellar mass increases, central galaxies are more likely to occupy halos with a larger variety of masses (looking at the dashed line, for higher mass bins there is shallowing of the slope at which the central occupation decreases with halo mass). This behavior can come from a quenching of massive galaxies – as their stellar mass growth stops, the halo they inhabit continues to grow in mass. Finally, in clusters (Mh >  1013M), low-mass satellites dominate the number of galaxies in the halo. This can also be seen as a consequence of quenching: satellites stop their growth because of quenching in the halo and remain less massive, while the halo can grow by merging with other halos containing more satellites of low masses.

To investigate the redshift evolution, in Fig. 7 we show the mean occupation distribution for galaxies with 9.0 <  log M*/M <  9.5 (left panel), 9.5 <  log M*/M <  10 (middle panel) and 10.0 <  log M*/M <  10.5 (right panel) as a function of redshift at three different halo masses log Mh/M = [12.0,  12.5,  13.0]. Dashed and dotted lines show the central and satellite mean halo occupations, while the points connected with transparent solid line show the total. The panels show that the total ⟨N(Mh/M*)⟩ of Mh ≤ 1012M halos is dominated by centrals at all redshifts, whereas satellite dominate at higher halo masses at all redshifts. An exception are galaxies with 9.0 <  log M*/M <  9.5 (left panel) which are found as satellites in halos of Mh ≥ 1012M and at all redshifts. In each panel and for every halo mass, we detect little-to-no evolution of the mean occupation number, in accordance with previous findings based on N-body simulations (e.g., Kravtsov et al. 2004). At z >  2.5 there are variations toward higher mean halo occupation number but with overly large uncertainties to be significant. These results indicate that (statistically) in terms of mean occupation numbers, galaxies populate DM halos in the same way throughout cosmic time.

thumbnail Fig. 7.

Mean halo occupation in halos of a given mass as a function of redshift. In each panel, we show the ⟨Ntot(Mh)⟩ at log Mh/M = [12.0,  12.5. 13.0]. The three different panels show the mean occupations by galaxies in three different stellar mass bins log M*/M = {[9.0,  9.5], [9.5,  10.0], [10.0,  10.5]}. The dashed and dotted lines show the central and satellite occupations, while the points connected with transparent solid line show the total.

5.4. Satellite fraction and its evolution with redshift

Dark matter halos are usually inhabited by a massive central galaxy and a number of smaller satellite galaxies orbiting the potential well of the halo. At a fixed stellar mass, a galaxy can be either a central in a relatively low-mass halo or a satellite in a massive one. The number of satellite galaxies in a halo and its evolution with redshift reflects the halo’s evolutionary history in terms of its hierarchical merger tree, but it also reflects the physical processes and environmental effects that can affect the assembly of satellites. Using our constraints on the HOD in the broad redshift span up to z ∼ 5, we can study the evolution of the satellite fraction and get insights into the halos’ evolutionary history.

Within the HOD framework, we can compute the fraction of satellite galaxies, summed over all halos and with masses above a given stellar mass threshold; then, using our best-fit parameters in the ten z-bins reconstruct its evolution with redshift. To compute the satellite fraction, we perform the following integration:

(19)

where, as before, dn/dMh is the halo mass function, g is the mean number density of galaxies with and is the mean occupation function for centrals with the best-fit parameters.

Our results on satellite fraction of galaxies with masses above log M*/M >  [8.7,  9.0,  9.3,  9.7,  10.2,  10.7,  11.0] as a function of redshift are shown in Fig. 8. The general trend at all mass thresholds is an increase of the satellite fraction as cosmic time flows. For example, galaxies with masses log M*/M >  9.7 see an increase from about 10% at z ∼ 3 to ∼18% at z ∼ 1.5 all the way up to ∼25% at z ∼ 0.9. In the lowest bin 0.2 <  z <  0.5, fsat appears to systematically drop by about 3–4% for all stellar mass thresholds. This is likely results from a feature in the data, since the survey is not optimized for low redshifts. The fraction of satellites depends on the stellar mass threshold – at all redshifts there are more low-mass satellites than high-mass ones. Furthermore, the increase with redshift is different with respect to the stellar mass threshold – the fraction of high mass satellites increases more slowly, only to reach ∼8% at z ∼ 0.6. We note that the fsat in 2.5 <  z <  3.0 have all very similar values, which is an artifact arising from systematic errors in the HOD parameters. We investigated that this is mainly driven by βcut parameter which is poorly constrained.

thumbnail Fig. 8.

Fraction of satellite galaxies with masses above a given threshold as a function of redshift (solid lines). In each of the three panels, we compare with the satellite fractions measured in the hydrodynamical simulations (dashed lines) TNG100 (left), HORIZON-AGN (center), and EAGLE (right). The redshift evolution of the satellite fraction is shown for galaxies with masses above log M*/M >  [8.7,  9.0,  9.3,  9.7,  10.2,  10.7,  11.0].

The satellite fraction fsat as a function of stellar mass threshold rises sharply from very massive to intermediate-mass satellites but then reaches a plateau for intermediate to low-mass thresholds, especially at low z. This can be explained by the fact that low-mass galaxies are preferentially central galaxies in smaller halos rather than being satellites in more massive halos. This can be understood considering the halo mass function and the halo occupation function: even though the number of satellites increases as a power law with halo mass, there are simply more low-mass halos that can host a lower mass central; furthermore, the exponential high-mass cut-off of the halo mass function means that high mass halos that can host many low-mass satellites are very rare. Therefore, at a fixed low redshift the satellite fraction increases with decreasing stellar mass threshold and reaches a plateau at about 30%.

5.5. Inferred SHMR for centrals

The SHMR and M*/Mh ratio for centrals are shown in Fig. 9 at all z in the top and bottom panels, respectively. The shaded region envelops the 16th and 84th percentiles of the distribution of M* at a given Mh that is obtained by plugging in the parameters of the MCMC chain in Eq. (9). The solid line corresponds to the 50th percentile of this distribution. In the remainder of the paper, the 1σ confidence intervals are always computed in this way, unless stated otherwise. On the right-hand side of the y-axis, we show the corresponding halo star-formation efficiency (SFE) in percentages.

thumbnail Fig. 9.

Stellar-to-halo mass relation (top) and M*/Mh ratio (bottom) in the ten redshift bins. The solid lines and shaded regions show our inferred SHMR and 1σ confidence interval color-coded according to the redshift bin.

The SHMR increases monotonically with halo mass, changing slope at Mh ≈ 1012M and M* ≈ 5 × 1010M. Below this pivot mass, the SHMR increases steeply with a slop that remains constant with redshift. Above the pivot mass, the slope suddenly decreases and the stellar mass increases more slowly with halo mass. The SHMR is higher at low-z for masses below the pivot, and lower at low-z for masses above the pivot.

The M*/Mh ratio, which can be considered as the star-formation efficiency integrated over the halo’s lifetime, strongly depends on halo mass. The SFE can be defined as to quantify how efficiently baryons are converted into stars in galaxies residing in halos of a given mass – it is essentially a ratio between the star-formation rate and halo growth rate multiplied by the universal baryonic fraction. Our results, in line with previous findings, show that at all halo masses and at redshifts at least up to z ∼ 3 the SFE is lower than 20%, indicating a globally inefficient galaxy formation process. In the last three z-bins above z >  3, our results become very uncertain – the large error bars on the fitted parameters propagate into large uncertainties on the SHMR that make the interpretation difficult. This can be due to increasingly smaller sample, especially of high-mass galaxies, as well as uncertainties in the physical parameters and possible cosmic variance effects.

The SFE peaks at 17% occurs at halo masses of about Mh = 2 × 1012M. It then decreases rapidly at lower and higher halo masses – about a 15% decrease in SFE for a decrease of 1 dex in halo mass, and a ∼10% decrease for an increase of 1 dex in halo mass. This behaviour indicates that the majority (around two-thirds) of star-formation occurs in a relatively narrow range of halo masses around this peak (see e.g., Behroozi et al. 2013, 2019). This peak corresponds to stellar mass of about M* = 5 × 1010M, which is the typical M* mass scale of Milky Way-like galaxies (∼6 × 1010M, Licquia & Newman 2015). With respect to redshift, the peak SFE shows only a mild evolution, generally toward lower values with increasing redshift. This is further discussed in the next subections.

5.6. Redshift evolution of the peak mass quantities

The peak halo mass () in the M*/Mh ratio represents the mass at which the galaxy formation process, integrated over the entire history of the halo, has been most efficient. Since the feedback mechanisms also depend on halo mass, the redshift evolution of informs us about the halo mass scales at which different feedback mechanisms become more important throughout cosmic time. We compute the peak SFE from the central M*/Mh.

Redshift evolution of the peak halo mass. Figure 10 shows the redshift evolution of inferred from our analysis, compared to a compilation of measurements from the literature. To obtain and its error bars we compute the peak mass for each parameter set of the MCMC samples; then from this distribution we compute the median, 16th and 84th percentile. Our results show that the peak halo mass increases with redshift from at z = 0.35 at z = 2.75. The peak halo mass continues to increase up to in our highest bin at z = 5. At z >  3, the uncertainty of the peak position increases due to the large uncertainties in the SHMR. While at low redshifts the peak value and evolution is in agreement with the literature, at z >  3 there is a large scatter in the literature with values ranging from Mh ∼ 1012M, as found by Behroozi et al. (2013), to Mh ∼ 2.5 × 1013M, as found by Legrand et al. (2019).

thumbnail Fig. 10.

Evolution of the peak halo (left panel) and peak stellar mass (right panel) with redshift. The results from our analysis are shown in purple points. For comparison, we show measurements from the literature rescaled to match our chosen value for H0 = 70 km s−1 Mpc−1. The literature measurements include Legrand et al. (2019), Leauthaud et al. (2012, L+12), Coupon et al. (2012, C+12), Coupon et al. (2015, C+15), Cowley et al. (2018, C+18), Moster et al. (2013, M+13), Behroozi et al. (2013, B+13), Behroozi et al. (2019, B+19), and from the hydrodynamic simulations HORIZON-AGN, TNG100, and EAGLE (references in the main text).

Redshift evolution of the peak stellar mass. The right panel of Fig. 10 shows the evolution of the peak stellar mass . At the peak stellar mass, galaxies can be considered to have been the most efficient in converting baryons to stars. We find an increase of from 3.1 × 1010M at z = 0.35 to 8.7 × 1010M at z = 2.75. This increase means that at earlier times more massive galaxies have been more efficient in the star-formation process; as time elapses, this efficiency moves toward lower mass galaxies. However, the co-evolution of both the peak halo and stellar mass leaves the M*/Mh ratio nearly constant with time (further discussed in Sect. 5.8).

The trends of increase with increasing redshift of both the peak halo and peak stellar mass is a signature of the downsizing scenario. Downsizing, in its most general sense, refers to the decrease with time of some mass scale parameter that is related to stellar growth or star-formation (Cowie et al. 1996). In our case, these mass scale parameters are the peak halo and stellar mass, which are related to the efficiency of the star-formation process. Their increase with redshift means that higher mass halos and galaxies were more efficient in converting the baryon reservoir to stars at higher redshifts. Consequently, the feedback mechanisms, especially the ones active at the massive end, were less efficient in the past.

Literature comparison. From the literature compilation, we remark on comparisons with the following works. Legrand et al. (2019) used the previous iteration of the photometric catalog in the COSMOS field – COSMOS2015 – to infer the SHMR by fitting the same functional form using parametric sub-halo abundance matching. Their results (shown in green squares) are in close agreement with our results. Next, in three z-bins up to z <  1, we include the results from Leauthaud et al. (2012), which serve as our main reference for the theoretical modeling. Their analysis is based on a joint abundance, clustering, and galaxy-galaxy lensing fit on measurements done in the COSMOS field; their results are shown in dark blue circles. Unsurprisingly, our results are in agreement with the trend. The higher value in the z ∼ 0.65 bin found by our work, Leauthaud et al. (2012) and Legrand et al. (2019), is likely a feature of the COSMOS field. Indeed, using 10 000 spectroscopic redshifts, Kovač et al. (2010) reported a very large overdense structure in COSMOS at these redshifts. Behroozi et al. (2013, 2019) used empirical modeling where galaxies are populated in dark matter halos and are traced within their halos over time; the models are constrained to match a set of observables such as the SMF, luminosity function, and cosmic star formation rate, among others. Their results are in general agreement at low redshifts, but show that the peak halo mass turns over and starts to decrease at z ∼ 2 or z ∼ 3. Our results, as well as Legrand et al. (2019), on the other hand, show a peak halo mass increasing up to z ∼ 4. It is possible that this effect is driven by effectively the same COSMOS dataset used in ours and Legrand et al. (2019) analyses. Behroozi et al. (2019), for example, constrained the SHMR at high redshifts using GSMF derived from Song et al. (2016) data. It is important to note that they used UV-to-M* conversions to estimate stellar masses, which comes with some caveats. The SHMRs at z ∼ 3−4 of Behroozi et al. (2019) method are sensitive to the GSMF at z >  4, and clearly sensitive to the choice of observational constraints (see e.g., Behroozi & Silk 2015).

Finally, in the left panel of Fig. 10, we compare the peak stellar mass with the literature. Our results are in good agreement with Leauthaud et al. (2012). Interestingly, no hydrodynamical simulations show a clear downsizing trend in the peak stellar mass.

In summary, our results show that the peak halo mass increases monotonically with redshift in agreement with the literature up to z <  3, including findings in hydrodynamic simulations. At higher redshift, the literature suggests a turnover and decrease of the peak halo mass, which is not captured by our analysis. The increase of both peak halo and peak stellar mass with redshifts is in accordance with the downsizing scenario, where the more massive halos and galaxies were more efficient in star-forming earlier in the universe and the peak efficiency shifts toward lower mass halos.

5.7. Total SHMR

As the mass of the halo increases, the number of satellite galaxies that occupy it also increases and, naturally, their contribution to the total stellar mass budget of the halo becomes important. In massive halos, stellar mass is assembled from in-situ star-formation and from mass accretion via merging of halos, while the growth is regulated by various quenching mechanisms. The ratio between the total stellar mass and the halo mass (total SHMR) can then inform us about the efficiency of the combination of both effects. The model adopted in this study allows us to compute the total stellar mass contained in a halo of a given mass using Eq. (14). To obtain the total stellar mass, an integration is carried out over the stellar masses with lower and upper mass limits. Ideally, we would integrate over the whole range of possible stellar masses, but in our case that would mean an extrapolation of the models beyond the stellar masses probed in our analysis. This can introduce inaccuracies, especially for computing the satellite contribution. However, we checked that most of the contribution to the total stellar mass content in Mh >  1012M halos comes from satellites in the mass range of 1010M <  M* <  1011M, well within the mass scales probed by our analysis; this is also stated in Leauthaud et al. (2012). For our purposes, to compute , we set the lower integration limit to M* = 108.5M. This lower stellar mass limit is below our completeness limit at z >  2, but we expect that the extrapolation at lower masses is not inaccurate enough to to bias the results.

Figure 11 shows our results on the total stellar content as a function of the host halo mass in the top panel and ratio in the bottom panel in 0.2 <  z <  0.5 and for the other nine z-bins in Fig. 12. We show the central SHMR in dashed purple, while the dotted purple line shows the satellite contribution. The shaded purple area envelops the 1σ uncertainty in the total SHMR. In Fig. 12, both the central and satellite are displayed in solid purple with the 1σ envelope of the total ratio. The break in the lines and in the shaded region indicates the upper stellar mass limit probed by our analysis, as well as the extrapolation to higher masses is shown in transparent purple lines and hatched region. On the right-hand side y-axis, we show the integrated SFE. We recall here that by definition, the integrated SFE describes how efficiently stellar mass has been assembled in halos integrated over the halo’s lifetime. This inevitably mixes various assembly paths, such as stellar mass formed in low-mass halos (where different mechanisms regulate growth) that are later accreted as satellites in massive halos.

thumbnail Fig. 11.

Total stellar-to-halo mass relation (top panel) and total M*/Mh (bottom panel) at 0.2 <  z <  0.5 compared to hydrodynamical simulations. The purple dashed and dotted lines show our central and satellite contribution to the total SHMR respectively, with the shaded region showing 1σ confidence interval of the sum of the two. The break in solid purple lines and shaded regions indicate the highest stellar mass probed in our analysis, which we take to be the highest mass bin in the SMF. The transparent purple lines and hatched region is an extrapolation at higher masses. We overplot the SHMRs measured in the hydrodynamical simulations HORIZON-AGN in teal, TNG100 in red, and EAGLE in dark yellow, where the dashed dotted and solid lines show the central, satellite, and total SHMR.

thumbnail Fig. 12.

Total M*/Mh in the different redshift bins compared to hydrodynamical simulations. The purple lines show our inferred central and satellite contribution to the total SHMR with the shaded region showing a 1σ confidence interval of the sum of the two. The break in solid purple lines and shaded regions indicate the highest stellar mass probed in our analysis, which we take to be the highest mass bin in the SMF. The dashed purple lines and hatched region is an extrapolation at higher masses. The comparison includes total M*/Mh found in the hydrodynamical simulations HORIZON-AGN in teal, TNG100 in red and EAGLE in dark yellow. The dashed and dotted lines for the simulations indicate the central and satellite contributions, respectively.

At masses below Mh ≲ 1013M, the total stellar mass content is completely dominated by central galaxies. It rises sharply up to the peak halo mass at around Mh ∼ 1012M, and then falls more gradually with increasing halo mass. The peak of the total ratio is completely set by the centrals, meaning that the physical processes shaping the peak have to be related to the quenching of the central galaxy. Since the satellite contribution only becomes important at halo masses almost one order of magnitude higher, the accumulation of stellar mass in satellites (instead of the central growth or mergers), cannot be responsible for setting the peak efficiency.

For masses higher than Mh ∼ 1012M, the satellite contribution to rise as the central contribution drops and a transition occurs at about Mh = 2 × 1013M where satellites start to dominate. Going to higher masses, the satellite starts to flatten out, meaning that the stellar mass keeps up with the growth rate of halos in group- and cluster-scale halos. Excluding the 2.5 <  z <  3.5 and 4.5 <  z <  5.5 bins, the total ratio always remains below the peak set by the centrals at an SFE below 5%. This suggests that even if all the satellites were to merge into the central, the SFE would still be lower than the peak, indicating strong environmental quenching mechanisms in massive halos. In other words, the accumulation of mass in satellites is not responsible for the low M*/Mh.

5.8. Redshift dependence of M*/Mh at fixed halo mass

The M*/Mh ratio (i.e., the integrated star-formation efficiency) might evolve with redshift differently depending on the mass of the halo. This dependence, especially when considering the contributions from both centrals and satellites, could shed light on the importance of the feedback mechanisms that regulate star-formation at different halo mass-scales. Figure 13 shows the M*/Mh ratio as a function of redshift at different halo masses log Mh/M = [11.50,  12.00,  13.00,  13.60], as well as for . For the most massive halos at Mh = 1013.6M, we also decomposed it in the contributions from satellites and centrals. We restricted our analysis to z <  2.5 since the large uncertainties at higher z prohibit a meaningful quantitative analysis.

thumbnail Fig. 13.

Redshift dependence of M*/Mh at fixed halo mass. We show the total (centrals and satellites) M*/Mh fixed at log Mh/M = [11.50,  12.00,  13.00,  13.60, ], while the M*/Mh for log Mh/M = 13.60 is also shown for centrals and satellites separately with star and diamond symbols, respectively.

Low-mass halos, below the peak halo mass (Mh <  1012M), show their SFE steadily increasing from z ∼ 2.3 up to present day. For Mh = 1011.5M the SFE goes from ∼2.5% at z = 2.3 to ∼7% at z = 0.3. For slightly more massive halos of Mh = 1012M, the SFE increases even faster from ∼8% at z = 2.3 to ∼17% at z = 0.3. In contrast, for high-mass halos above the peak halo mass (Mh >  1013M), the SFE remains almost constant, with a slight decrease with decreasing redshift. Furthermore, more massive halos show an even lower SFE by several percents.

These trends are a signature of the downsizing scenario, which was already mentioned in Sect. 5.6. Downsizing refers to the observation that, contrary to the hierarchical formation scenario in which small halos are formed first and subsequently grow by merging and accretion, more massive and early-type galaxies have stellar populations that are formed earlier (De Lucia et al. 2006; Thomas et al. 2010). This downsizing is observed in our results from the increase of the SFE of low-mass halos and the slight decrease in the efficiency of high-mass halos with decreasing redshift. Low-mass halos having been more efficient in forming stars at later times means that the stellar populations of lower mass galaxies inhabiting them would also be younger.

The SFE at the peak halo mass shows only a weak evolution with redshift, remaining constant to z = 0.95 then dropping from 20% to 16% at z = 2.5. Previous findings also point to a peak efficiency constant with redshift (Behroozi et al. 2019; Moster et al. 2018). This behavior suggests that the can be considered as the halo mass scale at which the integrated SFE history remains constant with redshift.

Satellites dominate the mass budget in high-mass halos, as can be seen in Fig. 13 for Mh = 1013.6M where we also show M*/Mh for centrals and satellites separately. The satellite dominance is stronger at lower redshift, whereas from z >  1.2, the central contribution catches up and both remain comparable. This is unsurprising given the fact that the satellite fraction increases with decreasing redshift.

6. Discussion

6.1. Physical mechanisms that regulate the stellar mass assembly of centrals

The shape of the SHMR as a function of halo mass can provide us with qualitative information on the stellar mass growth mechanisms. For example, at Mh <  1012M, the steep increase in the SHMR with halo mass tells us that the stellar mass growth in galaxies is mainly driven by in situ star formation, as opposed to mergers (Conroy & Wechsler 2009; Leauthaud et al. 2012). The explanation comes from the fact that halos grow by merging with lower mass halos, in which the stellar mass drops as a power law (very low mass halos can even be devoid of stars). Since merging will increase the halo mass more significantly compared to the stellar mass, the stellar mass growth in low-mass halos has to be driven by star formation in order to obtain the steep rise of the SHMR.

The shape of the SFE with its characteristic peak is mostly driven by different feedback processes, whose efficiencies are dependent on halo mass. At low masses, stellar feedback from supernovae (SNe), stellar winds, radiation pressure, and photoheating can heat up and prevent the baryons from collapsing into stars (e.g., Dekel & Silk 1986; Dubois & Teyssier 2008; Hopkins et al. 2014; Kimm et al. 2015). As the halo mass increases, its potential well deepens and cold flows feed the halo with gas that fuels star-formation. The SNe are not powerful enough to prevent these cold flows, so the SFE increases with halo mass. As the halo mass increases even further, virial shock heating and AGN activity heat up the gas and quench star formation, thus driving the SFE back down (see e.g., Silk et al. 2014, for a review). Vogelsberger et al. (2013) show in their hydrodynamical simulations that the mechanical, radio-mode AGN feedback is the most responsible for the suppression of star formation in these massive halos (see also e.g., Dubois et al. 2010). Similar conclusions have been reached from AGN observations in the COSMOS field, where Vardoulaki et al. (2021) show that the radio-mode AGN feedback plays a significant role in shaping the SHMR at lower redshifts. Additionally, Gabor & Davé (2015) in their hydrodynamical simulations, these authors show that galaxies above the peak mass indeed tend to be quenched. These quenched centrals can still continue to grow in stellar mass via dry merging, but halo mass growth being faster – since halos keep accreting from their large-scale environment, while the gas accreted by galaxies is not efficiently turned into stars any longer due to feedback – leads the M*/Mh ratio to decrease.

In Fig. 10, we showed that more massive halos are more efficient in forming stars at higher redshifts – the peak halo mass increases with z. This peak in the central SHMR is shaped by the interplay between cold flows being more efficient as the halo mass increases and shock heating in the hot halo. The increase of the peak halo mass is a signature that cold flows become more important in driving star formation at z >  1.5 in Mh >  1012M (Dekel et al. 2009; Oser et al. 2010; Dubois et al. 2013).

With respect to redshift, we observe in Fig. 12 that the peak generally flattens out at higher redshift. This is directly connected to the fact that the SMF at higher redshifts increases its low mass slope and the knee is smoothed (i.e., the SMF and HMF become more similar in shape). This indicates that the high mass halos decrease in their star-forming efficiency with time and that the AGN feedback has a larger impact at later times.

6.2. Physical mechanisms that regulate the stellar mass assembly in satellites

In Sect. 5.8 and Fig. 13, we presented how the total SHMR ratio evolves with redshift at fixed halo mass. Comparing the contributions of centrals and satellites to the total stellar mass content at different Mh can shed some light on the relative importance of the quenching mechanisms acting on the central and satellite galaxies in a halo. Figure 13 shows that in high-mass halos of about Mh = 1013.6M, the M*/Mh is about 0.2 dex lower than M*/Mh at the peak halo mass. At least up to z <  2.5, the total M*/Mh ratio for high-mass halos is always lower than M*/Mh at the peak halo mass. This suggests that the low M*/Mh ratio for centrals is not due to the stellar mass being accumulated in satellites instead of the central because, even if all satellites were to merge with the central, the total stellar mass content would still be lower than at the peak. This is a clear indication of powerful feedback present in massive halos that prevents gas from cooling and forming stars, and these feedback processes are acting on both the central and the satellites.

In halos more massive than ∼1012M, which includes group- and cluster-scale halos, the virial shock heats the halo gas (Birnboim & Dekel 2003). Naturally, then, satellite galaxies that reside in these massive halos have their SF impeded by the hot gas. This hot-halo mode is likely to be responsible for the satellite quenching. This mechanism is also called environment quenching (Peng et al. 2010; Gabor & Davé 2015) and is found to be independent of stellar and halo mass – galaxies of all stellar masses reside in the densest environments in halos of Mh >  1012. Additional environmental effects such as strangulation (Larson et al. 1980; Balogh et al. 2000), ram-pressure stripping (Gunn & Gott 1972), and harassment (Moore et al. 1996) further prevent the stellar mass assembly in satellites. Altogether, these processes could explain the low SFE at high halo masses in our results. Our results at high masses remain too uncertain to unveil a consistent picture of the redshift evolution of the efficiency of the satellite quenching mechanisms. Incorporating lensing measurements to better constrain the satellite contribution and separating active versus passive populations could shed more light on this aspect. This is a natural extension to this work that we plan to carry out in the future.

6.3. Comparison with hydrodynamical simulations

We compare our results with several state-of-art cosmological hydrodynamical simulations of galaxy formation: HORIZON-AGN7 (Dubois et al. 2014; Kaviraj et al. 2017), EAGLE8 (Crain et al. 2015; Schaye et al. 2015), and TNG100 of the ILLUSTRISTNG project9(Nelson et al. 2018, 2019; Springel et al. 2018; Pillepich et al. 2018b; Marinacci et al. 2018; Naiman et al. 2018). A brief recap of the main features of these simulations is given in Appendix C. Once central and satellite galaxies are matched with DM halos as described in Appendix C, the SHMR is directly measured in halo mass bins at different redshifts, as the mean of the ratio between the stellar and the halo mass. Baryonic physics has a small but real impact on the halo mass function. Beltz-Mohrmann & Berlind (2021) showed that the halo mass function in DM-only simulations is overestimated (with respect to hydrodynamical simulations) by a factor of about ∼1.1–1.2 depending on the mass range at z = 0, though this effect disappears at high redhsift and high mass (see also e.g., Bocquet et al. 2016; Desmond et al. 2017; Castro et al. 2020). This aspect should be kept in mind, knowing that the modeling presented in Sect. 4 relies on DM-only halo mass function, while for the hydrodynamical simulations, the mass of the DM host halos is naturally impacted by baryonic physics. This could lead to a small systematic underestimation of the SHMR in the observational data with respect to the simulated ones.

Redshift evolution of the peak halo mass. In Fig. 10, we show the evolution of the peak halo mass in the three hydrodynamic simulations HORIZON-AGN, TNG100, and EAGLE. We directly compute the peak halo mass using the measured SHMR in the simulations. The errorbars are obtained by computing the peak of the lower and upper 1σ values of the M*/Mh relation as measured in the simulations. Here, TNG100 shows the best agreement with both our results and the literature, whereas EAGLE shows slightly higher peak halo masses up to z ∼ 3. We also note that in both simulations, the peak halo mass decreases from z ∼ 3 to z ∼ 4, which is a priori in disagreement with our results. That being said, the ratio M*/Mh at z >  3 is also highly uncertain (as seen in Fig. 12) because the simulations start to miss massive halos with masses comparable to the peak due to their fixed small volume and the rarity of these halos. HORIZON-AGN shows the biggest discrepancy with our results and the literature. In this simulation, the peak efficiency is reached in much lower mass halos: at z ∼ 0.35 that decreases up 1011.55M to z ∼ 2 and then increase at higher redshifts. This discrepancy can be explained by the inefficiency of the SN feedback in low mass galaxies as implemented in the simulation (see e.g., Hatfield et al. 2019 or Kaviraj et al. 2017 for further discussions).

The simulated SHMR in the low-mass regime. We present the for HORIZON-AGN (in teal), TNG100 (red), and EAGLE (dark yellow color) in Figs. 11 and 12. At the low-mass end, where increases with halo mass, TNG100 and EAGLE show a reasonable agreement with our results with a similar normalization, slope and position of the peak, albeit the peak SFE is lower by a few percent. HORIZON-AGN, on the other hand, has a higher normalization of the and a peak efficiency at lower halo masses (as already noted above) probably driven by insufficient SN feedback especially at high redshift. This issue is also responsible for the overall overestimation of the stellar mass in low-mass galaxies. The good agreement of TNG100 and EAGLE with our results in this regime suggests that the implementation of feedback from SNe and stellar winds in these simulations is well calibrated to reproduce observational data. We note that we cannot conclude from this measurement that the star-formation model at the subgrid scale is realistic nor that stellar feedback in these simulations is physically meaningful (some feedback mechanisms could be missing, such as radiation from young stars that suppresses star formation in low-mass galaxies, and others could be overestimated). What we can conclude is simply that the cumulative strength of the implemented feedback processes leads to realistic SHMR relations (see the discussion on subgrid models in Appendix C.2).

Simulated SHMR in the high-mass regime. The most significant discrepancies between the observational data and all simulations appear for masses above the peak. At least up to z = 3.5 the central for all the simulations shows a peak, meaning that AGN feedback are powerful enough to quench the central galaxy growth compared to the halo growth. However, at all redshifts in TNG100 and EAGLE, the central M*/Mh decreases with increasing halo mass more gradually (shallower slope) than in observations, while HORIZON-AGN shows this tendency only at z <  1. Furthermore, the contribution of satellites relatively to central at all redshifts is higher in the simulations. These two facts contribute to a flattening of the peak and a higher SFE with increasing mass in the simulations, especially in TNG100 and EAGLE.

Redshift evolution of the satellite fraction. In Fig. 8 we also compare the satellite fractions against the three hydrodynamical simulations. The dashed lines in the left, middle, and right panels show the fsat in TNG100, HORIZON-AGN, and EAGLE, where the different color-coded lines correspond to the mass thresholds indicated in the legend. The satellite fraction in all simulations was computed with the same criterion as in Eq. (19) – the ratio between the number of M* >  Mth galaxies that are satellites and the total number of M* >  Mth galaxies for all the halos in the simulation. The satellite fraction in all the simulations increases linearly with decreasing redshift, a trend which is in agreement with our results. However, at fixed stellar mass thresholds, the simulations exhibit higher numbers of satellites at all redshifts compared to our analysis. In TNG100 and EAGLE, fsat goes from ∼20% at z ∼ 3 to ∼42% at z ∼ 0.3 for satellites with log M*/M >  9.0. Compared to our results, the satellite fraction is higher by about a factor of 2 for the lowest mass thresholds and a factor of 1.2 for the highest mass thresholds. In HORIZON-AGN, the satellite fraction is lower by about 0.5 compared to the other two simulations, but still a factor of 1.4 higher than our results for the lowest mass thresholds.

The difference in the satellite fractions can be due to a degeneracy between the satellite fraction and the spatial distribution of satellites within halos (Xu et al. 2016) – steeper satellite profiles at small scales increase the correlation, thus decreasing the need for high satellite fraction. This raises the possibility that the higher satellite fraction in the simulations is due to flatter satellite profiles compared to our models which assume an NFW. Bose et al. (2019) show that the mean radial number density of luminous satellites in TNG100 matches the NFW profile. Furthermore, we compared the radial satellite density distribution in all three simulation and find that they follow the NFW profiles reasonably well, with TNG100 showing the best match, while the EAGLE and HORIZON-AGN show flatter profiles that are consistent with NFW but for lower concentration parameters. The agreement in TNG100 suggests that satellite profiles cannot serve as the explanation behind the different satellite fractions. In EAGLE and HORIZON-AGN, on the other hand, the flattening of the profiles can be interpreted as a consequence of AGN feedback (see e.g. Peirani et al. 2017, for comparison of DM halo profiles with and without AGN).

Excess of satellites in halos “above fixed mass thresholds” would translate into a higher level of small-scale clustering, which is measurable by the two-point correlation function of low-mass and red galaxies. Artale et al. (2017) compares the small-scale correlation function of galaxies in EAGLE and the GAMA survey at z = 0.1. They confirm that low-mass red galaxies have a considerably higher correlation at small scales, which is consistent with high fractions of satellite galaxies. Although they did not investigate this feature specifically, Springel et al. (2018) mentioned that the present-day clustering of red galaxies less massive than 1010M was overestimated in TNG100 with respect to the SDSS. This difference in the satellite fractions corroborates our measurement of the M*/Mh in Fig. 12, where the satellite contribution in the simulations is usually higher than our measurements. A higher fraction of satellites above a mass threshold could be the consequence of satellite galaxies having grown more massive due to lack of quenching, at least during a certain period of cosmic time. Due to the power-law increase of the SMF at the low-mass end, inefficient satellite quenching would shift the satellite SMF toward higher masses. Then, above a stellar mass threshold, there would be more satellites. This indicates that satellite quenching in the hydrodynamical simulations is inefficient, or that it happens at later times in the halo lifetime, or both.

It is worth noting that Donnari et al. (2021) made a thorough comparison of the quenching fraction between TNG100 simulation and observations up to z = 0.65. They concluded that TNG quenched fractions of centrals and satellites are qualitatively in agreement with observations. The scenario behind our findings might therefore lie at higher redshift: if strong satellite quenching occurs too late in the simulations, these galaxies would have had time to grow more massive than in the observations, leading to the overestimation of the satellite SHMR at high halo mass. We note that, at least for TNG100 and EAGLE, this quenching inefficiency must be specific to galaxies that have already been accreted as satellites: indeed the SHMR of low-mass central galaxies at high redshift is in good agreement with our observational measurements10. It is possible to argue that this issue is related to the resolution limit of the simulations. The coarse resolution reached in the circum-galactic medium might be insufficient to correctly model high gas temperature inherited from virial shock heating, preventing an overly strong quenching of satellite galaxies in massive halos (Gabor & Davé 2015). However, this interpretation does not hold to explain the too large contribution of satellites at high redshift in halos less massive than Mh >  1012M (which are less prone to shock heating). For those very small satellites, we could suspect that other environmental mechanisms (e.g., ram-pressure and tidal stripping, harassment) are not correctly modeled, also due to the lack of resolution in the simulations. Interestingly, Costa et al. (2019) also shows that the efficiency of satellite tidal stripping depends on the degree of pre-processing of low-mass galaxies by radiative stellar feedback (which are not implemented in the simulations studied here; see also Katz et al. 2020 for a discussion of stellar radiation during reionization). Therefore, the lack of satellite quenching at high redshift in the cosmological simulations studied here could also be a consequence of the absence of pre-processing by radiative stellar feedback.

6.4. Sources of uncertainties and the effects of model assumptions on the inferred SHMR

The measurements presented here can suffer from a number of systematic errors. Firstly, due to the relatively small volume probed by COSMOS, the effects of sample variance can lead to biased measurements. This might be the case at z ∼ 1, where several works indicate an overabundance of rich structures (McCracken et al. 2015, and references therein). Such an overdensity increases the normalization of the SMF and adds extra power on intermediate and large scales in w(θ). The effect of these overabundances on the inferred SHMR would be to decrease this ratio – indicating an even lower efficiency of converting baryons to stars – and a shift of toward higher masses.

A second source of systematic error are uncertainties in the estimation of physical parameters, for instance, stellar masses and redshifts11. Stellar mass uncertainties propagate into the Eddington bias at the high-mass end of the SMF. The effect of the Eddington bias on the inferred parameters is an increase in the value of σLog M* which sets the scatter in stellar mass at fixed halo mass.

In Fig. 14, we show the uncertainties in the SMF as a function of stellar mass and redshift. We show the fractional error (σ/Φ) from the cosmic variance (CV) and SED fitting. We do not show the Poisson errors because they are subdominant. Cosmic variance dominates the error budget by about 0.8 dex at M* <  1011M compared to the SED fitting. Both increase with mass, but the SED fitting uncertainties increase exponentially at M* >  1011M and become dominant over CV at the most massive end (M* >  6 × 1011M). The most massive galaxies are rare and reside in the densest regions, therefore, a small survey is more likely to get a biased view of these objects. However, due to photometric errors and degeneracies between different SEDs, their stellar masses come with uncertainties. Due to the exponential shape of the SMF at the high mass end, these uncertainties in the M* translate into large uncertainties in the SMF (Eddington bias) and become dominant. The SED fitting uncertainties in the M* uncertainties being amplified by the Eddington bias into large σSED means that it will be very difficult, even for future surveys, to improve on these uncertainties.

thumbnail Fig. 14.

Uncertainties in the SMF expressed as a fractional error σ/Φ. We show the uncertainties from cosmic variance with solid line-connected points and the uncertainties from SED fitting in dashed line-connected diamonds. We show the uncertainties for four redshift bins.

Additional systematic errors arise from the number of assumptions in the model. The ingredients that go into the HOD-based model of the observables are the halo abundance (HMF), their clustering properties (halo bias, bh), the radial distribution of dark matter and galaxies within halos, and halo mass definitions. The literature abounds with prescriptions for all of these ingredients, which can all lead to different results for the SHMR. This often makes the comparison with the literature difficult. Coupon et al. (2015) has investigated the effect on various model assumptions such as the σ8 value, bh(Mh) relation, assembly bias, mass concentration relation, and halo profiles, on the error budget of several HOD parameters. Their conclusion is that the model systematics can lead to errors comparable to the statistical errors.

7. Summary and conclusions

This paper presents the redshift evolution of the stellar-to-halo mass ratio to z ∼ 5, derived from measurements of the galaxy stellar mass function and angular correlation function in the COSMOS2020 catalog fitted to a phenomenological model. The advantages of our work is the use of a single, homogenous dataset to perform all the measurements. Additionally, the HOD-based modeling allows us to consistently probe the contribution to the total stellar mass budget in halos of both central and satellite galaxies over this large redshift range, for the first time in the literature.

Our principal results are as follows:

  • The mean halo occupation shows little-to-no evolution with redshift at fixed stellar mass, suggesting that galaxies occupy dark matter halos similarly throughout cosmic history.

  • The M*/Mh ratio for central galaxies – which may be interpreted as the integrated star formation efficiency (SFE) – strongly depends on halo mass, increasing up to a peak at halo masses of around 2 × 1012 and then decreasing again as the halo mass increases. The SFE shows little-to-no evolution with redshift and remains lower than 20% at least up to z ∼ 3, indicating a globally inefficient galaxy formation process. The peak levels off with increasing redshift, consistent with a scenario in which AGN feedback in higher mass halos is less important at earlier times.

  • The halo mass and stellar mass scale at which the SFE peaks, and , increase continuously with redshift at least to z ∼ 4. This stands in contrast to other works (e.g., Behroozi et al. 2019) where the peak halo mass decreases beyond z ≳ 3. However, given our errors the peak halo mass evolution at these high redshifts remains uncertain.

  • The total stellar mass content of halos, , shows that at Mh ≲ 1013M central galaxies completely dominate the stellar mass budget of the halo at all redshifts. The peak of the ratio is set by the central galaxy, indicating that the physical processes that shape the peak efficiency are related to the quenching of central galaxies.

  • Satellite galaxies start to dominate the total stellar mass budget at Mh >  2 × 1013M and the ratio flattens out as the halo mass increases. The fact that in the satellite-dominated regime, the ratio is lower than the peak means that strong quenching mechanisms must be present in massive halos that quench the mass assembly of satellites.

  • For all samples, the satellite fraction decreases at higher redshifts. Moreover, there are always more low-mass satellites than high-mass ones. The satellite fraction increases more steeply for lower-mass satellites and reaches ∼25% at low redshifts, whereas the most massive galaxies reach up to ∼15%.

We compared our SHMR measurements with three state-of-the-art hydrodynamical simulations HORIZON-AGN, TNG100, and EAGLE. For low halo masses, our results are in general agreement with TNG100 and EAGLE in terms of slope and peak of the ratio. However, there is a significant discrepancy with HORIZON-AGN that can be explained by insufficient stellar feedback. The most significant discrepancies with the simulations are for masses above the peak. The ratio in the simulations flattens out at the peak and has a larger value at higher masses compared to our results, which is mainly driven by a higher satellite contribution at all redshifts in the simulations – this excess in the satellite contribution relative to the central at the high-end is higher in TNG100 and EAGLE than in HORIZON-AGN. Furthermore, the simulations show higher fractions of satellite galaxies at all redshifts and all masses by about a factor of two. Both these findings at the high-mass end – excess of and high satellite fractions – suggest that the environmental quenching of satellites is less efficient in the simulations or that quenching occurs much later. This can be due to either an inefficient hot halo quenching mode, or from other environmental effects, such as ram-pressure or tidal stripping or harassment from other satellites, which are not well captured in the simulations – probably due to resolution effects. Lack of pre-processing by stellar radiative feedback could also have an impact.

To date, the COSMOS photometric redshift catalogs have provided the only homogeneous dataset to consistently study the evolution of the SHMR over a large redshift range (since z ∼ 5) for a statistically representative area of the sky. However, the 2 deg2 of the survey does not eliminate the effects of cosmic variance; at z >  2.5, the uncertainty in our results makes it difficult to provide a convincing interpretation and comparison with simulations.

Our analysis is based on a phenomenological model that cannot provide any further insight into the physical processes governing the shape of the SHMR, especially the relative contribution of different feedback modes. But comparisons with hydrodynamical simulations, where the effects of feedback on the SHMR can be traced back, could provide insights into the relative importance of the different feedback mechanisms acting at different mass scales and environments and, ultimately, can help bridge the gap between observations and simulations.

One of the main issues raised by our comparisons with hydrodynamical simulations concerns satellite quenching. To investigate how satellites are quenched with respect to their DM halo throughout cosmic time, we could further separate the sample into star-forming and quenched galaxies. Upcoming works and surveys will make this possible. In the imminent future, COSMOS-Web will provide JWST observations in four NIR bands down to AB ∼ 27 over 0.6 deg2 in COSMOS. This unprecedented resolution and depth in NIR will enable lensing measurements (such as galaxy-galaxy lensing) to z ∼ 2.5 and make it possible to measure the SHMR dependence on star-formation activity and even color gradients. At a slightly longer timescale, the Cosmic Dawn Survey will carry out deep NUV to MIR observations in ∼50 deg2 in the Euclid Deep Fields. Accurate photometric redshift and stellar mass measurements from this survey will probe the most massive end of the SHMR while greatly reducing cosmic variance.


1

σNMAD = 1.48 × median[(|Δz − median(Δz)|)/(1 + zspec)]; Δz = zphot − zspec.

2

Defined as |Δz|> 0.15(1 + zspec).

3

Defined as median(Δz).

4

In the catalog this is selected using the keyword FLAG_COMBINED.

5

The correlation functions are computed using the TreeCorr code (Jarvis et al. 2004).

10

The issue is different for HORIZON-AGN, as previously noted. In this simulation, all low-mass galaxies (central and satellites) are concerned by inefficient quenching, which suggests a different type of failure or mis-modeling in the simulation.

11

For a detailed analysis of the uncertainties affecting the SHMR see Behroozi et al. (2010).

12

A light cone has also been extracted on-the-fly in order to build realistic mocks (Laigle et al. 2019; Gouin et al. 2019). In particular, based on these mocks, Hatfield et al. (2019) showed that the propagation of statistical and systematic uncertainties inherited from redshift and mass photometric estimates lead to an underestimation of the clustering amplitude by ∼0.1 dex. In the present work, we use galaxy and halo catalogs extracted from the snapshot data, but we checked that using the light cone data does not significantly change the measured SHMR.

15

While this implementation of stochastic feedback enhance the impact of feedback in those (GADGET) SPH simulation, it makes no difference with standard thermal release when simulated with (RAMSES) AMR (Rosdahl et al. 2017).

Acknowledgments

This research is partly supported by the Centre National d’Études Spatiales (CNES). MS acknowledges Elena Sarpa, Lukas Furtak and Louis Legrand for useful discussions. Part of this study benefitted from earlier, unpublished work by J. Coupon. MS acknowledges the partial thesis funding by Euclid-CNES. HJMcC and CL acknowledge support from the Programme National Cosmology et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES. ID has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska- Curie grant agreement No. 896225. This work used the CANDIDE computer system at the IAP supported by grants from the PNCG, DIM-ACAV and CNES and maintained by S. Rouberol. This work is based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under ESO program ID 179.A-2005 and on data products produced by CALET and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium. This work is based in part on observations made with the NASA/ESA Hubble Space Telescope, obtained from the Data Archive at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.

References

  1. Anglés-Alcázar, D., Faucher-Giguère, C.-A., Quataert, E., et al. 2017, MNRAS, 472, L109 [Google Scholar]
  2. Applebaum, E., Brooks, A. M., Quinn, T. R., & Christensen, C. R. 2020, MNRAS, 492, 8 [CrossRef] [Google Scholar]
  3. Arnouts, S., Moscardini, L., Vanzella, E., et al. 2002, MNRAS, 329, 355 [Google Scholar]
  4. Artale, M. C., Pedrosa, S. E., Trayford, J. W., et al. 2017, MNRAS, 470, 1771 [NASA ADS] [CrossRef] [Google Scholar]
  5. Aubert, D., Pichon, C., & Colombi, S. 2004, MNRAS, 352, 376 [Google Scholar]
  6. Balogh, M. L., Navarro, J. F., & Morris, S. L. 2000, ApJ, 540, 113 [Google Scholar]
  7. Behroozi, P. S., & Silk, J. 2015, ApJ, 799, 32 [NASA ADS] [CrossRef] [Google Scholar]
  8. Behroozi, P., & Silk, J. 2018, MNRAS, 477, 5382 [Google Scholar]
  9. Behroozi, P. S., Conroy, C., & Wechsler, R. H. 2010, ApJ, 717, 379 [Google Scholar]
  10. Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57 [Google Scholar]
  11. Behroozi, P., Wechsler, R., Hearin, A., & Conroy, C. 2019, MNRAS, 488, 3143 [NASA ADS] [CrossRef] [Google Scholar]
  12. Beltz-Mohrmann, G. D., & Berlind, A. A. 2021, ApJ, 921, 112 [NASA ADS] [CrossRef] [Google Scholar]
  13. Berlind, A. A., & Weinberg, D. H. 2002, ApJ, 575, 587 [Google Scholar]
  14. Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393 [Google Scholar]
  15. Birnboim, Y., & Dekel, A. 2003, MNRAS, 345, 349 [Google Scholar]
  16. Bocquet, S., Saro, A., Dolag, K., & Mohr, J. J. 2016, MNRAS, 456, 2361 [Google Scholar]
  17. Booth, C. M., Agertz, O., Kravtsov, A. V., & Gnedin, N. Y. 2013, ApJ, 777, L16 [CrossRef] [Google Scholar]
  18. Bose, S., Eisenstein, D. J., Hernquist, L., et al. 2019, MNRAS, 490, 5693 [Google Scholar]
  19. Bower, R. G., Benson, A. J., Malbon, R., et al. 2006, MNRAS, 370, 645 [Google Scholar]
  20. Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503 [Google Scholar]
  21. Brodwin, M., Lilly, S. J., Porciani, C., et al. 2006, ApJS, 162, 20 [Google Scholar]
  22. Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000 [Google Scholar]
  23. Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80 [Google Scholar]
  24. Castro, T., Borgani, S., Dolag, K., et al. 2020, MNRAS, 500, 2316 [Google Scholar]
  25. Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]
  26. Chisari, N. E., Alonso, D., Krause, E., et al. 2019, ApJS, 242, 2 [Google Scholar]
  27. Conroy, C., & Wechsler, R. H. 2009, ApJ, 696, 620 [NASA ADS] [CrossRef] [Google Scholar]
  28. Cooray, A., & Sheth, R. 2002, Phys. Rep., 372, 1 [Google Scholar]
  29. Costa, T., Rosdahl, J., & Kimm, T. 2019, MNRAS, 489, 5181 [CrossRef] [Google Scholar]
  30. Coupon, J., Kilbinger, M., McCracken, H. J., et al. 2012, A&A, 542, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Coupon, J., Arnouts, S., van Waerbeke, L., et al. 2015, MNRAS, 449, 1352 [Google Scholar]
  32. Cowie, L. L., Songaila, A., Hu, E. M., & Cohen, J. G. 1996, AJ, 112, 839 [Google Scholar]
  33. Cowley, W. I., Caputi, K. I., Deshmukh, S., et al. 2018, ApJ, 853, 69 [NASA ADS] [CrossRef] [Google Scholar]
  34. Cowley, W. I., Caputi, K. I., Deshmukh, S., et al. 2019, ApJ, 874, 114 [NASA ADS] [CrossRef] [Google Scholar]
  35. Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 450, 1937 [Google Scholar]
  36. Croton, D. 2013, PASA, 30 [CrossRef] [Google Scholar]
  37. Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11 [Google Scholar]
  38. Dalla Vecchia, C., & Schaye, J. 2012, MNRAS, 426, 140 [Google Scholar]
  39. Dashyan, G., & Dubois, Y. 2020, A&A, 638, A123 [EDP Sciences] [Google Scholar]
  40. Dashyan, G., Choi, E., Somerville, R. S., et al. 2019, MNRAS, 487, 5889 [NASA ADS] [CrossRef] [Google Scholar]
  41. Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70 [Google Scholar]
  42. Davis, M., Efstathiou, G., Frenk, C. S., & White, S. D. M. 1985, ApJ, 292, 371 [Google Scholar]
  43. De Lucia, G., Springel, V., White, S. D. M., Croton, D., & Kauffmann, G. 2006, MNRAS, 366, 499 [NASA ADS] [CrossRef] [Google Scholar]
  44. Dekel, A., & Silk, J. 1986, ApJ, 303, 39 [Google Scholar]
  45. Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Nature, 457, 451 [Google Scholar]
  46. Desmond, H., Mao, Y.-Y., Wechsler, R. H., Crain, R. A., & Schaye, J. 2017, MNRAS, 471, L11 [NASA ADS] [CrossRef] [Google Scholar]
  47. Despali, G., Giocoli, C., Angulo, R. E., et al. 2016, MNRAS, 456, 2486 [Google Scholar]
  48. Diemer, B. 2018, ApJS, 239, 35 [Google Scholar]
  49. Donnari, M., Pillepich, A., Nelson, D., et al. 2021, MNRAS, 506, 4760 [NASA ADS] [CrossRef] [Google Scholar]
  50. Dubois, Y., & Teyssier, R. 2008, A&A, 477, 79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2010, MNRAS, 409, 985 [Google Scholar]
  52. Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2012, MNRAS, 420, 2662 [Google Scholar]
  53. Dubois, Y., Gavazzi, R., Peirani, S., & Silk, J. 2013, MNRAS, 433, 3297 [CrossRef] [Google Scholar]
  54. Dubois, Y., Pichon, C., Welker, C., et al. 2014, MNRAS, 444, 1453 [Google Scholar]
  55. Dubois, Y., Volonteri, M., Silk, J., et al. 2015, MNRAS, 452, 1502 [Google Scholar]
  56. Dubois, Y., Beckmann, R., Bournaud, F., et al. 2021, A&A, 651, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Duffy, A. R., Schaye, J., Kay, S. T., & Vecchia, C. D. 2008, MNRAS, 390, L64 [NASA ADS] [CrossRef] [Google Scholar]
  58. Engler, C., Pillepich, A., Joshi, G. D., et al. 2020, MNRAS, 500, 3957 [NASA ADS] [CrossRef] [Google Scholar]
  59. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
  60. Gabor, J. M., & Davé, R. 2015, MNRAS, 447, 374 [NASA ADS] [CrossRef] [Google Scholar]
  61. Gouin, C., Gavazzi, R., Pichon, C., et al. 2019, A&A, 626, A72 [Google Scholar]
  62. Groth, E. J., & Peebles, P. J. E. 1977, ApJ, 217, 385 [Google Scholar]
  63. Gunn, J. E., & Gott, J. R., III 1972, ApJ, 176, 1 [Google Scholar]
  64. Habouzit, M., Volonteri, M., & Dubois, Y. 2017, MNRAS, 468, 3935 [NASA ADS] [CrossRef] [Google Scholar]
  65. Hatfield, P. W., Laigle, C., Jarvis, M. J., et al. 2019, MNRAS, 490, 5043 [NASA ADS] [CrossRef] [Google Scholar]
  66. Hopkins, P. F., Quataert, E., & Murray, N. 2012, MNRAS, 421, 3522 [Google Scholar]
  67. Hopkins, P.F., Kere, D., Oñorbe, J., et al. 2014, MNRAS, 445, 581 [NASA ADS] [CrossRef] [Google Scholar]
  68. Hopkins, P. F., Grudić, M. Y., Wetzel, A., et al. 2020, MNRAS, 491, 3702 [CrossRef] [Google Scholar]
  69. Ilbert, O., Arnouts, S., McCracken, H. J., et al. 2006, A&A, 457, 841 [Google Scholar]
  70. Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55 [Google Scholar]
  71. Ilbert, O., de la Torre, S., Martinet, N., et al. 2021, A&A, 647, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Ishikawa, S., Kashikawa, N., Tanaka, M., et al. 2020, ApJ, 904, 128 [NASA ADS] [CrossRef] [Google Scholar]
  73. Jarvis, M., Bernstein, G., & Jain, B. 2004, MNRAS, 352, 338 [Google Scholar]
  74. Joachimi, B., & Bridle, S. L. 2010, A&A, 523, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Jose, C., Baugh, C. M., Lacey, C. G., & Subramanian, K. 2017, MNRAS, 469, 4428 [NASA ADS] [CrossRef] [Google Scholar]
  76. Kaiser, N. 1992, ApJ, 388, 272 [Google Scholar]
  77. Katz, H., Ramsoy, M., Rosdahl, J., et al. 2020, MNRAS, 494, 2200 [NASA ADS] [CrossRef] [Google Scholar]
  78. Kauffmann, G., Colberg, J. M., Diaferio, A., & White, S. D. M. 1999, MNRAS, 307, 529 [NASA ADS] [CrossRef] [Google Scholar]
  79. Kaviraj, S., Laigle, C., Kimm, T., et al. 2017, MNRAS, 467, 4739 [NASA ADS] [Google Scholar]
  80. Kennicutt, R. C. Jr. 1998, ApJ, 498, 541 [Google Scholar]
  81. Kimm, T., Cen, R., Devriendt, J., Dubois, Y., & Slyz, A. 2015, MNRAS, 451, 2900 [CrossRef] [Google Scholar]
  82. Koekemoer, A. M., Aussel, H., Calzetti, D., et al. 2007, ApJS, 172, 196 [Google Scholar]
  83. Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18 [Google Scholar]
  84. Kovač, K., Lilly, S. J., Cucciati, O., et al. 2010, ApJ, 708, 505 [CrossRef] [Google Scholar]
  85. Kravtsov, A. V., Berlind, A. A., Wechsler, R. H., et al. 2004, ApJ, 609, 35 [Google Scholar]
  86. Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24 [Google Scholar]
  87. Laigle, C., Davidzon, I., Ilbert, O., et al. 2019, MNRAS, 486, 5104 [Google Scholar]
  88. Landy, S. D., & Szalay, A. S. 1993, ApJ, 412, 64 [Google Scholar]
  89. Larson, R. B., Tinsley, B. M., & Caldwell, C. N. 1980, ApJ, 237, 692 [Google Scholar]
  90. Leauthaud, A., Tinker, J., Behroozi, P. S., Busha, M. T., & Wechsler, R. 2011, ApJ, 738, 45 [NASA ADS] [CrossRef] [Google Scholar]
  91. Leauthaud, A., Tinker, J., Bundy, K., et al. 2012, ApJ, 744, 159 [Google Scholar]
  92. Legrand, L., McCracken, H. J., Davidzon, I., et al. 2019, MNRAS, 486, 5468 [Google Scholar]
  93. Licquia, T. C., & Newman, J. A. 2015, ApJ, 806, 96 [Google Scholar]
  94. Limber, D. N. 1953, ApJ, 117, 134 [Google Scholar]
  95. Mandelbaum, R., Seljak, U., Kauffmann, G., Hirata, C. M., & Brinkmann, J. 2006, MNRAS, 368, 715 [Google Scholar]
  96. Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113 [NASA ADS] [Google Scholar]
  97. Marinoni, C., & Hudson, M. J. 2002, ApJ, 569, 101 [NASA ADS] [CrossRef] [Google Scholar]
  98. Martín-Navarro, I., Burchett, J. N., & Mezcua, M. 2019, ApJ, 884, L45 [CrossRef] [Google Scholar]
  99. Matthee, J., Schaye, J., Crain, R. A., et al. 2017, MNRAS, 465, 2381 [NASA ADS] [CrossRef] [Google Scholar]
  100. McAlpine, S., Helly, J. C., Schaller, M., et al. 2016, Astron. Comput., 15, 72 [Google Scholar]
  101. McCracken, H. J., Peacock, J. A., Guzzo, L., et al. 2007, ApJS, 172, 314 [NASA ADS] [CrossRef] [Google Scholar]
  102. McCracken, H. J., Milvang-Jensen, B., Dunlop, J., et al. 2012, A&A, 544, A156 [Google Scholar]
  103. McCracken, H. J., Wolk, M., Colombi, S., et al. 2015, MNRAS, 449, 901 [NASA ADS] [CrossRef] [Google Scholar]
  104. Mead, A. J., Peacock, J. A., Heymans, C., Joudaki, S., & Heavens, A. F. 2015, MNRAS, 454, 1958 [Google Scholar]
  105. Mead, A. J., Heymans, C., Lombriser, L., et al. 2016, MNRAS, 459, 1468 [Google Scholar]
  106. Mead, A., Brieden, S., Tröster, T., & Heymans, C. 2021, MNRAS, 502, 1401 [Google Scholar]
  107. Meneux, B., Guzzo, L., de la Torre, S., et al. 2009, A&A, 505, 463 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  108. Moneti, A., McCracken, H. J., Shuntov, M., et al. 2022, A&A, 658, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  109. Moore, B., Katz, N., Lake, G., Dressler, A., & Oemler, A. 1996, Nature, 379, 613 [Google Scholar]
  110. Moster, B. P., Somerville, R. S., Maulbetsch, C., et al. 2010, ApJ, 710, 903 [Google Scholar]
  111. Moster, B. P., Somerville, R. S., Newman, J. A., & Rix, H.-W. 2011, ApJ, 731, 113 [Google Scholar]
  112. Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121 [Google Scholar]
  113. Moster, B. P., Naab, T., & White, S. D. M. 2018, MNRAS, 477, 1822 [Google Scholar]
  114. Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206 [Google Scholar]
  115. Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [Google Scholar]
  116. Nayyeri, H., Hemmati, S., Mobasher, B., et al. 2017, ApJS, 228, 7 [NASA ADS] [CrossRef] [Google Scholar]
  117. Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624 [Google Scholar]
  118. Nelson, D., Springel, V., Pillepich, A., et al. 2019, Comput. Astrophys. Cosmol., 6, 2 [Google Scholar]
  119. Nuñez-Castiñeyra, A., Nezri, E., Devriendt, J., & Teyssier, R. 2021, MNRAS, 501, 62 [Google Scholar]
  120. Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713 [Google Scholar]
  121. Oser, L., Ostriker, J. P., Naab, T., Johansson, P. H., & Burkert, A. 2010, ApJ, 725, 2312 [Google Scholar]
  122. Peacock, J. A., & Smith, R. E. 2000, MNRAS, 318, 1144 [Google Scholar]
  123. Peirani, S., Dubois, Y., Volonteri, M., et al. 2017, MNRAS, 472, 2153 [NASA ADS] [CrossRef] [Google Scholar]
  124. Peng, Y., Lilly, S. J., Kovac, K., et al. 2010, ApJ, 721, 193 [NASA ADS] [CrossRef] [Google Scholar]
  125. Pillepich, A., Nelson, D., Hernquist, L., et al. 2018a, MNRAS, 475, 648 [Google Scholar]
  126. Pillepich, A., Springel, V., Nelson, D., et al. 2018b, MNRAS, 473, 4077 [Google Scholar]
  127. Planck Collaboration XVI. 2014, A&A, 571, A16 [Google Scholar]
  128. Planck Collaboration XIII. 2016, A&A, 594, A13 [Google Scholar]
  129. Pozzetti, L., Bolzonella, M., Zucca, E., et al. 2010, A&A, 523, A13 [Google Scholar]
  130. Prgomet, M., Rey, M. P., Andersson, E. P., et al. 2022, MNRAS, 513, 2326 [CrossRef] [Google Scholar]
  131. Roche, N., & Eales, S. A. 1999, MNRAS, 307, 703 [NASA ADS] [CrossRef] [Google Scholar]
  132. Rosdahl, J., Blaizot, J., Aubert, D., Stranex, T., & Teyssier, R. 2013, MNRAS, 436, 2188 [Google Scholar]
  133. Rosdahl, J., Schaye, J., Dubois, Y., Kimm, T., & Teyssier, R. 2017, MNRAS, 466, 11 [NASA ADS] [CrossRef] [Google Scholar]
  134. Salem, M., & Bryan, G. L. 2014, MNRAS, 437, 3312 [CrossRef] [Google Scholar]
  135. Salpeter, E. E. 1955, ApJ, 121, 161 [Google Scholar]
  136. Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521 [Google Scholar]
  137. Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103 [Google Scholar]
  138. Schmidt, M. 1968, ApJ, 151, 393 [Google Scholar]
  139. Scoccimarro, R., Sheth, R. K., Hui, L., & Jain, B. 2001, ApJ, 546, 20 [NASA ADS] [CrossRef] [Google Scholar]
  140. Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1 [Google Scholar]
  141. Seljak, U. 2000, MNRAS, 318, 203 [Google Scholar]
  142. Shankar, F., Lapi, A., Salucci, P., De Zotti, G., & Danese, L. 2006, ApJ, 643, 14 [NASA ADS] [CrossRef] [Google Scholar]
  143. Silk, J., & Mamon, G. A. 2012, Res. Astron. Astrophys., 12, 917 [Google Scholar]
  144. Silk, J., Di Cintio, A., & Dvorkin, I. 2014, Proc. Int. School Phys. Enrico Fermi, 186, 137 [NASA ADS] [Google Scholar]
  145. Somerville, R. S., & Davé, R. 2015, ARA&A, 53, 51 [Google Scholar]
  146. Song, M., Finkelstein, S. L., Ashby, M. L. N., et al. 2016, ApJ, 825, 5 [NASA ADS] [CrossRef] [Google Scholar]
  147. Springel, V. 2005, MNRAS, 364, 1105 [Google Scholar]
  148. Springel, V. 2010, ARA&A, 48, 391 [Google Scholar]
  149. Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289 [Google Scholar]
  150. Springel, V., White, S. D. M., Tormen, G., & Kauffmann, G. 2001, MNRAS, 328, 726 [Google Scholar]
  151. Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676 [Google Scholar]
  152. Steinhardt, C. L., Jespersen, C. K., & Linzer, N. B. 2021, ApJ, 923, 8 [NASA ADS] [CrossRef] [Google Scholar]
  153. Szalay, A. S., Connolly, A. J., & Szokoly, G. P. 1999, AJ, 117, 68 [Google Scholar]
  154. Tasitsiomi, A., Kravtsov, A. V., Wechsler, R. H., & Primack, J. R. 2004, ApJ, 614, 533 [NASA ADS] [CrossRef] [Google Scholar]
  155. Teyssier, R. 2002, A&A, 385, 337 [Google Scholar]
  156. Thomas, D., Maraston, C., Schawinski, K., Sarzi, M., & Silk, J. 2010, MNRAS, 404, 1775 [NASA ADS] [Google Scholar]
  157. Tinker, J. L., Robertson, B. E., Kravtsov, A. V., et al. 2010, ApJ, 724, 878 [Google Scholar]
  158. Tinker, J. L., Leauthaud, A., Bundy, K., et al. 2013, ApJ, 778, 93 [Google Scholar]
  159. Tweed, D., Devriendt, J., Blaizot, J., Colombi, S., & Slyz, A. 2009, A&A, 506, 647 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  160. Vale, A., & Ostriker, J. P. 2004, MNRAS, 353, 189 [Google Scholar]
  161. Vardoulaki, E., Gozaliasl, G., Finoguenov, A., Jiménez-Andrade, E. F., & Team T. C. 2021, ArXiv e-prints [arXiv:2104.04288] [Google Scholar]
  162. Vogelsberger, M., Genel, S., Sijacki, D., et al. 2013, MNRAS, 436, 3031 [Google Scholar]
  163. Vogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518 [Google Scholar]
  164. Weaver, J. R., Kauffmann, O. B., Ilbert, O., et al. 2022, ApJS, 258, 11 [Google Scholar]
  165. Wechsler, R. H., & Tinker, J. L. 2018, ARA&A, 56, 435 [Google Scholar]
  166. White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341 [Google Scholar]
  167. Xu, H., Zheng, Z., Guo, H., Zhu, J., & Zehavi, I. 2016, MNRAS, 460, 3647 [NASA ADS] [CrossRef] [Google Scholar]
  168. Zehavi, I., Weinberg, D. H., Zheng, Z., et al. 2004, ApJ, 608, 16 [NASA ADS] [CrossRef] [Google Scholar]
  169. Zehavi, I., Zheng, Z., Weinberg, D. H., et al. 2011, ApJ, 736, 59 [Google Scholar]
  170. Zheng, Z., Berlind, A. A., Weinberg, D. H., et al. 2005, ApJ, 633, 791 [NASA ADS] [CrossRef] [Google Scholar]
  171. Zheng, Z., Coil, A. L., & Zehavi, I. 2007, ApJ, 667, 760 [Google Scholar]

Appendix A: Impact of N(z) on w(θ)

As discussed in Section 2.2, the assumed redshift distribution of the sample can lead to differences in the modeled w(θ) that are considerably larger than the measurement uncertainties that may bias the inferred model parameters. We acknowledge, following the discussion in Ilbert et al. (2021), that the N(z) constructed by stacking the posteriors 𝒫(z) for each individual galaxy is more representative of the true underlying distribution of the sample, compared to simply stacking the likelihoods ℒ(z) as outputted from the template fitting code LePhare. In Fig. A.1, we show the N(z) for the sample in the bin 0.8 <  zphot <  1.1 obtained using ℒ(z) in red, 𝒫(z) in blue, and the histogram of the sources that have a spectroscopic redshift available in green. The figure shows that indeed stacking the individual 𝒫(z) agrees better with the spec-z histogram, especially in the tails which contribute the most in the w(θ) amplitude. Figure A.2 shows the model w(θ) of the two N(z) (upper panel). The relative difference between the two is about 35% as shown in the bottom panel. Furthermore, we explore the effect of a bias in the mean redshift as inferred from the N(z) on the w(θ). For this we use a Gaussian distribution centered at the mean redshift of the bin and then move the mean by −0.02. This results in relative difference in less than 3% which is safely within the error bars of the measurements (which are about 10%).

thumbnail Fig. A.1.

Redshift distributions based on spec-z histogram and stacking of photo-z likelihood and posterior distributions. All distributions are normalized to the maximum of the spec-z histogram.

thumbnail Fig. A.2.

Effect of the redshift distribution on the model correlation function

Appendix B: Details of the HOD-model derivation of the correlation functions

B.1. Clustering correlation function

In brief, one central component of the model is the galaxy-galaxy power spectrum, which can be separated in contributions from clustering of galaxies within the same halo (one-halo term) and between different halos (two-halo term):

(B.1)

These two components of the power spectrum of galaxies can be modeled under the HOD framework (Cooray & Sheth 2002; Berlind & Weinberg 2002) where the one-halo and two-halo terms are given by:

(B.2)

In these equations, is the mean number density of galaxies, bh(Mh, z) is the large-scale halo bias which we chose the one given by Tinker et al. (2010). us(k) is the Fourier transform of the over-density profile of satellite galaxies, for which we assume that it follows the NFW profile Navarro et al. (1997) with a mass-concentration relation as calibrated by Duffy et al. (2008). The important thing to note here is that the power spectrum is defined in terms of the occupation distributions of centrals and satellites ⟨Nc⟩ and ⟨Ns⟩ specified by Eqs. (10) and (11) and is where our parametrization of the SHMR enters the model. The linear power spectrum Plin(k, z) enters the two-halo term and dominates large scales. The halo model-based power spectrum is known to be inaccurate on quasi-linear scales at about k ∼ 0.5 Mpc−1 at z = 0 which falls at the transition between the one- and two-halo terms (Mead et al. 2015, 2016; Chisari et al. 2019). We correct this by following Eq. (23) of Mead et al. (2021), where a time-dependent function α(z) is applied to the one- and two-halo power spectrum terms.

Having the galaxy-galaxy power spectrum, one can then compute the angular power spectrum via the Limber equation (Limber 1953; Kaiser 1992; Joachimi & Bridle 2010)

(B.3)

where H(z) is the expansion rate at redshift z, χ(z) is the radial comoving distance at z, and Ni(z) and Nj(z) are the redshift distribution of the galaxies in the samples i and j. Since we’re interested in the angular auto-correlation function, both i and j are the same sample. Finally, to arrive at the angular correlation function w(θ), we perform an inverse Fourier transform of C numerically through the Hankel transforms under the flat sky approximation:

(B.4)

where J0(θ ℓ) is the 0-th order Bessel function (for more details see Chisari et al. 2019).

Appendix C: Description of the hydrodynamical simulations

C.1. Main characteristics and catalog extraction

HORIZON-AGN has been produced with the adaptive mesh refinement code RAMSES (Teyssier 2002), using WMAP7 cosmological parameters (Komatsu et al. 2011) in a box of size 100 h−1 Mpc a side 12 and a DM mass resolution of MDM, res = 8 × 107M. The simulation includes gas heating, cooling, star-formation, feedback from stellar winds, type Ia and type II supernovae with mass, energy, and metal release in the interstellar-medium. The simulation also follows the formation of black holes and energy release from AGN in a quasar or radio mode depending on the accretion rate. Full details on the subgrid implementation are given in Dubois et al. (2014).

Galaxies and halos are identified using the ADAPTAHOP structure finder (Aubert et al. 2004; Tweed et al. 2009) applied to the distribution of star and DM particles respectively, as described in previous works (e.g., Dubois et al. 2014). Galaxy masses are obtained by summing the masses of all individual stellar particles within twice the effective radius of the galaxies, while halo masses are obtained by summing all DM particles within the halos. Galaxies and halos are positionally matched. The central galaxies of a halo is defined as the most massive galaxies within 0.1 times the halo virial radius.

We note that stellar mass losses in HORIZON-AGN were implemented assuming a Salpeter IMF (Salpeter 1955). This can lead to ∼0.1 dex more stellar masses at later time than when assuming a Chabrier IMF (Chabrier 2003). The simulated masses were therefore rescaled accordingly for a consistent comparison with both the other simulations and the observational data (which assumes a Chabrier IMF when deriving the mass from the photometry).

The simulation is in relatively good agreement with observations up to z = 4. Known discrepancies include the overestimation of galaxy masses, especially at low mass (Kaviraj et al. 2017). These low-mass galaxies are on overall too quenched, an indication that star formation at high redshift was not regulated enough and the gas in the environment of these galaxies was too rapidly consumed. In addition, the bimodality is not as well-marked as in observations due to residual star-formation in massive galaxies, possibly because of a slightly insufficient strength for the AGN feedback.

TNG100-1 has been produced with the moving-mesh code AREPO (Springel 2010) using cosmological parameters from Planck Collaboration XIII (2016) in a box of size 75 h−1 Mpc a side and a DM mass resolution of MDM, res = 7.5 × 106M. It follows magnetic fields in addition to the hydrodynamical processes above-mentioned for the HORIZON-AGN simulation (with different subgrid implementations, see the discussion below). Full details on the subgrid implementation are given in Pillepich et al. (2018b), Springel et al. (2018).

Halos, subhalos, and galaxies have been identified using the FoF (Davis et al. 1985) and SUBFIND (which, within a halo identified with the FoF technique, relies on all particle species to identify the galaxy, see Springel et al. 2001) algorithms, as described in, for instance, Pillepich et al. (2018a). We downloaded the halo and galaxy catalogs from the public website13. In the following, galaxy masses are estimated by summing all stellar mass particles within twice the stellar effective radius, while halo masses are taken as being the sum of all individual dark matter particles in the identified halos, which matches the definition of halo and galaxy masses in HORIZON-AGN. Galaxies that do not reside within R200 of a larger halo are identified as centrals. Galaxy clustering has been measured in Springel et al. (2018) and the SHMR has been presented in Pillepich et al. 2018a (see also Engler et al. 2020). In particular, these authors investigated how the chosen estimate for stellar mass modifies the measured relation. While we are aware of this discussion, we chose in this work to measure stellar mass within twice the effective radius for consistency between different simulated datasets. Such measurement is also closer to our observational estimate of stellar masses (which derive from SED-fitting on the total magnitudes – and not on aperture magnitudes).

EAGLE has been produced using a modified version of the N-Body Tree-PM smoothed particle hydrodynamics (SPH) code GADGET-3 (Springel 2005), adopting cosmological parameters from Planck (Planck Collaboration XVI 2014). In this work, we use the reference run called Ref-L0100N1504,corresponding to a box of 100 comoving Mpc a side with a DM mass resolution of MDM, res = 9.7 × 106M. EAGLE follows gas heating and cooling, star-formation, feedback from stellar winds and AGB stars, along with type Ia and type II supernovae and from AGN. Full details on the subgrid implementation are given in Crain et al. (2015) and Schaye et al. (2015).

As for TNG100-1, the halos, subhalos, and galaxies have been identified using the FoF (Davis et al. 1985) and SUBFIND (Springel et al. 2001) algorithms, as described in McAlpine et al. (2016). We downloaded halo and galaxy catalogs from the public website14. Galaxy and halo masses were obtained by summing the masses of all stellar and dark matter particles, respectively, which are part of the objects. The central galaxy is taken as the one which contains the particle with the lowest value of the gravitational potential. The redshift evolution of the SHMR for central galaxies in EAGLE simulation was briefly presented in Matthee et al. (2017).

C.2. Implementation of subgrid recipes

Our comparisons of the SHMR in observations and simulations has highlighted two main discrepancies: (1) in HORIZON-AGN, the SHMR is always overestimated with respect to observations and the two other simulations, especially at low masses; (2) in all simulations, the relative contribution of satellites with respect to centrals is too large compared to observations.

With regard to (1), we note first that in EAGLE, the free parameters of the subgrid (stellar and AGN) feedback model were calibrated so that the simulations reproduce the galaxy stellar mass function and SHMR, galaxy sizes, and the empirical relation between black hole mass and stellar mass all at z = 0 (Schaye et al. 2015). In TNG100-1, the calibration was performed against the star formation rate density as a function of cosmic time and the stellar mass function and SHMR both at z = 0 (Pillepich et al. 2018b). In HORIZON-AGN however, the calibration was less constraining since only the efficiency of AGN feedback was tuned so that the black hole–galaxy scaling relation at z = 0 was reproduced (see Dubois et al. 2012, for details). It is therefore not surprising that both EAGLE and TNG100-1 better reproduce, overall, the stellar mass function (Schaye et al. 2015; Pillepich et al. 2018a) and SHMR (Fig. 12) at low masses, because their calibration specifically constrain stellar feedback, which is not the case for HORIZON-AGN.

Without pretending to exhaustively discuss the differences in stellar feedback and star-formation implementation between HORIZON-AGN and the two other simulations, we instead opt to highlight a few aspects. In the EAGLE simulation (as explained in Crain et al. 2015), the amount of injected energy from feedback depends on the local properties of the gas (it decreases with metallicity and increases with gas density), the calibration being adjusted in order to reproduce the local stellar mass function. This tuning contributes to increase the supernova feedback energy at high redshift beyond the energy available for their adopted Chabrier IMF (Crain et al. 2015). Furthermore, the stochasticity of energy deposition of supernovae is artificially increased to enhance the impact of their feedback in terms of wind mass loading and quenching of star formation15 (Dalla Vecchia & Schaye 2012). Finally, the EAGLE star-formation law does not follow the standard Kennicutt-Schmidt prescription (Kennicutt 1998, as adopted in HORIZON-AGN); rather, it depends on pressure instead of density and includes a metallicity-dependent density threshold (against a simple density threshold in HORIZON-AGN), which tends to reduce star-formation in metal-rich regions. Those features are likely to contribute to a more efficient quenching of star-formation in EAGLE compared to HORIZON-AGN.

Star-formation in TNG100 follows the empirical Kennicutt-Schmidt relation, however, feedback implementation is sensibly different, as described in Pillepich et al. (2018b). More specifically, the energy transfer from SNe to large-scale galactic winds is very efficient because of the hydrodynamical decoupling of the launched wind gas from the dense star-forming gas, until they recouple hydrodynamically with the circum-galactic gas (Springel & Hernquist 2003; Vogelsberger et al. 2013). These authors used a wind velocity that is proportional to the local DM velocity dispersion and cosmic time so that the winds are faster in more massive halos and at lower redshift. Similarly to the EAGLE model, the given supernovae energy to the gas is higher for energy deposit in lower metallicity gas, with the same scaling with metallicity as in EAGLE. Taken together, these features contribute to make stellar feedback more effective at high redshift in TNG100 compared to HORIZON-AGN.

Finally, we note that many missing mechanisms could be naturally added to the HORIZON-AGN model in order to increase the strength of stellar feedback in low-mass galaxies at high redshift in a physically motivated way (without necessarily requiring empirical tuning of parameters). These processes include radiation from stars (see the discussion below), cosmic ray feedback from supernovae (e.g., Booth et al. 2013; Salem & Bryan 2014; Dashyan & Dubois 2020), a gravo-turbulent model for star formation (e.g., Nuñez-Castiñeyra et al. 2021; Dubois et al. 2021), or adopting an IMF varying with redshift, stellar density, or metallicity (e.g., Applebaum et al. 2020; Prgomet et al. 2022).

With regard to (2), however, the efficient stellar feedback implemented in TNG100 and EAGLE do not prevent the too large satellite fraction (Fig. 8) and excessive contribution of satellites to the total SHMR (Fig. 9). As previously noted in Sec. 6.3, such a discrepancy must be due to a lack of satellite quenching at high redshift, making these satellites to grow too massive (with respect to observations) before being eventually quenched. Galaxy pre-processing by stellar radiation feedback (e.g., Rosdahl et al. 2013; Hopkins et al. 2020) resulting in less tightly-bound galaxies could be a prerequisite to efficiently quench satellites (once they have been accreted) through tidal stripping, as suggested by Costa et al. (2019).

Finally, some studies mention the possible role of AGN feedback from the central galaxy in quenching its satellites (e.g., Dashyan et al. 2019; Martín-Navarro et al. 2019): the lack of quenching in high-redshift satellites could therefore also be the consequence of an imperfect AGN feedback implementation at high redshift, although hydrodynamical simulations show that the activity of SNe in low mass galaxies quench the growth of black holes and their associated AGN feedback (e.g., Dubois et al. 2015; Habouzit et al. 2017; Anglés-Alcázar et al. 2017).

Appendix D: Data

The data used to make all the plots in this paper – measurements of the clustering and GSMF, best-fit models, derived quantities, and so on, are shared publicly with the community in tabular form on the following github repository: https://github.com/mShuntov/SHMR_COSMOS2020

Appendix E: Posterior probabilities of the parameters

thumbnail Fig. E.1.

Posterior distributions for each redshift fit

thumbnail Fig. E.2.

Posterior distributions for each redshift fit

Appendix F: Best-fit values of the model parameters

Table F.1.

Best-fit values of the model parameters in the ten redshift bins

All Tables

Table 1.

Sample selection in redshift and stellar mass thresholds.

Table 2.

Adopted ingredients in the halo model.

Table F.1.

Best-fit values of the model parameters in the ten redshift bins

All Figures

thumbnail Fig. 1.

Sample selection in the stellar mass-redshift plane. The solid grid lines show the mass threshold for each sample. The solid violet curve indicates the stellar mass completeness limit. We note that the histogram does not correspond to the final selection since further selection criteria (e.g., based on PDF(z) width) are applied.

In the text
thumbnail Fig. 2.

Redshift distribution of the ten galaxy samples used for the clustering and abundance measurements. The redshift distribution is obtained by stacking the posterior photo-z distributions for all the sources in a given bin as described in Sect. 2.2.

In the text
thumbnail Fig. 3.

Best-fit models for clustering and abundance compared with measurements at 0.5 <  z <  0.8. Top: clustering measurements for the 6 mass threshold-selected samples (empty circles with errorbars) along with the best-fit models in solid lines in corresponding colors. Bottom: measurements of the stellar mass function together with the best fit model. The dashed lines show the SMF in the same redshift bin obtained by Davidzon et al. (2017) for comparison.

In the text
thumbnail Fig. 4.

Best-fit models of clustering and abundance plotted over the measurements for all z-bins apart from 0.5 <  z <  0.8, which is shown in Fig. 3.

In the text
thumbnail Fig. 5.

Correlation of galaxies with M* >  1010M (left panel) and GSMF (right panel) for all ten redshift bins. The green dashed lines in the right panel correspond to the GSMF of Davidzon et al. (2017).

In the text
thumbnail Fig. 6.

Mean number of galaxies with stellar masses in a given mass bin as a function of the mass of the halo that they occupy. We show the mean halo occupation function for galaxies in 4 stellar mass bins (color coded accordingly) at 0.2 <  z <  0.5. The thick solid lines show the total ⟨Ntot⟩ and the dashed and dotted lines show the centrals and satellites.

In the text
thumbnail Fig. 7.

Mean halo occupation in halos of a given mass as a function of redshift. In each panel, we show the ⟨Ntot(Mh)⟩ at log Mh/M = [12.0,  12.5. 13.0]. The three different panels show the mean occupations by galaxies in three different stellar mass bins log M*/M = {[9.0,  9.5], [9.5,  10.0], [10.0,  10.5]}. The dashed and dotted lines show the central and satellite occupations, while the points connected with transparent solid line show the total.

In the text
thumbnail Fig. 8.

Fraction of satellite galaxies with masses above a given threshold as a function of redshift (solid lines). In each of the three panels, we compare with the satellite fractions measured in the hydrodynamical simulations (dashed lines) TNG100 (left), HORIZON-AGN (center), and EAGLE (right). The redshift evolution of the satellite fraction is shown for galaxies with masses above log M*/M >  [8.7,  9.0,  9.3,  9.7,  10.2,  10.7,  11.0].

In the text
thumbnail Fig. 9.

Stellar-to-halo mass relation (top) and M*/Mh ratio (bottom) in the ten redshift bins. The solid lines and shaded regions show our inferred SHMR and 1σ confidence interval color-coded according to the redshift bin.

In the text
thumbnail Fig. 10.

Evolution of the peak halo (left panel) and peak stellar mass (right panel) with redshift. The results from our analysis are shown in purple points. For comparison, we show measurements from the literature rescaled to match our chosen value for H0 = 70 km s−1 Mpc−1. The literature measurements include Legrand et al. (2019), Leauthaud et al. (2012, L+12), Coupon et al. (2012, C+12), Coupon et al. (2015, C+15), Cowley et al. (2018, C+18), Moster et al. (2013, M+13), Behroozi et al. (2013, B+13), Behroozi et al. (2019, B+19), and from the hydrodynamic simulations HORIZON-AGN, TNG100, and EAGLE (references in the main text).

In the text
thumbnail Fig. 11.

Total stellar-to-halo mass relation (top panel) and total M*/Mh (bottom panel) at 0.2 <  z <  0.5 compared to hydrodynamical simulations. The purple dashed and dotted lines show our central and satellite contribution to the total SHMR respectively, with the shaded region showing 1σ confidence interval of the sum of the two. The break in solid purple lines and shaded regions indicate the highest stellar mass probed in our analysis, which we take to be the highest mass bin in the SMF. The transparent purple lines and hatched region is an extrapolation at higher masses. We overplot the SHMRs measured in the hydrodynamical simulations HORIZON-AGN in teal, TNG100 in red, and EAGLE in dark yellow, where the dashed dotted and solid lines show the central, satellite, and total SHMR.

In the text
thumbnail Fig. 12.

Total M*/Mh in the different redshift bins compared to hydrodynamical simulations. The purple lines show our inferred central and satellite contribution to the total SHMR with the shaded region showing a 1σ confidence interval of the sum of the two. The break in solid purple lines and shaded regions indicate the highest stellar mass probed in our analysis, which we take to be the highest mass bin in the SMF. The dashed purple lines and hatched region is an extrapolation at higher masses. The comparison includes total M*/Mh found in the hydrodynamical simulations HORIZON-AGN in teal, TNG100 in red and EAGLE in dark yellow. The dashed and dotted lines for the simulations indicate the central and satellite contributions, respectively.

In the text
thumbnail Fig. 13.

Redshift dependence of M*/Mh at fixed halo mass. We show the total (centrals and satellites) M*/Mh fixed at log Mh/M = [11.50,  12.00,  13.00,  13.60, ], while the M*/Mh for log Mh/M = 13.60 is also shown for centrals and satellites separately with star and diamond symbols, respectively.

In the text
thumbnail Fig. 14.

Uncertainties in the SMF expressed as a fractional error σ/Φ. We show the uncertainties from cosmic variance with solid line-connected points and the uncertainties from SED fitting in dashed line-connected diamonds. We show the uncertainties for four redshift bins.

In the text
thumbnail Fig. A.1.

Redshift distributions based on spec-z histogram and stacking of photo-z likelihood and posterior distributions. All distributions are normalized to the maximum of the spec-z histogram.

In the text
thumbnail Fig. A.2.

Effect of the redshift distribution on the model correlation function

In the text
thumbnail Fig. E.1.

Posterior distributions for each redshift fit

In the text
thumbnail Fig. E.2.

Posterior distributions for each redshift fit

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.