Unveiling the (in)consistencies among the galaxy stellar mass function, star formation histories, satellite abundances, and intracluster light from a semi-empirical perspective

ABSTRACT

In a hierarchical, dark matter-dominated Universe, stellar mass functions (SMFs), galaxy merger rates, star formation histories (SFHs), satellite abundances, and intracluster light (ICL), should all be intimately connected observables. However, the systematics affecting observations still prevent universal and uniform measurements of, for example, the SMF and the SFHs, inevitably preventing theoretical models to compare with multiple data sets robustly and simultaneously. We here present our holistic semi-empirical model decode (Discrete statistical sEmi-empiriCal mODEl) that converts via abundance matching dark matter merger trees into galaxy assembly histories, using different SMFs in input and predicting all other observables in output in a fully data-driven and self-consistent fashion with minimal assumptions. We find that: (1) weakly evolving or nearly constant SMFs below the knee (⁠|$M_\star \lesssim 10^{11} \, \mathrm{M}_\odot$|⁠) are the best suited to generate SFHs aligned with those inferred from MaNGA, SDSS, GAMA, and, more recently, JWST; (2) the evolution of satellites after infall only affects the satellite abundances and SFHs of massive central galaxies but not their merger histories; (3) the resulting SFR–|$M_\star$| relation is lower in normalization by a factor of |$\sim 2$| with respect to observations, with a flattening at high masses more pronounced in the presence of mergers; (4) the latest data on ICL can be reproduced if mass-loss from mergers is included in the models. Our findings are pivotal in acting as pathfinder to test the self-consistency of the high-quality data from, e.g. JWST and Euclid.

1 INTRODUCTION

In a Lambda cold dark matter (⁠|$\Lambda$|CDM) scenario galaxies form within the potential well of dark matter haloes. They are believed to grow their stellar mass via both in situ and ex situ processes. In situ star formation originates from cooling of the infalling gas, while ex situ growth is predominantly driven by mergers with other galaxies. Several works have shown that most of the galaxies with |$M_\star \lt 10^{11}\, \mathrm{M}_\odot$| grow primarily via star formation, with mergers becoming increasingly relevant in more massive galaxies (e.g. van Dokkum et al. 2010; Leitner 2012; Shankar et al. 2015; Buchan & Shankar 2016; Groenewald et al. 2017; Matharu et al. 2019, Paper I; Eisert et al. 2023).

In traditional models of galaxy evolution, galaxies are evolved following the hierarchical growth of their host dark matter haloes, with a tendency to build up via mergers more massive systems at later epochs (e.g. Guo & White 2008; Oser et al. 2010; Cattaneo et al. 2011; Lackner et al. 2012; Rodriguez-Gomez et al. 2016; Pillepich et al. 2018b; Monachesi et al. 2019; Grylls et al. 2020, Paper I). However, in stellar archaeological findings and star formation histories (SFHs) analysis via spectral energy distribution (SED) fitting (e.g. Sánchez et al. 2019; Thomas et al. 2019; Bellstedt et al. 2020b), a different pattern emerges. Galaxies are observed to follow a downsizing trend, wherein more massive galaxies form at earlier times, possibly triggered by large bursts of star formation, while less massive galaxies form at later times with progressively lower star formation rates probably induced by an overall ‘cosmic starvation’ caused by the general reduction in available cold gas available to feed galaxies (e.g. Gunn & Gott 1972; Cowie & Songaila 1977; Larson, Tinsley & Caldwell 1980). Many galaxies, especially those with a prominent bulge component (e.g. Bell 2008; Bluck et al. 2014, 2022; Lang et al. 2014; Dimauro et al. 2022), tend to show signs of a more rapid quenching of their star formation. The primary physical processes responsible for rapid quenching are still unclear, with leading theories proposing a combination of internal and external processes, such as stellar and AGN feedback, halo quenching, morphological quenching, mergers, ram-pressure stripping (e.g. Granato et al. 2004; Dekel & Birnboim 2006, 2008; Shankar et al. 2006; Martig et al. 2009; Lilly et al. 2013; Schawinski et al. 2014; Lapi et al. 2018; Xu & Peng2021).

Many theoretical models and techniques have been developed in the last decades to study galaxy formation and evolution at different levels of detail, most notably, hydrodynamical simulations, semi-analytical models (SAMs), and semi-empirical models (SEMs). Hydrodynamical simulations, as a first example, are an extremely useful tool for studying galactic evolution, as they evolve dark matter and stellar particles simultaneously, allowing for a comprehensive overview of their co-evolution (Vogelsberger et al. 2014; Nelson et al. 2019). However, these simulations generally require high computational power and are affected by mass/volume resolution limitations. SAMs and SEMs, as complementary models, require a lower computational cost. SAMs typically initialize gas at early epochs and subsequently apply physical assumptions and parametrization for the latter (De Lucia et al. 2006; González et al. 2011; Guo et al. 2011; Cattaneo et al. 2020; Ayromlou et al. 2021b; Koutsouridou & Cattaneo2022)

, but could be susceptible to the degeneracy resulting from the high number of parameters. SEMs, instead, initialize galaxy stellar mass by matching the galactic properties (such as the stellar mass, star formation rate, or luminosity) to the dark matter halo properties (such as the halo mass, accretion rate, or peak velocity) via their relative abundances or via analytical parametrizations (e.g. Shankar et al. 2006, 2014; Hopkins et al. 2010a; Moster, Naab & White 2013, 2018; Tollet et al. 2017; Behroozi et al. 2019; Grylls et al. 2019; Boco et al. 2023), and evolve galaxies ensuring that at each redshift the main statistical observational properties of galaxies (e.g. star formation, number densities, and clustering) are reproduced. By design, SEMs strongly rely on observational data as input and, therefore, having a robust determination of the statistical properties of galaxies and scaling relations is fundamental for building a successful data-driven model.

The need for well calibrated, uniform, and self-consistent observational data sets is indeed of paramount importance for all types of theoretical galaxy evolution models, not just semi-empirical ones. Besides the quality of the data, it is fundamental to decipher the consistency of the different data sets used as term of comparison. Failures in simultaneously fitting distinct observational probes may be certainly ascribed to shortcomings in the underlying modelling, but they could also arise from disagreements in different data sets. For instance, the galaxy stellar mass function (SMF) contains essential information on the galaxy stellar mass growth. However, stellar masses are commonly inferred via SED fitting, which can introduce potential biases due to the SED fitting algorithm, observational bands, and cosmic variance. Indeed, several studies have reported different SMFs, suggesting various slopes at the bright end and evolution with redshift (e.g. Bernardi et al. 2013, 2016, 2017; Shankar et al. 2014; Davidzon et al. 2017; Kawinwanichakij et al. 2020; Weaver et al. 2023). Furthermore, observed star formation rates (SFRs) might not align consistently with the observed stellar masses. Typically, SFRs predicted by theoretical models and simulations tend to be lower than the observationally estimated values, leading to discrepancies in the observed stellar mass growths and mass functions (see e.g. Bernardi et al. 2010; Leja et al. 2015, 2022; Lapi et al. 2017; Grylls et al. 2020). For example, Leja et al. (2015) demonstrated that the observed SMF could not be reconciled with the observed main sequence. Grylls et al. (2019) showed that by using the observed SFRs (from e.g. Tomczak et al. 2016) the predicted satellites abundances are not compatible with those inferred from observations (e.g. Bernardi et al. 2017), and also showed that the predicted merger rates are significantly different when adopting different SMFs as inputs in SEMs (also see e.g. Hopkins et al. 2010b; O’Leary et al. 2021).

SEMs, by empirically linking different observables, prove to be extremely useful for testing possible inconsistencies among distinct data sets, which can often arise due to significant systematics in, e.g. measurements of stellar masses or star formation rates. For example, studying how galaxies build up their stellar mass across cosmic time will yield valuable insights into the properties and systematics of the SMF. In the previous paper of this series (Fu et al. 2022, hereafter referred to as Paper I) we presented decode, the Discrete statistical sEmi-empiriCal mODEl (DECODE) and we showed its suitability in addressing these tasks. In this work, we use decode to test the systematic inconsistencies that may be present among distinct observables, namely the SMF, the SFHs, the satellite abundances, and the intracluster light (ICL). To achieve this goal, we will use the most up to date estimates of the SFHs in galaxies from SED fitting, as well as the latest determinations of the SMF, of the satellite abundances and the ICL. We will show that there are still significant systematic inconsistencies among these data sets, which must be seriously considered when attempting to provide a holistic model of galaxy evolution.

This paper is structured as follows. In Section 2, we describe the data sets that we use in this work. In Section 3, we provide the details of our methodology and test the self-consistency of decode in predicting SFHs against hydrodynamical simulations. In Section 4, we present decode’s predictions for the SFHs of central galaxies, main sequence, satellite abundances, and ICL. In Sections 5 and 6, finally, we discuss the results that we found and draw our conclusions. In what follows we adopt the |$\Lambda$|CDM cosmological model with parameters from Planck Collaboration VI (2018) best-fitting values, i.e. |$(\Omega _{\rm m}, \Omega _{\Lambda }, \Omega _{\rm b}, h, n_{\rm S}, \sigma _8) = (0.31, 0.69, 0.049, 0.68, 0.97, 0.81)$|⁠, and we use a Chabrier (2003) stellar initial mass function.

2 DATA

In this work, we use (1) simulated data from the TNG100 simulation to validate our methodology and test the performance of decode for predicting SFHs; (2) the GalICS SAM to compare decode to an ab initio galaxy formation model; (3) observational data from SDSS (Sloan Digital Sky Survey), MaNGA (Mapping Nearby Galaxies at Apache Point Observatory), and GAMA (Galaxy And Mass Assembly) to test decode’s predictions for SFHs and satellite abundances. Below, we provide more details on each of these simulated and observed data sets.

2.1 The TNG simulation

In this study, we employ the TNG100 simulation, a component of the IllustrisTNG hydrodynamical simulation suite (TNG henceforth; Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018; Pillepich et al. 2018b; Springel et al. 2018). The TNG simulations are performed utilizing the moving-mesh arepo code (Springel 2010), which provides solutions for a combination of gravity and magnetohydrodynamical equations (Pakmor, Bauer & Springel 2011; Pakmor & Springel 2013). TNG facilitates the simulation of galaxies via subgrid modelling of critical galaxy formation processes, such as gas cooling, star formation, stellar evolution, feedback from stars, and supermassive black hole associated processes, including seeding, merging, and feedback (see Weinberger et al. 2017; Pillepich et al. 2018a for a comprehensive model description).

The TNG simulations encompass varying box sizes and mass resolutions. In this paper, we utilize the TNG100 simulation, executed within an approximately 100 Mpc box, which is the most suitable simulation for the objectives of our study. This approach enables us to circumvent any potential non-trivial resolution effects (see Pillepich et al. 2018a for further details).

For the identification of haloes, the TNG simulation employs the Friends of Friends algorithm (FOF; Davis et al. 1985), which locates a group of dark matter particles whose mutual distance is less than a linking length |$b=0.2$|⁠. Subsequently, the subhaloes are identified as gravitationally bound groups of particles, by executing the SUBFIND algorithm (Springel et al. 2001). Each FOF halo possesses a central subhalo, typically the most massive subhalo of the halo, with all other subhaloes identified as satellite subhaloes. Galaxies correspond to subhaloes with a non-zero stellar mass. FOF haloes generally lack a well-defined shape. However, it is customary to consider |$R_{200}$| (⁠|$R_{500}$|⁠), the radius within which the mean density is 200 (500) times the critical density of the Universe, as the halo radius.¹ Given that FOF haloes can extend beyond the |$R_{200}$|⁠, satellites of a FOF halo can exist and evolve both within and beyond the halo boundary (e.g. see Ayromlou et al. 2021a; Rohr et al. 2023).

In order to evaluate our abundance matching technique against the TNG simulation, we calculate the SMHM relation from the TNG data, by measuring the ratio between the total stellar mass of all galaxies (centrals and satellites) and the total mass within |$R_{200}$|⁠, for all haloes in the simulation.

We then employ subhalo merger trees, created using the sublink algorithm (Rodriguez-Gomez et al. 2015), to trace the evolution of individual galaxies over time. For each galaxy, we start from the present time (⁠|$z=0$|⁠, snapshot = 99) and trace back the evolution of its main progenitor through the 100 snapshots/redshifts of the TNG100 simulation. We continue this process until reaching the first appearance of the galaxy in the simulation. This approach enables us to extract several properties of all individual galaxies across cosmic time, including stellar mass growths and SFHs. We also segregate the stellar masses of galaxies into in situ, formed within the main progenitor branch of a galaxy, and ex situ, originating from minor and major mergers of the galaxy with other objects (see Rodriguez-Gomez et al. 2016). This allows the careful study of the evolution of galaxies, as presented in the subsequent sections.

2.2 GALICS

We make use of galics to compare the stellar mass assemblies and SFHs predicted by decode. galics (Galaxies In Cosmological Simulations; Hatton et al. 2003) is a semi-analytical model which follows the evolution of baryons in halo merger trees extracted from cosmological N-body simulations. In this work, we consider the galics 2.0 framework that includes information on dark matter substructures within each halo (Cattaneo et al. 2017). Central galaxies are at the centre of their host haloes. Satellite galaxies are associated with subhaloes. A formula based on hydrodynamic simulations by Jiang et al. (2008) is used to determine how long a satellite galaxy should take to merge with the central galaxy of its host system. As soon as a galaxy becomes a satellite, a merging countdown timer starts ticking.

In galics 2.1 (Cattaneo et al. 2020) there was introduced a physical criterion to determine when the gas that accretes on to a halo streams on to the central galaxy in cold flows and when it is shock-heated. In absence of feedback from active galactic nuclei (AGNs), the shock-heated gas in a halo or subhalo can cool on to the associated galaxy. Its cooling provides a second mechanism for gas accretion. Environmental effects such as ram pressure and tidal stripping reduce the importance of cooling in satellite galaxies (Koutsouridou & Cattaneo 2019, 2022, for details).

Gas accreted through cold accretion and cooling settles into discs, the sizes of which are determined by the conservation of angular momentum, and forms stars on a time-scale equal to 25 orbital times at one exponential scale-length from the galactic centre (Cattaneo et al. 2017). The conversion of gas into stars is much faster in merger-driven starbursts. Mergers also cause galaxies to grow bulges. The modelling of the effects of galaxy mergers is based on hydrodynamic simulations by Kannan et al. (2015) and described in Koutsouridou & Cattaneo (2022; GalICS 2.2). At each merger, a stellar mass equal to |$20{{\ \rm per\ cent}}$| of the stellar mass transferred to the bulge is scattered into the halo. We note that this assumption is very similar to the one made in decode for major mergers, which transfer most of the stars to the bulge component.

The main new feature of the GalICS 2.2 version is AGN feedback, described with the empirical model of Chen et al. (2020). Black holes grow and deposit feedback energy into the surrounding gas, until this feedback energy is larger than a fraction or multiple of its gravitational binding energy, at which point the gas is unbound or, more realistically, heated to high entropy, so that its cooling time becomes long. By doing so, AGN feedback quenches star formation and prevents its reactivation. The shutdown of star formation induced by AGN quenching is assumed to be permanent because when cooling restarts, it is self-limited. As soon as some gas cools, it reactivates the central AGN, which shuts down cooling again (maintenance mode).

2.3 Sloan Digital Sky Survey

To test the stellar mass function of satellites predicted by decode at |$z\sim 0.1$|⁠, we use the photometric catalogue of Meert, Vikram & Bernardi (2015, 2016) which is based on the Sloan Digital Sky Survey Data Release 7 (DR7; Abazajian et al. 2009) images. As also done in Paper I (section 2.4 therein), we use the Meert et al. catalogue of stellar masses derived from the best-fitting Sérsic-Exponential / Sérsic photometry on r-band observations, with mass-to-light ratios by Mendel et al. (2014). Furthermore, we adopt the truncation of the light profile as prescribed in Fischer, Bernardi & Meert (2017). The Meert et al. catalogues are matched with the Yang et al. (2007, 2012) group catalogues, which allow us to identify the satellite galaxies.

We also compare the stellar mass growths from the Mendel et al. (2014) sample with an average redshift of |$z\sim 0.1$|⁠, selecting 100 galaxies in each stellar mass bin of interest. To calculate the SFHs, we combine the MILES stellar population synthesis models of Vazdekis et al. (2010) with the Penalized PiXel-Fitting (pPXF) inversion algorithm. MILES models provide single stellar population (SSP) predictions in the optical range (⁠|$3540 {\small --} 7410$| Å) for a wide range of ages (0.03–14 Gyr) and metallicities (⁠|$-2.27 \lt [M/H] \lt 0.26$|⁠) at a 2.51 Å resolution based on the BaSTI isochrones (Pietrinferni et al. 2004, 2006) and assuming a Milky Way-like IMF. We then feed the pPXF with these models, finding the best-fitting linear combination of MILES SSPs in order to reproduce the observed SDSS spectra. The output of pPXF provides the weight to each SSP model with given age and metallicity. Then, SFHs are computed as the sum of the weights over the metallicities as a function of age (see Cappellari 2017 for more details).

2.4 MaNGA galaxies and their SFH

Mapping Nearby Galaxies at Apache Point Observatory (MaNGA) survey (Bundy et al. 2015) is part of the fourth generation of SDSS (Section 2.3; see also Blanton et al. 2017). The data set provides optical-IFU spectroscopy for |$\sim 10,000$| galaxies at low redshift (⁠|$0.01\lt z\lt 0.15$|⁠). We use the observationally inferred SFHs of star-forming MaNGA galaxies from Bertemes & Wuyts (in preparation), which were derived based on a full spectral fitting procedure following the stellar population synthesis (SPS) method. In more detail, the integrated MaNGA spectra were fit in the optical wavelength range of |$3700 {\small --} 7400$| Å simultaneously with associated photometry from the NASA-Sloan Atlas (Blanton et al. 2011) using the Bagpipes code (Carnall et al. 2018). Bagpipes is based on the 2016 version of the Bruzual & Charlot (2003) models, and assumes a Kroupa & Boily (2002) IMF.

In this paper, we use piecewise constant SFHs, which consist of 7 segments of constant SFR in 7 age bins, and thus leave a significant amount of freedom to the evolution of the SFR in time. To favour smooth SFHs in the absence of constraining information, the jumps in SFR between the age bins are drawn from a prior corresponding to a Student’s t-distribution (Leja et al. 2019a). A Calzetti et al. (2000) extinction law was assumed (with the dust attenuation |$A_{V}$| being a free parameter), and all stars were assumed to share a single metallicity, which is left to vary between |$0.1 \, Z_\odot$| and |$2 \, Z_\odot$|⁠. All emission lines were subtracted for the fitting process by using the emission line only (EMLINE) model cube produced by the MaNGA DAP (Westfall et al. 2019). Further, spectra were scaled to an |$\rm S/N$| of 30 prior to fitting to allow the procedure to sufficiently explore the parameter space.

2.5 GAMA

The Galaxy And Mass Assembly (GAMA) survey (Driver et al. 2011; Liske et al. 2015) is a spectroscopic survey which provides redshifts for |$\sim$|300 000 galaxies over 5 regions with a total sky area of 230 square degrees. The spectroscopic data are accompanied by 20-band photometry from the ultraviolet to the infrared, providing a vast multiwavelength data set which allows us to measure SFHs, SFRs, and stellar masses of the detected galaxies via SED fitting. As of the final data release DR4² (Driver et al. 2022), which is based on a new derivation of the underlying photometry (Bellstedt et al. 2020a) using the source-extracting software profound (Robotham et al. 2018), the survey is 95 per cent complete down to an r-band magnitude of 19.65 mag.

We use the SFHs from Bellstedt et al. (2020b), who studied a sample of 6688 galaxies with redshift |$z\lt 0.06$| and |$r\lt 19.5 \, {\rm mag}$| in the three equatorial fields. The 20-band photometry were processed with the SED fitting code prospect (Robotham et al. 2020) to derive the SFHs, using the parametric massfunc_snorm_trunc function to describe the SFH (corresponding to a skewed Normal function with a truncation in the early Universe), and an evolving metallicity implementation, where the final metallicity is fitted as a free parameter. For further details on the SED fitting implementation, we refer the reader to Bellstedt et al. (2020b).

3 THE DECODE MODELLING

To probe the consistency among the different data sets, we make use of the decode presented in Paper I. decode relies on statistical merger trees generated stochastically from analytical halo mass functions (HMF) and subhalo mass functions (SHMF). Halo accretion tracks and merger histories are converted transparently into stellar mass growths via input SMHM relations.

In this Section, we describe our recipes to model the evolution of satellite galaxies after the time of infall, such as star formation and stellar stripping, and to predict the mean SFH of central galaxies. Our aim in this work is to model the evolution of central and satellite galaxies with minimal assumptions and input parameters, to be guided as much as possible by the data in extracting constraints on the evolution of galaxies in a cosmological context.

The main steps that we follow to compute SFHs can be summarized as follows:

• generation of halo merger trees via decode (Paper I);

• assignment of galaxies to dark matter haloes (Section 3.1);

• evolution of satellite galaxies after infall (Section 3.2);

• prediction of merger histories (Paper I);

• computing SFHs (Section 3.3).

Fig. 1 shows a cartoon sketching the idea behind decode, giving a basic description of how different shapes and evolution of the input SMF (therefore SMHM relation) can impact the output SFHs. In brief, our methodology is based on the dark matter merger trees generated following the recipe presented in Paper I, and converted into galaxy mergers to compute the mean cumulative mass assembly history of mergers. The total stellar mass growth of a galaxy is computed by converting the total halo mass assembly via the SMHM relation. Then, the stellar mass growth history from star formation is derived by computing the difference between the total mass growth and the merger history, as we will further detail in Section 3.3. In this work, we include the evolution of satellites after their infall, unlike in Paper I where we focused only on a frozen model where the mass of satellites is assumed to be constant. While a frozen model represents a good first-order approximation, it is not realistic to ignore any further evolutionary process affecting individual satellites after infall, although we will see that in some instances, a frozen model provides a very good approximation, as in the comparison with the TNG100 simulation (see Section 3.4). In this new version of decode, we therefore also allow for the possibility for the satellites to form stars, quench their star formation and to undergo stellar tidal stripping. The mass stripped from satellites or the mass lost during mergers will contribute to the formation of the ICL. Depending on the evolution of the merging and surviving satellites, the merger and SFH of central galaxies will also change consequently. Throughout, we assume that |$20{{\ \rm per\ cent}}$| of the stellar mass is lost during each merger, following the results of several SAMs and hydrodynamical simulations (e.g. Contini et al. 2014; Fontanot et al. 2015; Kannan et al. 2015; Koutsouridou & Cattaneo 2022). We note that by increasing such mass-loss fraction up to |$\sim 40{{\ \rm per\ cent}}$|⁠, which is among the highest values adopted in cosmological models in the literature (e.g. Moster et al. 2018), the resulting galaxy merger rate is reduced by a factor of 0.1 dex, even for the most massive galaxies with stellar mass |$M_\star \sim 10^{12}\, \mathrm{M}_\odot$|⁠, but the rest of the relevant observables is barely impacted, including the satellite mass function and the SFHs.

Ideally, SEMs, developed on top of a dark matter-dominated hierarchical Universe, are designed to be self-consistent between their input and output observables. This means, in our case, that the input observed SMF should produce an SMHM relation which predicts SFHs, morphologies, ICL, satellite abundances and so on, self-consistent among each other when compared simultaneously with observational data. Therefore, simultaneously matching robust and homogeneous data sets via semi-empirical hierarchical models like ours, will not only be able to test the input assumptions on, e.g. infall time-scales and amount of stellar stripping, but also will probe the degree of reliability of a |$\Lambda$|CDM-based galaxy evolution model in reproducing the real Universe. However, these are challenging goals for a number of reasons. First of all, observations are not fully consistent among themselves, for instance integrated SFHs have often resulted in the literature to overproduce the mass locked up in galaxies when compared to the integrated SMFs at any given epoch, as further discussed below. Furthermore, several observables predicted by decode depend not only on the input SMF, but also on the way we model satellites evolution, which can be affected by a number of, sometime degenerate, physical processes (e.g. star formation, quenching, and stripping). Nevertheless, we will discuss below that some observables are more dependent than others on post-infall satellite evolution and it is thus possible to pin down the effectiveness of a given SMF and SMHM relation to reproduce different data sets.

As shown in Fig. 1, the working flow of decode, from the left to the right, starts by applying abundance matching between the SMF and central halo mass function to define a redshift and mass-dependent SMHM relation. The SMHM relation is in turn used to translate halo assembly histories into galaxy growth histories, which are then decomposed into mergers and SFHs as detailed below.

3.1 Abundance matching

The connection between the galaxy stellar mass and the host dark matter halo mass is one of the main ingredients of decode. Our SMHM relation is computed via the abundance matching technique presented in Paper I (equation 8) according to the formalism of Aversa et al. (2015)

where |$\tilde {\sigma }_{\log M_*} = \sigma _{\log M_*} / \mu$|⁠, with |$\sigma _{\log M_*}$| being the Gaussian scatter in stellar mass at given halo mass and |$\mu = \mathrm{d} \log M_* / \mathrm{d} \log M_h$| the derivative of the SMHM relation. This method allows to compute the mean stellar mass at fixed halo mass taking as input the dark matter HMF and galaxy SMF, and an input scatter in stellar mass. We make use of the total HMF which includes also surviving unstripped subhaloes, which is computed via the correction to the parent HMF as described in Appendix B in Paper I. For the SMF, we employ several models that are described in Section 4.1. We then make use of the output SMHM relation to populate parent haloes and infalling unstripped subhaloes with galaxies. We remind the reader that although decode generates a statistical ensemble of merger trees, it only robustly predicts mean galaxy assembly, merger, and SFHs, as detailed in Appendix C of Paper I. We checked that all the sanity checks carried out in Paper I still hold even after inclusion of different satellite evolutionary histories.

3.2 Satellite evolution

In Paper I, we assumed that the mass of satellite galaxies remains mostly constant after infall. This assumption allows to rapidly predict the galaxy merger rates, satellite abundances and star formation rates, with minimal computational time. Although the mass of each individual satellite is expected to evolve after infall via several physical processes such as stripping, star formation, and quenching (e.g. Cattaneo et al. 2011; Wetzel et al. 2013; Fillingham et al. 2016; Smith et al. 2016; Shi et al. 2020; Wright et al. 2022; Engler et al. 2023), we will show in Section 4.3 that a frozen model is sufficiently accurate for our purposes. In this work, we show how these processes impact the prediction of the satellite abundances, merger rates, and SFH of central galaxies. Below, we provide a brief description of how we model these additional key physical processes.

3.2.1 Stellar stripping

We apply the recipe for stellar stripping according to the results from the N-body simulations by Smith et al. (2016)

where |$f_{\rm DM}$| is the ratio between the mass of the subhalo at given redshift z and its peak mass, which we assume to be in good approximation equal to the mass of the subhalo at infall, and |$f_{\rm str}$| is the fraction of the satellite’s stellar mass at redshift z after stripping.

3.2.2 Satellites star formation and quenching

For each model rendition we first generate the mean main-sequence SFR–stellar mass relation (see Section 4.2), which we then use to assign an SFR to the infalling satellites. We assume that the satellites have enough gas reservoirs to sustain this SFR until a given quenching time that we define below. We stress that in decode, the SFR for centrals is a prediction, while for the satellites it is an input assumption. Obviously, the prediction of the former slightly depends on the latter according to the prescription of Section 3.3, since the SFHs of centrals depend on the merger histories, and therefore on the way we treat satellites evolution. We also allow for the possibility to include a scatter in the SFR–|$M_\star$| relation as well. However, we found that by adding a scatter up to 0.2 dex, the predicted mean quantities of interest here do not vary appreciably.

We quench the star formation of satellites following the recipe from Wetzel et al. (2013) of the ‘delayed-then-rapid’ quenching, according to which satellite galaxies continue to form stars for a period equal to the delay time-scale |$\tau _{\rm delay}$| (see Section 4.3 therein), and after which they undergo a rapid quenching where the SFR is truncated exponentially.

3.3 Star formation histories of centrals

Star formation is an additional process, complementary to mergers, significantly contributing to the mass assembly history of a galaxy. In order to compute the SFH of central galaxies, we simply assume that the main components of the stellar mass growth of a galaxy are in first approximation (1) the stellar mass coming from all mergers with other satellite galaxies and (2) the in situ formation of stars. The stellar mass formed in situ is derived by the difference between the total stellar mass growth history, which is inferred by converting the halo mass evolutionary track to stellar mass via the SMHM relation, and that from the cumulative merger history of the galaxy

|$M_*^{\rm SF}$| indicates the stellar mass formed via star formation including also (1) any star formation triggered by wet mergers (e.g. Hopkins et al. 2008), (2) any star formation from secular process. Our SFRs, when integrated to yield the fractional stellar mass originating from star formation, have to then be corrected at each time step by the cumulative amount of stellar mass returned in gas to the interstellar medium (ISM), described by the global gas mass loss rate (GMLR). In this paper, we make use of the recipe presented in section 3.2 of Leitner & Kravtsov (2011) (see also discussion therein) to account for the gas loss component to derive the SFRs from stellar masses, based on the equation

where |$\dot{M}_\star ^{\rm SF} = \mathrm{d} M_\star ^{\rm SF} / \mathrm{d} t$| is the uncorrected SFR, and |$f_{\rm ML}$| is the fraction of mass loss given by equation (1) of Leitner & Kravtsov (2011).

3.4 Validating DECODE’s self-consistency

Before progressing into generating all decode’s predictions for the mass assembly histories of galaxies, SFHs, merger rates, ICL fractions, using observed SMFs as input, we proceed by testing the inner self-consistency of decode by comparing with the TNG100 simulation. We specifically aim to test if, by starting from the SMF generated by the TNG100 simulation, decode is able to predict mass assembly histories, SFHs, and merger rates consistent with those in the TNG100 simulation. We stress that the aim here is not to use the TNG data to calibrate our modelling, but to test the self-consistency and performance of decode in predicting the aforementioned quantities by taking TNG’s inputs, i.e. SMHM relation and scatter in stellar mass.

First of all, we test our abundance matching prescription described in Section 3.1, which takes as inputs the galaxy SMF and the HMF extracted from the TNG simulation. We find that the SMHM relation implied by equation (1) is in good agreement with the one directly computed from the simulation itself at all redshifts, demonstrating the robustness of our abundance matching procedure in determining the right mapping between central galaxies and host haloes. We show the comparison in Fig. A1 in Appendix  A for redshifts |$z=0, 1 \, {\rm and} \, 2$|⁠, where the upper panels show the SMHM relation from our abundance matching compared to those from the TNG, and the lower panels show the logarithmic difference between the two of them, with a small residual of less than 0.1 dex.

Secondly, we test the self-consistency of our methodology in calculating the SFH of galaxies with stellar mass today of |$M_\star (z=0) \lesssim 10^{11} \, \mathrm{M}_\odot$|⁠. In Fig. 2, we show the mean stellar mass accreted in situ (blue lines) and ex situ (orange lines) for four different mass bins of central galaxies at redshift |$z=0$|⁠, as labelled. The SFHs and merger histories predicted from decode using TNG’s SMHM relation in input are in good agreement with those extracted from the simulation itself. This test proves the validity of our methodology which we aim to extend to observationally driven SMHM relations in Section 4.2.

We now turn to test the applicability of our methodology to self-consistently predict the SFHs, or at least the mass assembly histories, of more massive galaxies. For galaxies of stellar mass above |$M_\star (z=0) \gtrsim 10^{11} \, \mathrm{M}_\odot$|⁠, the contribution from mergers becomes progressively more important, and it should be carefully accounted for when predicting SFHs. Since the total mass assembly history of galaxies in our modelling are given by the SMHM relation and the assembly histories of the host dark matter haloes, it is vital to carefully determine the contribution of mergers to derive reliable estimates of the implied SFHs (as described in Section 3.3). In order to accurately predict galaxy merger histories, it is important to have a good understanding of the evolution of merging satellite galaxies, whose mass increase or decrease will lead to different merger rates, as already discussed in Section 3. The upper panel of Fig. 3 shows the fractional mass evolution of satellite galaxies from the TNG100 simulation. The curves show the fraction of stellar mass growth of the TNG satellites as a function of their time since infall, for satellite galaxies of different masses selected at redshift |$z=0$|⁠. The aim of the plot is to highlight the stellar mass evolution of the TNG satellites since infall, regardless of their different values of infall redshift. For stellar mass above |$M_\star \gtrsim 10^{10} \, \mathrm{M}_\odot$|⁠, some galaxies have growing mass due to star formation, some others have decreasing mass due to ram-pressure stripping. However, regardless of the physical mechanism that controls the evolution of the single satellite galaxy, we find that in the TNG most of the satellites have, on average, a roughly constant stellar mass growth history since their epoch of infall. Therefore, over a large sample, the frozen model is a good approximation for the evolution of satellites in the TNG simulation. Only the very least massive satellites continue to form stars after infall but are then rapidly quenched (e.g. Shi et al. 2020; Ding et al. 2024), resulting in an average increase in their final stellar mass of roughly |$0.1{\small --}0.2$| dex since infall.

The lower panel of Fig. 3 reports the satellite SMF in the TNG simulation at |$z=0$| (dashed–dotted, purple line), compared to the predictions from decode when adopting the TNG’s SMHM relation in input, as well as the logarithmic difference between the two of them. We found that a frozen satellites model can relatively well reproduce the simulation’s satellite abundances, except for the less massive satellites which continue to form some mass via star formation after their time of infall. When including both stellar stripping and satellite star formation after infall, decode is able to reproduce the stellar mass function of satellite galaxies with high accuracy across almost the full stellar mass range considered here, further validating decode’s methodology in successfully, flexibly, and self-consistently reproducing the outputs of a comprehensive galaxy evolution model like the TNG, by only using its SMHM relation in input. It is interesting to note from the bottom panel of Fig. 3 that the detailed evolutionary histories of satellites after infall, within the modelling explored in this work, have a relatively minor effect on the resulting SMF at |$z=0$| of less than |$\sim 0.1$| dex, for satellites with stellar mass |$M_\star \gtrsim 10^{9} \, \mathrm{M}_\odot$|⁠, which is comparable or even less than the statistical uncertainties in observational estimates (e.g. Grylls et al. 2019). The only appreciable difference would apply to lower mass satellites, which could be matched by allowing a fraction of galaxies to increase their stellar mass via star formation before their cold gas being removed via tidal stripping or ram-pressure stripping (see e.g. Rohr et al. 2023; Ding et al. 2024).

Having tested the accuracy of decode in reproducing the observed local SMF of satellite galaxies, we are now in a position to proceed to the comparison between predicted and observed stellar mass assembly histories of massive, central galaxies. The upper panels of Fig. 4 show the cumulative stellar mass accreted from mergers for two mass bins of central galaxy (⁠|$M_\star (z=0) = 10^{11.5} \, {\rm and} \, 10^{12} \, \mathrm{M}_\odot$|⁠), as labelled in the figure, computed via decode using TNG’s SMHM relation as input compared to the merger histories computed from the simulation directly. We found that decode is capable of reproducing fairly well (within |$\sim 0.1$| dex) the merger histories of the TNG simulation when taking as input its SMHM relation and adopting a frozen model for the mass of satellites. In both stellar mass bins the scenario where satellites lose their mass through stripping tends to slightly underestimate the TNG’s merger histories. On the other hand, a model, in which satellites evolve either via star formation only or both star formation and stellar stripping, tends to somewhat overestimate the TNG’s merger histories. At the level in which it is modelled in decode, the stellar stripping tends to have less impact than the SFR in shaping satellites after infall, and thus less influence on the merger histories.

The lower panels of Fig. 4 show the SFHs for the same galaxy stellar mass bins. As expected, at fixed stellar mass bin, a satellite evolutionary model that produces a higher merger history predicts lower SFHs and vice versa. Similarly as for the merger histories, the frozen satellites model allows to reproduce the TNG’s SFHs fairly in agreement. We note that, for the stellar mass bin of |$M_\star (z=0) = 10^{12} \, \mathrm{M}_\odot$|⁠, overall decode is slightly less performant in reproducing TNG’s star formation history for a couple of reasons. First of all, the limited number of clusters (⁠|$\lt 10$|⁠) in the TNG100 simulation could introduce some bias in the sample. More crucially, the bright end of the SMHM relation computed from the simulation could be affected by a non-negligible numerical uncertainty, with tiny variations in the slope causing large discrepancies in the mass assembly histories in decode.

4 RESULTS

Having tested the efficacy of decode in self-consistently reproducing the mean galaxy assembly histories of the TNG simulation in terms of star formation, mergers, and satellite abundances, by using in input only their SMHM relation, we are now ready to turn to real data. The aim of this section is to use the latest multiple (but different) determinations of the SMFs at both low and high redshifts, extract via abundance matching the respective SMHM relations implied by these SMFs in a |$\Lambda$|CDM cosmological context, and then derive the expected SFHs, mergers, and satellite abundances against the corresponding, independently derived data sets. Our aim here is two-fold, we aim (1) to test whether the currently available data sets are inherently self-consistent with each other and (2) to discern the best models capable of simultaneously reproducing most of the data sets considered here, including ICL and satellite abundances. This section is structured as follows. We will first present the determinations of the expected SMHM derived from current data on the SMFs in Section 4.1. We will then show the implied SFHs for low-mass and high-mass central galaxies in Section 4.2, the corresponding populations of surviving satellites and the implied ICL in Sections 4.3 and 4.4, respectively.

4.1 Stellar mass-halo mass models

As already discussed by many works in the literature (e.g. Shankar et al. 2006 ; Moster et al. 2010, 2013; Guo et al. 2011 ; Grylls et al. 2020, Paper I), the SMHM relation is directly determined by the input observed SMF, the shape and evolution of which are still a matter of hot debate, even at low redshifts. On one hand, some groups suggested a negligible time evolution of the SMF. For example, Moustakas et al. (2013), Bernardi et al. (2013, 2016, 2017), and Leja et al. (2020) showed that there is no evolution in the high-mass end of the SMF in PRIMUS, SDSS, CMASS, and 3D-HST up to |$z \simeq 1$|⁠. Kawinwanichakij et al. (2020) suggested that there is no evolution in the SHELA survey’s SMF even up to |$z=1.5$|⁠. On the other hand, other groups found more substantial evolution (e.g. Tomczak et al. 2014; Davidzon et al. 2017; Leja et al. 2020; Weaver et al. 2023). At any given redshift, even at |$z=0$|⁠, the shape of the SMF is still far from being firmly established. Bernardi et al. (2010, 2013, 2016) extensively discussed how different choices of background, photometry, mass-to-light ratios, have a profound impact on the SMF, especially at high stellar masses. Other groups have derived very different shapes at both low and high masses and even around the knee of the SMF at all measured redshifts (e.g. Tomczak et al. 2014; Davidzon et al. 2017; Kawinwanichakij et al. 2020; Weibel et al. 2024).

In Paper I, we found that an SMHM relation implied by an SMF with flatter and evolving bright end is more suitable to describe the galaxy major merger rates, elliptical fractions and bulge-to-total distributions, compared to a model derived from a constant SMF up to |$z=1.5$| at least for |$M_\star \gtrsim 10^{11} \, \mathrm{M}_\odot$|⁠. On the other hand, the constraints on the shape of the input SMHM relation become gradually significantly less tight at stellar masses below |$10^{11} \, \mathrm{M}_\odot$|⁠, for example, major mergers are relatively rare at lower stellar masses (e.g. Hopkins et al. 2010a, b; Rodriguez-Gomez et al. 2015; O’Leary et al. 2021; Huško, Lacey & Baugh 2022, Paper I). In what follows, we will underline the role of SFHs and their close link with the shape of the SMF in particular at low masses. Flatter SFHs are on average better mapped into slowly evolving SMFs, as expected also from basic continuity equation arguments (e.g. Leja et al. 2015; Grylls et al. 2020). To highlight the connections between SFHs and SMFs, we will consider two models for the SMHM relation that broadly bracket the current measurements of the SMF at different redshifts:

• Constant SMF model: SMF constant up to redshift |$z \sim 1.5$| (Bernardi et al. 2017) and gradually dropping in normalization above this redshift according to equation (11) of Paper I, |$\log _{10} \phi (M_\star (z)) \simeq f(z) \cdot \log _{10} \phi (M_\star (z=0.1))$|⁠, where |$f(z)=1$| for |$z\le 1.5$| and |$f(z) = (0.99 + 0.13(z-1.5))$| for |$z\gt 1.5$|⁠.

• Evolving SMF model: SMF characterized by a weakly evolving low-mass end (below the knee) at low redshifts |$z \lesssim 1.5$| following Davidzon et al. (2017), as plotted in the top panel of Fig. 5.

We note that Davidzon et al. (2017)’s SMF predicts less abundant number of massive galaxies with respect to Bernardi et al. (2017), which could be attributed to the choice of different mass-to-light ratios and/or different photometries. As we are here mainly interested in the differences in the evolution of the SMFs, and their impact in the SFHs, we have added a scatter of 0.3 dex to the original fit of Davidzon et al. (2017)’s SMF to align their fit to Bernardi et al. (2017)’s, but retained their original redshift evolution. We also note that Davidzon et al. (2017)’s double power-law fit to the SMF predicts slightly less abundant numbers of galaxies around the knee with respect to Bernardi et al. (2017)’s estimate for the SMF (solid blue line and data in the upper panel of Fig. 5), but this relatively small difference has little impact on any of the results presented below. Furthermore, we notice that the SMF built in this way is well consistent in shape and normalization at high redshifts (⁠|$z\sim 4{\small --}5$|⁠) with the latest inferred JWST SMF from Weibel et al. (2024), who reported a steep low-mass end of the SMF at those redshifts with moderate evolution.

The upper panel of Fig. 5 shows the two models of SMF described above and the lower panel shows the redshift evolution of the SMHM relation computed using the evolving SMF model in the redshift range |$0 \lt z \lt 3$|⁠. The resulting SMHM relation from the evolving model is characterized by an increase in stellar mass at fixed halo mass of 0.29 dex from |$z=0$| to |$z=1$| in the high-mass range (⁠|$M_{\rm h} \sim 10^{14}\, \mathrm{M}_\odot$|⁠), and by a decrease of 0.24 dex in the low-mass range (⁠|$M_{\rm h} \sim 10^{11.5}\, \mathrm{M}_\odot$|⁠). For both models, we choose an input scatter in stellar mass at fixed halo mass of 0.15 dex, noticing that varying this input parameter, within reason, does not affect any of our main conclusions. In Paper I, indeed, we showed that a redshift- and/or mass-dependent scatter does not significantly alter the merger rates and satellite abundances, mainly because varying the scatter mainly impacts the steepening of the high-mass end of SMHM relation, and thus all the observables below |$10^{11}\, \mathrm{M}_\odot$|⁠, which is the main focus of this work, remain largely unaltered.

4.2 Star formation histories

As already discussed above, the high-mass end slope of the SMHM relation has a profound imprint on the major merger histories of galaxies, steeper relations would preferentially correspond to less major mergers (e.g. Hopkins et al. 2010a; Grylls et al. 2020; O’Leary et al. 2021, Paper I). On the other hand, at fixed total stellar mass, a larger contribution of (major) mergers will inevitably correspond to less mass growth via star formation, and thus the shape of the input SMHM relation is expected to also have a significant impact on the SFH of a galaxy.

In Fig. 6, we compare decode’s predictions on the SFHs of galaxies of different stellar masses at |$z=0$|⁠, with the SFHs extracted from SED fitting from different surveys and other theoretical models, as labelled. We here focus only on galaxies with stellar mass |$M_\star \lesssim 10^{11}\, \mathrm{M}_\odot$| at |$z=0$|⁠. We restrict to this mass limit in this figure because we want to isolate the full effect of SMF variations on the SFHs. Most theoretical models indeed support the view that galaxy growth is fully SF dominated below |$10^{11}\, \mathrm{M}_\odot$| (e.g. Oser et al. 2010; Lackner et al. 2012; Rodriguez-Gomez et al. 2016; Pillepich et al. 2018b; Davison et al. 2020; Cannarozzo et al. 2023; Trujillo-Gomez et al. 2023), and decode indeed predicts the same behaviour (see also Paper I). More specifically, we compare decode’s outputs (red dashed–dotted and brown solid lines) with the TNG100 hydrodynamical simulation (green-dotted lines), GalICS SAM (blue dotted lines), and the inferred (integrated) SFHs from the MaNGA and GAMA surveys (orange and cyan solid lines and shaded areas). There are some differences between the MaNGA and GAMA SFHs, especially at stellar masses |$M_\star \gt 10^{10} \, \mathrm{M}_\odot$|⁠, which are possibly due to the differences between the two data sets themselves. First of all, the two surveys cover different wavelengths, with MaNGA focusing only in the optical regime, while GAMA covers a broad wavelength range from the UV to the IR. Furthermore, due to the smooth truncation of the skewed normal function shape at the beginning of the Universe, the GAMA SFHs might have higher SFRs at high redshifts with respect to MaNGA, where no truncation is used. Additionally, the usage of constant metallicity without modelling the chemical evolution in the MaNGA fit might introduce some slight bias towards younger galaxies compared to those from GAMA. Additionally, we note that there might be a tiny shift in redshift between the SFHs from theoretical models and observations, due to the different mean redshifts of the samples from the latter, being |$z\lesssim 0.1$| for SDSS and MaNGA, and |$z\lt 0.06$| for GAMA. We have checked that this effect would be negligible compared to the effects introduced by the systematics of the surveys.

The mean SFHs in decode are computed by selecting in the stellar bin of interest a relatively large number of galaxies, which are then evolved in mass following the recipes described in Section 3.3. Here, we show the SFHs from decode computed via the frozen stellar mass model for satellites. Indeed, similarly to the tests on the TNG in Section 3.4, we have checked that the SFHs and SFRs for central galaxies of |$M_\star \lesssim 10^{11}\, \mathrm{M}_\odot$| are not altered by the satellites evolution model, varying by less than 0.05 dex (see Appendix  B).

The constant input SMF model generates SFHs that are broadly consistent with the SFHs from the TNG simulation, GalICS and MaNGA at |$M_\star (z=0) = 10^{9.5} \, {\rm and} \, 10^{10} \, \mathrm{M}_\odot$|⁠, while it produces flatter SFHs for the other two mass bins, being more consistent with the GAMA and SDSS surveys. On the other hand, the SMHM relation implied by an evolving SMF produces mean SFHs steeper than those predicted by a constant SMF model, as expected, given that there is more stellar mass build-up in the former model. The evolving SMF model produces SFHs that are steeper than all the data at |$M_\star \lesssim 10^{10}\, \mathrm{M}_\odot$| (top panels), while at larger masses |$M_\star \sim 10^{10} {\small --} 10^{11} \, \mathrm{M}_\odot$| they are consistent with what inferred from MaNGA, but significantly steeper than those from GAMA and SDSS (bottom panels). The SFHs from the evolving SMF model are broadly consistent with those predicted by the TNG100 at |$M_\star \sim 10^{10} {\small --} 10^{11} \, \mathrm{M}_\odot$|⁠, a trend expected as the SMF of the TNG100 has a weaker evolution below the knee but predicting a stronger evolution at higher stellar masses (see e.g. Pillepich et al. 2018b).

While at lower masses, the total stellar mass growth of a galaxy is in good approximation mostly contributed by star formation alone, for more massive galaxies the comparison with observations is less straightforward as mergers can provide a non-negligible contribution to the overall mass assembly histories of galaxies, especially for SMHM relations characterized by a flatter high-mass slope (e.g. Grylls et al. 2020, Paper I), as discussed above. On the other hand, the SFHs retrieved from SED fitting record the mass formed in stars but cannot distinguish between the mass formed ex situ or in situ. For this reason, in what follows we do not compare directly with measured SFHs, but rather show the predicted mass grow histories of massive galaxies predicted by decode using our constant and evolving SMF models, distinguishing between the mass growth via in situ star formation and via mergers. In what follows we will mainly present the results using the evolving SMF model, but highlight the differences in our results when switching to a constant SMF model. In addition, for completeness, we also discuss the predicted mass accretion histories for different types of satellite evolution after infall.

Fig. 7 shows the mean merger histories (upper panels) and the integrated SFHs (central panels) for galaxies with stellar masses of |$M_\star = 10^{11.5} \,\, {\rm and} \,\, 10^{12} \, \mathrm{M}_\odot$| at |$z=0$|⁠, as well as the specific SFRs as a function of time (lower panels). The results in the plot are shown for the evolving SMF as input, and we checked that the stellar mass growths change by less than |$\sim 0.15$| dex in normalization but maintain the same redshift evolution when using the constant SMF as input. We find that, all the merger histories predicted by decode are similar to each other, irrespective of the chosen model for the satellite evolution (orange dotted, blue solid, green long dashed lines, as labelled), with only a slight marginal increase of |$\sim 10{{\ \rm per\ cent}}$| in stellar mass growth when satellites continue forming stars after infall in the most massive centrals (orange dotted lines, top right panel). In this respect, the shape of the SMHM relation has a larger impact than the satellite evolution in shaping the merger histories of central galaxies.

The predicted stellar mass growth via in situ star formation (central panels) is relatively smaller than the one from mergers by a factor of |$\sim 2$| and |$\sim 3$| in the two stellar mass bins, respectively, for both the cases of evolving and constant SMF. It is interesting to note that the SFHs, being overall contributing less to the stellar mass growth in these massive galaxies, are proportionally more dependent on the satellite evolution, in particular the model with star-forming satellites tends to predict an SFH lower by a factor of |$\sim 2$| than that other two models (orange dotted versus blue solid and green long dashed lines). We stress that the total stellar mass growth remains the same in all models at fixed SMHM relation, i.e. at fixed SMF model, but only the relative proportions of merging and star formation change depending on the satellite evolution. Moreover, the stronger dependence of the in situ stellar mass growth on the satellite evolution with respect to the ex situ growth is an apparent visual effect, due to the high ex situ fraction with respect to the in situ fraction across cosmic time. Indeed, differences of the order of |${\rm few} \times 10^{-2} \, {\rm dex}$| in the ex situ growth will translate into differences of 10 times larger in the in situ growth, causing the apparent more sensible dependence on the satellite evolution. The specific star formation rates (sSFR) plotted in the bottom panels of Fig. 7 all show, irrespective of the evolution of satellites, a clear decreasing trend with cosmic time, steepening in more massive galaxies, mimicking the overall gas starvation that all galaxies are expected to undergo given the strong fall off in the cosmic SFR below |$z\sim 1{\small --}2$| (e.g. Madau & Dickinson 2014). It is interesting to note that in our data-driven model decode, which does not contain any recipe for gas exhaustion or quenching, it is ultimately the shape of the halo accretion rate that drives the stellar mass growth and the sSFR histories, again lending further support to the capital importance of the underlying halo assembly histories in driving the growth of central galaxies (e.g. Neistein & Weinmann 2010; Wechsler & Tinker 2018; Bose et al. 2019; Jiang et al. 2021; Boco et al. 2023; Lyu et al. 2023).

Having discussed the mean SFHs of galaxies of different mass along their evolutionary tracks, we now move to the study of the overall SFR–stellar mass relation or main sequence of galaxies as predicted by decode using our two reference SMF models with and without redshift evolution. In what follows, we focus on the lower mass range, |$M_\star \lesssim 10^{11.5} \, \mathrm{M}_\odot$|⁠, where the fraction of star-forming galaxies dominates the stellar population (e.g. Tomczak et al. 2014; Moutard et al. 2016; Davidzon et al. 2017; Moster et al. 2018; Pillepich et al. 2018b; Behroozi et al. 2019; Donnari et al. 2021; Weaver et al. 2023), guaranteeing a fairly consistent comparison with the observational estimates of the main sequence. Indeed, as also stressed in Paper I, in its present form, decode generates stochastic realizations of galaxies without distinction between star-forming or quenched galaxies. Nevertheless, we also checked that, weighting the mean SFR–|$M_\star$| relation predicted by decode only by the relative fraction of star-forming galaxies inferred by, e.g. Weaver et al. (2023), the resulting main-sequence relation is changed by less than 0.1 dex at high redshifts and up to 0.2 dex at |$z=0$| at |$M_\star \gtrsim 10^{11} \, \mathrm{M}_\odot$|⁠, with no impact on its shape, normalization and evolution.

The top panel of Fig. 8 shows the main sequence as predicted by the evolving SMF model at different redshifts, as labelled. The constant SMF would predict very similar shapes but lower in normalization by a factor of |$\sim 1.5{\small --}2$| at any given redshift up to |$z\sim 3$|⁠. We show two realizations of the model, the complete model (solid lines) and the model in which we assume galaxies only grow via star formation at all redshifts and masses (dashed lines), in other words implying that the infalling satellites always have very long dynamical time-scales. The latter rendition of the model is clearly possibly an oversimplification, and it would predict too many satellites in the local Universe (see Section 4.3), but it is included to provide an overview of the shape of the main sequence in a ‘maximal’ model dominated by star formation. We first note from the top panel of Fig. 8 that the SFR increases at the same pace in redshift at fixed stellar masses, especially below |$M_\star \lesssim 3\times 10^{10}\, \mathrm{M}_\odot$|⁠. At larger stellar masses, the increasing rate of the SFR is actually larger, especially when effective mergers are included, which are the main channel that grows the stellar mass in early-type galaxies.

It is also interesting to note that there is no sign of a downturn in the predicted main sequence at |$z\gt 1$| above |$M_\star \gtrsim 3\times 10^{10}\, \mathrm{M}_\odot$|⁠. A break or flattening in the main sequence starts appearing only at |$z\lt 1$| (broadly consistent with observation at |$z\sim 1{\small --}2$|⁠), and this break is visible even in the extreme model with only SFR. The flattening of the main sequence is a mere consequence of the double power law shape in the input SMHM relation, which implies less stellar mass growth when the central galaxy crosses the knee of the SMHM relation. In the full model inclusive of mergers, the flattening becomes even more evident at |$z\lt 1$|⁠, creating a downturn in the main sequence, because the majority of massive galaxies is predicted to be quenched with their mass growth being dominated by mergers. In this sense, SEMs represent a powerful tool to probe the effect of the evolution of the SMF, or SMHM relation, on the galaxy SFRs. From the observational point of view, the existence of the break has been widely discussed in the literature with contrasting results, with several works suggesting a single power law shape at all stellar masses and redshifts (e.g. Rodighiero et al. 2014; Speagle et al. 2014; Pearson et al. 2018), and others suggesting a break towards high stellar masses (e.g. Whitaker et al. 2014; Lee et al. 2015; Tomczak et al. 2016; Thorne et al. 2021; Leja et al. 2022; Popesso et al.2023).

The middle and lower panels of Fig. 8 show a close comparison between the main sequence predicted by the SMF evolving model from decode and three reference observational results from Leja et al. (2022), Popesso et al. (2023), and Koprowski et al. (2024), as well as, for completeness, the predictions from the TNG100 simulation. First off, it is important to note that this comparison is only at a qualitative level as the definitions of stellar masses, despite a tiny deviation due to the choice of the IMF,³ may still be systematically different in Popesso et al. (2023), Bernardi et al. (2017), and Davidzon et al. (2017), which are the ones on which our reference models are based on. At face value, decode tends to generate a mean SFR that is a factor of |$\sim 2{\small --}3$| systematically lower than in the data. Leja et al. (2022)’s results are based on neural networks to parametrize the galaxy population density in the main-sequence plane. With this method, Leja et al. (2022) find a mean SFR which is also lower than other direct estimates, in better agreement with our models, at least at |$z\sim 1$|⁠. Also Koprowski et al. (2024) and the TNG100 tends to predict on average lower SFRs than Popesso et al. (2023) and in better agreement with our predictions. Leja et al. (2022) also showed that their sSFR time evolution for galaxies of |$M_\star = 10^{10}\, \mathrm{M}_\odot$| well matches those from numerical simulations, such as EAGLE and IllustrisTNG, even though the main sequence still slightly differs at high redshifts. It has been discussed several times in the literature that some estimates of the SFR–|$M_\star$| relations may overproduce the stellar mass density recorded at any given epoch (e.g. Bernardi et al. 2010; Rodríguez-Puebla et al. 2017; Donnari et al. 2019; Hashemizadeh et al. 2021; Leja et al. 2022), which is also what we infer here using the evolving SMF model. On the other hand, other groups have found better agreement between the integrated SFR and SMFs (e.g. Bellstedt et al. 2020b). Our results further support the intimate link between observed SFRs and measured SMFs, and how each of these observations can provide valuable constraints on the other one.

4.3 Satellite abundances

As already discussed in other works (e.g. Moster et al. 2018; Behroozi et al. 2019; Grylls et al. 2019, Paper I), the abundances of surviving satellite galaxies provide a complementary and independent test on the shape of the input SMHM relation. Indeed, satellite abundances can be considered as the other side of the coin with respect to mergers, in a dark matter-dominated hierarchical Universe they are effectively ‘failed’ mergers at any given epoch, and can thus be used as independent test of the same merger model, and compared against direct observational estimates of the satellite abundances. Any successful cosmological model of galaxy evolution that matches the SFHs and merger histories of galaxies, must also provide the right amount of observed satellites. Any failure in this respect could be, as for any other observable considered so far, attributed to either a shortcoming of the model and/or inconsistencies in the different and independent data sets. Before presenting decode’s predictions on the satellite abundances, we stress once again that the satellite stellar mass function is a genuine prediction of the model, as it depends on both the input SMHM relation, which initialises the subhaloes at infall, and on the prescription for the evolution of the satellites after infall.

The top panel of Fig. 9 compares decode’s predicted satellite SMF (solid blue line) in the evolving SMF model with frozen satellites against the satellite SDSS galaxy data from Bernardi et al. (2017) as labelled in Yang et al. (2007)’s halo catalogue (blue filled circles). The results with the constant SMF as input do not change appreciably, since the satellites SMF is not very sensible to the input SMHM relation, as already shown in Paper I. It is apparent that the model, without any further fine-tuning, nicely matches the data, at least above |$M_{\rm \star ,sat} \gtrsim 3\times 10^{10}\, \mathrm{M}_\odot$|⁠, while there is a non-negligible shortfall of |$\sim 0.1$| dex at lower masses. When dividing the observed SDSS sample in satellites hosted in parent haloes above and below a chosen host halo mass of |$M_{\rm h,par} = 10^{13}\, \mathrm{M}_\odot$| (orange triangles and green squares), we see that the major shortfall occurs in lower mass haloes. Allowing for a longer dynamical friction time-scale (⁠|$f_{\rm dyn}=1$|⁠) with respect to decode’s reference time-scale (equation 9 of Paper I), improves the fit at low stellar masses both in lower mass parent haloes as well as in the total satellite SMF (dotted lines). A similar, nearly degenerate, effect of boosting in the number density of satellites is also found by allowing the satellites to continue forming stars after infall as described in Section 3.2.2 (dashed lines). We further show in the bottom panel, the effect on the total satellite mass function of including late, post-infall star formation (orange dashed line) and star formation and stellar stripping (red, dash–dotted line), which provides an improved fit to the data with respect to a frozen model at low stellar masses below |$M_{\rm \star ,sat} \lesssim 3\times 10^{10}\, \mathrm{M}_\odot$|⁠.

Although the two model renditions discussed above characterized by a lower dynamical friction time-scale or by the inclusion of star formation and stripping, provide both a nearly degenerate match to the local stellar mass function of satellites, we note that some stellar stripping in satellites is expected, and indeed observed, in groups and clusters (e.g. Poggianti et al. 2017; Franchetto et al. 2021; Akerman et al. 2024), and can also contribute to the ICL level measured in clusters, as discussed in the next Section.

4.4 Intracluster light

As a further application of decode and additional constraint for the assembly history of galaxies, we study the ICL, i.e. the faint diffuse light coming from stars which freely floats within each cluster’s gravitational potential and is not bound to any galaxy within the cluster (e.g. Montes 2022). Many observational works suggested that the origin of the ICL may be due to tidal stripping of satellites and galaxy–galaxy mergers (e.g. Gregg & West 1998; Mihos et al. 2005; Contini et al. 2014; Contini, Yi & Kang 2018). Implementing the production of ICL in our modelling is a valuable addition to constrain the viable models. Regardless of the satellites evolution scenario, a constant fraction of stellar mass is always transferred from mergers to the ICL. Furthermore, when including the stellar stripping in the evolution of satellites, which changes simultaneously the stellar masses of merged and surviving satellites, and merger rates and SFHs of central galaxies, the stellar mass stripped from satellites will contribute to the ICL. Indeed, we compute the ICL in decode by assuming that it is formed from the stellar mass stripped from the satellite galaxies (Section 3.2.1) and/or from the stellar mass lost during a galaxy merger. We assume that a fraction of |$\sim 20{{\ \rm per\ cent}}$| of the stellar mass during mergers involving the central Brightest Cluster Galaxy (BCG) is lost and transferred to the ICL, similarly to what is assumed in SAMs (see e.g. Contini et al. 2014). The amount of ICL might correlate with some cluster properties, such as the redshift and cluster halo mass. We discuss in this Section how the ICL evolves with the latter quantities.

The upper panel of Fig. 10 shows the redshift evolution of the fraction of ICL. Similarly to SAMs, in each cluster, we define the BCG as the most massive central galaxy that lives in the cluster halo (e.g. Tonini et al. 2012). We define the fraction of ICL as the ratio between the mass of the ICL and the total stellar mass of the cluster (i.e. the sum of the stellar masses of the BCG, surviving satellites and ICL itself). The lines in Fig. 10 show the predictions from decode for the evolving SMF model, where the ICL is formed by only stripping mechanisms (green dotted, this is an extreme scenario where mergers do not happen), only galaxy merger processes (orange dashed) and both stripping and mergers (blue solid). We compare our results with observational data from Rudick et al. (2011), Burke et al. (2012, 2015), Montes & Trujillo (2018), and Furnell et al. (2021). We selected haloes with mass between |$10^{14} \, \mathrm{M}_\odot \lesssim M_{\rm h} \lesssim 10^{15} \, \mathrm{M}_\odot$|⁠, in order to be as consistent as possible with the data sets we compare to. We stress that the ICL data are often derived in heterogeneous ways, using different methods and calibrations, but they are still all included in the same plot to provide a general term of comparison and broad guide to the models. What is evident from the full collection of the data is that despite the ICL fractions being sparse, they tend to point to averages of |$\sim 20{{\ \rm per\ cent}}$| at all cluster masses, which is a relevant constraint for models.

Consistently with the observations, we find in the top panel of Fig. 10 that the ICL gradually builds up noticeable mass below |$z\lesssim 0.7$|⁠. If only stellar stripping is considered, then the ICL is only visible at |$z\lt 0.3$|⁠. However, the predicted ICL in this model would still be limited to |$\sim 10{{\ \rm per\ cent}}$|⁠, which lies below the ensemble average of the data, possibly suggesting that the stellar mass stripped from infalling satellites by itself may not account for all the amount of observed ICL. Fig. 10 shows that the mergers contribution to the ICL is larger with respect to the stripping by a factor of |$\sim 2$|⁠. Overall, both the mergers only and stripping + mergers models better align with the general trend found in simulations and observations, i.e. a negligible amount of ICL at redshift |$z \gtrsim 1$| and gradually increasing up to |$f_{\rm ICL}\sim 20{{\ \rm per\ cent}}$| in the local Universe (see also Montes 2022, for a detailed discussion). The predicted ICL will increase/decrease by few percents and the cumulative mass growth from mergers will decrease/increase by a factor of |$\sim 0.05$| dex, by changing the fraction of mass-loss in mergers by a factor of |$10{{\ \rm per\ cent}}$|⁠. Here, once again we can appreciate the power of an SEM that can gradually build levels of complexity in a transparent way guided by the data.

The lower panel of Fig. 10 shows the fraction of ICL at redshift |$z=0$| as a function of the cluster total mass. The coloured lines show the predictions from decode, while the data points show the results from Da Rocha & Mendes de Oliveira (2005), Da Rocha et al. (2008), Kluge et al. (2021), and Poliakov et al. (2021). Exploring the correlation between the ICL and the total halo mass can give an overview on how differently the mass forms in groups and clusters of galaxies. Observations suggest that there is no correlation between the ICL and the halo mass, while simulations found contrasting results, i.e. ICL fraction constant, increasing or decreasing with halo mass (e.g. Lin & Mohr 2004; Murante et al. 2004; Purcell, Bullock & Zentner 2007; Dolag, Murante & Borgani 2010; Cui et al. 2014; Contini et al. 2018; Brough et al. 2024; Chun et al. 2024). Interestingly, our results for all the three models suggest an ICL fraction which is overall in agreement with the data, but slightly increasing from the low-mass haloes towards high-mass haloes. In our semi-empirical framework, this increase in ICL is due to the double power-law shape in the SMHM relation, which implies a higher fraction of major mergers in massive galaxies above the knee of the SMF, which are usually hosted in more massive haloes, as also discussed in section 4.2 of Paper I. Finally, the recent results of Brough et al. (2024), who computed the ICL from the Horizon-AGN, Hydrangea, IllustrisTNG, and Magneticum simulations, also suggests a fraction of ICL of |$\sim 10{\small --}20{{\ \rm per\ cent}}$| for haloes with mass |$14\lt \log _{10}(M_{\rm h}/\mathrm{M}_\odot) \lt 14.5$|⁠, broadly in line with our findings.

5 DISCUSSION

In this work, we have highlighted the intimate connection that exists in a hierarchical dark matter-dominated Universe among distinct observables, namely the SMFs, the SFHs, the merger rates, the abundances of satellites, and the level of ICL. All of these probes are all sides of the same coin, and a comprehensive model of galaxy evolution should aim at simultaneously predicting these observations. However, the challenge arises when these independent data sets may suffer from underlying inconsistencies which would prevent a safe comparison between models and data. We proved here that a semi-empirical, data-driven approach is the most suited to shed light on such possible discrepancies in the data, by using a subset of the data in input and others in output. decode, in particular, uses the SMF as input to predict the SMHM relation, and from there predicts the SFHs, merger rates, and ICL fractions using very minimal additional input assumptions and parameters.

We first highlighted in Section 4.1 that the SMF remains poorly known, even at redshift |$z=0$|⁠, with several works that proposed both very different shapes and redshift evolution for the SMF (see e.g. Tomczak et al. 2014; Bernardi et al. 2017; Davidzon et al. 2017; Leja et al. 2020; Kawinwanichakij et al. 2020; Weaver et al. 2023). The analysis in this work points out the strong dependence of the galaxy SFHs on the input SMHM relation, therefore on the SMF, as also found by Grylls et al. (2020). Our results suggest that by changing the mapping between stellar mass and halo mass, or the input SMF, the resulting mass growth history can alter significantly by several factors. In this paper, we explored two reference models, one with an evolving and another extreme one with a constant SMF up to |$z\sim 4$|⁠. A slowly evolving SMF at low masses predicts SFHs in good agreement with those retrieved from SED fitting from local data. For galaxies above the knee of the SMF instead, data suggest quite flat SFHs up to |$z\sim 3$|⁠, which would align with the SFHs from a constant SMF at high masses. We found a preference for an SMF characterized by a weak evolution in time at the faint end and by a bright end with flatter slope and significantly evolving. These features allow to predict SFHs for central galaxies with stellar mass today |$M_\star \lesssim 10^{11} \, \mathrm{M}_\odot$| that well match the observed ones (Section 4.2) and to produce enough major mergers to simultaneously reproduce the fraction of ellipticals, satellite abundances, and B/T distributions in models in which ellipticals are predominantly originating from major mergers (see Paper I and references therein).

Moreover, an SMF with a flatter and evolving high-mass end would be in agreement with the findings of significantly high volume densities of massive star-forming galaxies at |$z\gt 3$| from deep ALMA and radio surveys, likely to provide a conspicuous contribution to the high-z SFR density (HST-dark galaxies, e.g. Franco et al. 2018; Wang 2019; Gruppioni et al. 2020; Talia et al. 2021). The existence of similar numbers of high-z massive star-forming systems are a real challenge for the existing SAMs (e.g. Henriques et al. 2015) and hydrodynamical simulations (e.g. Pillepich et al. 2018b), which underestimate their number density by one to two orders of magnitude (see also Wang 2019). Recent observations with JWST extended these results to even higher redshifts, i.e. |$z=8$|⁠, finding heavily dust-obscured, massive (⁠|$M_\star \sim 10^{10} \, \mathrm{M}_\odot$|⁠), star-forming sources at |$z\sim 2{\small --}8$| with an surface densities of |$\sim 0.8 \, {\rm arcmin}^{-2}$| (Barrufet et al. 2023; Nelson et al. 2023). This suggests that an important fraction of massive galaxies may have been missing from our cosmic census at |$z\gt 3$| all the way to the epoch of reionization.

We also found that the main sequence as observed by several surveys is higher in normalization or flatter in slope at low masses than decode’s predictions, and are not always in line with their observed stellar masses (e.g. Speagle et al. 2014; Whitaker et al. 2014; Tomczak et al. 2016; Popesso et al. 2023). Leja et al. (2015) extensively discussed this problem and, following a continuity equation approach to evolve galaxy stellar mass growths forwards in time, also found that a main sequence with steeper low-mass end is more suitable to fit the observed SMFs. Similar results are put forward by Leja et al. (2019b), where they evolved the galaxy SMF backwards in time with the SFRs from the 3D-HST catalogues. Similarly, by deriving the SFRs directly from the stellar mass assemblies, which by design fit the mass function, we also found a main sequence which is lower in normalization and steeper in slope than those suggested by some observations cited above, whilst consistent with those from simulations (see e.g. Donnari et al. 2019). Moreover, we found that the main sequence presents a break and drops towards high stellar masses, especially at lower redshifts. In our semi-empirical modelling, the shape of the SFR–|$M_\star$| relation is a direct by-product of the SMHM relation. As the latter is a double power law, also the main sequence tends to present a similar shape. However, the break in the SFR–|$M_\star$| relation changes significantly with redshift, moving from |$2\times 10^{11} \, \mathrm{M}_\odot$| at |$z\ge 2$|⁠, to |$3\times 10^{10} \, \mathrm{M}_\odot$| at |$z\lt 1$|⁠. Therefore, the main sequence does not present any significant break or flattening up to high stellar masses of a few |$10^{11}\, \mathrm{M}_\odot$| at high redshifts, with a flattening appearing only below |$z\lt 1$|⁠, more or less pronounced depending on the shape of the input SMHM relation. This behaviour in the predicted main sequence is visible even in the extreme case where mergers are not included in the model (see Fig. 7). When mergers are self-consistently included, they become the dominant process in the galaxy stellar mass assembly leading to a further, more marked drop in the SFR, irrespective of the input SMF model. On the other hand, towards higher redshifts galaxies are still in an active star-forming phase where merger’s contribution is less with respect to lower redshifts, and the flattening in the main sequence is barely visible. Our results are in line with the recent works that found a flattening in the slope of the main sequence at high stellar masses, e.g. Leja et al. (2022) and Popesso et al. (2023), who suggested that such a curvature in the SFR–|$M_\star$| relation is most likely due to the suppression of the SFR in galaxies.

Moreover, from the tests on the TNG simulation and from our analysis in Section 4.3, we found consistently that the less massive satellite galaxies continue to form stars after their infall before being rapidly quenched (see Shi et al. 2020, for a detailed discussion). The addition of this feature allowed us to both reproduce the TNG’s and observed satellite abundances. Finally, in this work, we also found evidence for the need of both mergers in addition to the stellar stripping to form the ICL at all stellar masses. This finding is in line with what is found by several observational and theoretical works (see e.g. Puchwein et al. 2010; Rudick et al. 2011; Burke et al. 2012; Contini et al. 2014, 2018; Montes & Trujillo 2018, 2019; Contini 2021). Indeed, Contini et al. (2014) showed that mergers can form up to |$\sim 20{{\ \rm per\ cent}}$| of the mass of the ICL in massive clusters. Furthermore, many works, such as Montes & Trujillo (2018) and Contini & Gu (2020), suggested that the stellar stripping forms gradually more mass in the ICL during redshift evolution. The very recent work of Contini et al. (2023) showed that the fraction of ICL also slightly depends on the halo mass, increasing with the latter and then staying steady at |$\sim 0.35$|⁠, similarly to what we found. These findings are based on the assumption that |$20{{\ \rm per\ cent}}$| of the mass of the merging satellites is lost during mergers, which is the value assumed in most of the SAMs in the literature and allows a good compromise between mergers, SFHs and ICL.

6 CONCLUSIONS

In this paper, we have presented a holistic perspective encompassing galaxy SFHs, merger histories, satellite abundances, and ICL. All of these observables, should in principle be linked together in a hierarchical, dark matter-dominated Universe where galaxies grow via mergers and star formation, and live in multiple environments, from the field to clusters. To probe the connection and self-consistency among these distinct data sets, we made use of decode in a |$\Lambda$|CDM hierarchical framework, where dark matter assemblies and merger trees are converted into galaxy mass growth and merger histories via the input SMHM relation. The SFHs are then derived from the difference between the latter two quantities. We started by testing decode on a self-consistent and comprehensive model of galaxy formation and evolution, namely the TNG simulation. We found that by using in input TNG’s SMHM relation, we were able to simultaneously reproduce the mean star formation and merger histories of TNG’s central galaxies, as well as its satellite SMF. We then turned to apply decode to real data. To this purpose, we used in input two models of the SMHM relation based on the abundance matching between the SMF and the HMF. In one model, we assumed the SMF to be constant up to |$z=1.5$|⁠, and evolving after that, and in the other we assumed the SMF to evolve significantly from |$z=0$|⁠. The choice of these two models broadly bracket the most recent observations which are still not agreeing on the exact shape and evolution of the galaxy SMF, as discussed in Section 4.1. From these two models we derived via decode a variety of observables which we compared with a number of independent data sets.

Our main results can be summarized as follows:

• An SMHM relation characterized by a weakly evolving low-mass end produces SFHs in relatively good agreement with those inferred from SED fitting in local surveys such as MaNGA and GAMA, whilst fast evolving low-mass end of the SMF, suggested by some observational groups, are disfavoured (Fig. 6).

• For more massive galaxies above |$M_\star \gt 3 \times 10^{10}\, \mathrm{M}_\odot$| SED-based estimates of the SFHs do not currently agree with each other. Reproducing the approximately flat stellar mass growth histories inferred from GAMA requires a nearly constant SMF at all stellar masses and up to |$z\sim 2$|⁠, in line with the recent suggestions of an increased star formation efficiency in massive galaxies at |$z\gt 2$| from JWST data (e.g. Atek et al. 2023; Endsley et al. 2023; Labbe et al. 2023; Nelson et al. 2023; Weibel et al. 2024).

• The evolution of satellites does not affect the SFHs of central low-mass and intermediate-mass galaxies since the contribution from mergers is small.

• The merger histories of massive galaxies above |$M_\star \gt 10^{11} \, \mathrm{M}_\odot$| have a relatively weak dependence on the type of post-infall satellite evolution implemented in the model, while the integrated SFHs and specific star formation rates tend to change by a factor of at most |$\sim 2{\small --}3$| depending on the assumed satellite evolution after infall (Fig. 7).

• The main sequence SFR–|$M_\star$| relation implied by the two SMF models explored here tends to be lower in normalization by a factor |$\sim 2{\small --}3$| with respect to observations. It also shows clear signs of a flattening at |$M_\star \gt 3\times 10^{10} \, \mathrm{M}_\odot$| but only at |$z\lesssim 1$|⁠, which is independent of the merger rate but a natural byproduct of the break in the SMHM relation. Including mergers in the models tends to further sharpen the drop in the main sequence above |$3\times 10^{10}{\small --}10^{11} \, \mathrm{M}_\odot$|⁠, as more massive galaxies tend to grow proportionally less via star formation (Fig. 8).

• The satellites can be considered as failed mergers, and thus their abundances, which depend on the same SMHM in input, represent a valuable complementary probe to test galaxy evolution models. We found that the SDSS local satellites SMF can be nicely reproduced by both the SMHM relations explored here, and including residual star formation in satellites after infall improves the match with the low-mass end of the satellite SMF, especially in more massive parent haloes.

• The ICL is another probe of the assembly history of galaxies. Allowing in our reference models for |$\sim 20{{\ \rm per\ cent}}$| of stellar mass loss during mergers as well as stellar stripping in satellites within clusters, provides a good match to current constraints on the ICL at different redshifts and parent halo masses, irrespective of the input SMHM relation, and with minimal impact on the SFH/merger histories of galaxies (Fig. 10).

In this work, we have put forward a flexible and efficient data-driven approach to probe the self-consistency within a hierarchical, dark matter-dominated Universe, of the SMF with other key independent observational probes such as the star formation rates, the satellite abundance, and the ICL, which will be applied to the imminent data release from extragalactic missions. The advent of new high-quality data from ongoing surveys, such as JWST (Gardner et al. 2006) and Euclid (Amiaux et al. 2012), which will provide self-consistent determinations of the galaxy stellar mass function, will be extremely beneficial for building a more comprehensive and complete understanding of galaxy evolution.

ACKNOWLEDGEMENTS

This work received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 860744 for the BiD4BESt project (Coordinator: F. Shankar). MA acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG) through an Emmy Noether Research Group (grant number NE 2441/1-1).

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES

APPENDIX A: TESTING THE ABUNDANCE MATCHING

Fig. A1 shows the SMHM relation from our abundance matching using TNG’s SMF input, compared to that directly derived from the simulation (see discussion in Section 3.4).

APPENDIX B: EFFECT OF SATELLITES STAR FORMATION ON CENTRAL GALAXIES GROWTH

Fig. B1 shows the mean SFH in the stellar mass bin of |$M_\star 10^{10.5}\, \mathrm{M}_\odot$| at redshift |$z=0$| for the frozen and star-forming satellites scenarios, showing that the effect of the star formation in satellites on the SFH of the central galaxies is minimal.