A flexible subhalo abundance matching model for galaxy clustering in redshift space

MNRAS 508, 175–189 (2021) https://doi.org/10.1093/mnras/stab2560 Advance Access publication 2021 September 11 A flexible subhalo abundance matching model for galaxy clustering in redshift space S. Contreras ,1‹R. E. Angulo 1,2and M. Zennaro 1 1Donostia International Physics Center (DIPC), Manuel Lardizabal Ibilbidea, 4, E-20018 Donostia, Gipuzkoa, Spain 2IKERBASQUE, Basque Foundation for Science, E-48013 Bilbao, Spain Accepted 2021 September 3. Received 2021 August 24; in original form 2020 December 11 ABSTRACT We develop an extension of subhalo abundance matching (SHAM) capable of accurately reproducing the real and redshift-space clustering of galaxies in a state-of-the-art hydrodynamical simulation. Our method uses a low-resolution gravity-only simulation and it includes orphan and tidal disruption prescriptions for satellite galaxies, and a flexible amount of galaxy assembly bias. Furthermore, it includes recipes for star formation rate (SFR) based on the dark matter accretion rate. We test the accuracy of our model against catalogues of stellar-mass- and SFR-selected galaxies in the TNG300 hydrodynamic simulation. By fitting a small number of free parameters, our extended SHAM reproduces the projected correlation function and redshift-space multipoles for number densities 10−3−10−2h3Mpc−3,atz=1 and z=0, and for scales r∈[0.3 −20]h−1Mpc. Simultaneously, the SHAM results also retrieve the correct halo occupation distribution, the level of galaxy assembly bias, and higher order statistics present in the TNG300 galaxy catalogues. As an application, we show that our model simultaneously fits the projected correlation function of the SDSS in three disjoint stellar mass bins, with an accuracy similar to that of TNG300 galaxies. This SHAM extension can be used to get accurate clustering prediction even when using low and moderate-resolution simulations. Key words: galaxies: evolution – galaxies: formation – galaxies: haloes – galaxies: statistics – cosmology: theory – large-scale structure of Universe. 1 INTRODUCTION With a new generation of surveys soon to start collecting data, a new generation of theoretical predictions for galaxy clustering is needed. These ‘mock’ galaxies are essential to prepare for these observations, estimate sources of uncertainty, and exploit them scientifically. Due to the large volume and statistical power of upcoming surveys, this new generation of mocks will need to be accurate and computationally efficient. In the CDM paradigm, galaxies form from the gas that is captured in the gravitational potential well of dark matter haloes (White & Rees 1978). As this gas cools, it falls to the centre of haloes and settles into a disc, to then cool and fragment into stars. While the evolution of a halo is somewhat simple, growing in mass through accretion and mergers determined mostly by gravitational interactions, the evolution of galaxies is much more complex, with a large number of additional astrophysical processes involved. Hydrodynamic simulations are probably the most realistic method to model these processes, however, they can take up to hundreds of millions of CPU hours to execute as they require high spatial and mass resolution. This implies that most of these simulations have only reached volumes of the order of 100 h−1Mpc aside (e.g. EAGLE, Schaye et al. 2015; Illustris,Vogelsberger et al. 2014; HorizonAGN, Dubois et al. 2014; Simba, Dav´ eetal.2019, simulations). E-mail: [email protected] Other techniques which simplify the modelling of some astrophysical processes, like semi-analytical models of galaxy formation (SAMs; Henriques et al. 2015; Croton et al. 2016; Lacey et al. 2016; Stevens, Croton & Mutch 2016; Lagos et al. 2018)orsemiempirical models such as EMERGE (Moster, Naab & White 2018) or The Universe Machine (Behroozi et al. 2019) can be run in postprocessing over Gpc scale gravity-only simulations in hundreds of CPU hours (only run, not calibrate them). While much faster than hydrodynamic simulations, this is still too slow when thousands or even millions of mocks are needed to, e.g. scan high dimensional parameter spaces or compute covariance matrices. An additional limitation is that these approaches require underlying simulations with considerably high force and mass resolution. Although the evolution of galaxies and haloes is different, we expect a strong correlation between their properties, which motivates ‘empirical models’ as an alternative to creating fast and computationally efficient mocks. These models ‘paint’ galaxies in haloes or subhaloes using theoretically motivated relations between the galaxies and their host (sub)haloes. These methods can be run in a few minutes over large cosmological volumes, and have proved to successfully reproduce several galaxy observables, such as the galaxy two-point correlation function. One example of an empirical models is the SubHalo Abundance Matching technique (SHAM; Conroy, Wechsler & Kravtsov 2006; Reddick et al. 2013;Chaves- Montero et al. 2016; Lehmann et al. 2017; Dragomir et al. 2018). This technique is based on the idea that most massive subhaloes should host the most massive/luminous galaxies. The algorithm matches, with some scatter, the selected subhalo property to the expected C 2021 The Author(s) Published by Oxford University Press on behalf of Royal Astronomical Society Downloaded from https://academic.oup.com/mnras/article/508/1/175/6368853 by guest on 04 November 2025 176 S. Contreras, R. E. Angulo and M. Zennaro stellar mass/luminosity function. While simple, the model is capable of reproducing the galaxy clustering in complex galaxy formation models, such as hydrodynamic simulations (e.g. Chaves-Montero et al. 2016). The main limitation of SHAM is that it assumes a direct relationship between the fate of a dark matter satellite and the putative galaxy it might host. This assumption does not hold in detail, as the presence of a galaxy would modify e.g. dynamical friction time-scales and/or the resilience of a subhalo against tidal disruption. Furthermore, in relatively low resolution simulations, a subhalo might be stripped below the resolution limit of the simulation, but the hosted galaxy is still expected to survive as a satellite. Another shortcoming of SHAM is that it is limited to properties expected to correlate strongly with dark matter mass, so it cannot make predictions for star-formation rates, which would be relevant for upcoming galaxy surveys targeting emission line galaxies. In this paper, we develop an extension to the standard subhalo abundance matching technique that, with only three free parameters (five if also including assembly bias) and a low-resolution dark matter only (DMO) simulation, can reproduce the clustering of stellar massselected galaxies in real and redshift space with high precision. This extension includes the use of orphan satellite subhaloes/galaxies tracking structures below the resolution level of the simulation; a disruption mechanism to account for satellite galaxies that should have been disrupted to form part of the intracluster medium, or in which their stellar mass has decreased so much that they should no longer be part of the galaxy selection; and a flexible amount of galaxy assembly bias to account for any possible difference between correlation with the large-scale environment in SHAM and our target sample. To our knowledge, this is the first work that shows that abundance matching techniques can reach a high level of precision in redshift space. The only other work we know that has studied the performance of an abundance matching technique with redshift space clustering has been Guo et al. (2016a). There, the author looked at another extension of the abundance matching technique, which is capable of reproducing the observed monopole and quadrupole of the correlation function with a good degree of precision. By focusing their comparison on an observational sample, it is not possible to identify if the differences found are caused by the redshift space distortion modelling or by any other systematic, like, for example, by assuming a different cosmology to the one of the Universe. Along with these improvements, we develop a method capable of predicting the clustering of SFR-selected galaxies. This method assumes that the SFR of each galaxy is proportional to the dark matter accretion rate of its host halo, and that is modulated by an efficiency that depends on the mass of its host halo. Same as the stellar mass approach, this method is similar to that employed by the semi-empirical model EMERGE but adapted to not use the full merger tree on its computation but instead properties readily available in SHAM. These approaches should also be faster and use fewer free parameters while keeping a similar precision when reproducing the galaxy clustering in real and redshift space. To test the performance of our SHAM extension, we use the Illustris TNG300 magnetohydrodynamic simulation (Nelson et al. 2018), to our knowledge, the largest publicly cosmological hydrodynamic simulation available today (L=205 h−1Mpc). For the stellar mass-selected mocks, we find that our method is capable of reproducing almost perfectly the real and redshift-space clustering, as well as the galaxy assembly bias level, halo occupation distribution (HOD), and higher order statistics as quantified by the k-nearest neighbour cumulative distribution functions (kNN-CDF). Similarly, for the SFR-selected sample, we find that we are also able to reproduce the same galaxy statistics, albeit with somewhat less accuracy, particularly for the level of assembly bias. The differences in the galaxy clustering between the star-forming mock and the TNG300 are still much smaller than those found among different galaxy formation models (Contreras et al. 2013). We expect this new generation of mocks to be useful for the development and understanding of current and future galaxy surveys. The outline of this work is as follows: In Section 2, we present the simulations used in this work and introduce the SHAM technique. Our modifications to SHAM regarding the treatment of satellite galaxies are discussed in Section 3 and those regarding galaxy assembly bias in Section 4. In Section 5, we present our starformation-rate modelling. We compare the two-point clustering of our mock galaxies with that of the TNG300 in Section 6 and with observations in Section 7. To finalize, we present our conclusions in Section 8. Unless otherwise stated, the units in this paper are h−1Mfor masses, h−1Mpc for distances, km s−1for the velocities and Myr−1 for star formation rates. 2 SIMULATIONS AND THE STANDARD ABUNDANCE MATCHING MODEL 2.1 The TNG300 To test the accuracy of our mocks, we will compare their clustering against those of galaxies in the TNG300 magnetohydrodynamic simulation, the largest of ‘The Next Generation’ Illustris Simulations suite, and, to our knowledge, the largest volume high-resolution hydrodynamic simulation publicly available (Marinacci et al. 2018; Naiman et al. 2018;Nelsonetal.2018; Pillepich et al. 2018; Springel et al. 2018). The TNG300 simulation follows a periodic box of 205 h−1Mpc (∼300 Mpc) a side. It used 25003dark matter particles and gas cells, implying a baryonic mass resolution of 7.44 ×106h−1Mand of 3.98 ×107h−1Mfor dark matter. The simulations were carried out using the AREPO code (Springel 2010) adopting cosmological parameters consistent with recent analyses (Planck Collaboration XIII 2016).1While main results of this paper are only shown at z=0, we tested that they are also valid at z =1 (Appendix B). The outputs of the simulation are publicly available at the TNG project webpage.2 We built galaxies catalogues by selecting the most massive or starforming galaxies such that their abundance is equal to a number density of n=0.01,0.00316 ,&0.001 h3Mpc−3, equivalent to a minimum value of the stellar mass of 8.26 ×109, 2.93 ×1010, &6.46 ×1010 h−1Mand a minimum SFR of 0.468, 1.493, & 3.033 Myr−1, respectively. We defined the stellar mass of a galaxy as the sum of the mass of all particles within twice the stellar halfmass radius. We also test defining the stellar mass of the galaxies as the sum of the masses of all particles of a subhalo, finding similar results. We define the SFR of a galaxy as the sum of the individual star formation rate of all gas cells in its subhalo. We note that although the TNG300 simulation has shown to agree with many observables (Springel et al. 2018), it is not the aim of this work to create mocks that only resemble this simulation, but to create a flexible model that can describe a broader range of feasible galaxy formation physics and their correlation with 1dm =0.3089, b=0.0486, σ8=0.8159, ns=0.9667, and h=0.6774 2https://www.tng-project.org/ MNRAS 508, 175–189 (2021) Downloaded from https://academic.oup.com/mnras/article/508/1/175/6368853 by guest on 04 November 2025 A flexible SHAM model 177 cosmological parameters. We chose a hydrodynamic simulation as a benchmark since it models baryons and dark matter jointly, unlike other approaches that can reproduce galaxy clustering (such as HODs, semi-empirical models, or semi-analytical models). This adds an extra challenge for mocks created in DMO simulations, allowing us to perform a stringent test of the performance of our SHAM galaxies. In addition to the TNG300, we also use the TNG300-1-Dark, TNG300-2-Dark, and TNG300-3-Dark simulations. These are DMO simulations employing the same initial white noise field as the TNG300 run. The TN300-1-Dark has the same number of dark matter particles as the TNG300 (25003) whereas the TNG300-2- Dark and TNG300-3-Dark have factors of 23and 43fewer particles (i.e. 12503and 6253particles, respectively). For some comparisons, we will use the cross-match of TNG galaxies and subhaloes in the TNG300-1-Dark, available on the Illustris TNG website. 2.2 TNG300-3 mimic The main dark matter simulation we will employ is analogous to the TNG300-3-Dark (same volume and initial conditions as the TNG300, but 43fewer particles), but with additional properties computed for its subhaloes as we describe below. We refer to this simulation as the TNG300-3-mimic, and highlight that by using a simulation with the same initial condition as the TNG300 we can perform more accurate comparisons as it reduces any difference between the model caused by cosmic variance. The simulation was carried out with an updated version of L- Gadget3 (Angulo et al. 2012), a lean version of GADGET (Springel et al. 2005) used to run the Millennium XXL simulation and the Bacco Simulations (Angulo et al. 2020). This version of the code allows an on-the-fly identification of haloes and subhaloes using a Friend-of-Friend algorithm (FOF; Davis et al. 1985) and an extended version of SUBFIND (Springel et al. 2001). Our updated version of SUBFIND can better identify substructures by considering the information of its past history, while also measuring properties that are non-local in time such as the peak halo mass (Mpeak), peak maximum circular velocity (Vpeak), infall subhalo mass (Minfall), and mass accretion rate, among others. We note that we employ the same numerical accuracy parameters as the Bacco simulations (Angulo et al. 2020). And, as a validation test, we have compared the z=0 mass power spectrum of our simulation to that measured in the TNG300-3-Dark finding differences below 0.8 per cent on all scales. 2.3 Subhalo abundance matching In its most basic incarnation, SHAM assumes a one-to-one mapping between the mass of subhaloes and the stellar mass or luminosity of hypothetically hosted galaxies. More recent implementations of SHAM include scatter to this mapping and employ other subhalo properties, such as Minfall,Mpeak,Vmax,orVpeak, but still, it is assumed that the position and fate of galaxies are determined by the evolution of dark matter subhaloes. This, however, cannot be true in detail since subhaloes might be artificially disrupted in N-body simulations (because of numerical artefacts and/or limited mass resolution). More generally, the presence of a galaxy can modify the subhalo density profile, affecting its orbit, tidal stripping, and eventual disruption. This puts a limit to the accuracy of SHAM, which we will illustrate next. Fig. 1compares the stellar mass of TNG300 galaxies and the value of Vpeak, the peak value of the maximum circular velocity of a halo Figure 1. The stellar mass of galaxies in the TNG300 simulation as a function of Vpeak of the matching subhaloes in the TNG300-1-Dark simulation. The shaded regions represent 1, 2, & 3 σof the distribution, and the dashed red and blue lines show the median for samples of central and satellite galaxies, respectively. The horizontal dotted lines indicate the three number densities used in this work, as labelled. throughout its life (not to be confused with its value at the moment of infall, Vinfall), of the corresponding cross-matched subhalo in the TNG300-1-Dark. Note that Chaves-Montero et al. (2016) showed that in the EAGLE simulation (Schaye et al. 2015)Vpeak is one of the properties that best predict the value of stellar mass. Indeed, we can see that there is a strong correlation between these two properties in TNG300, with a mean relation being almost identical for centrals and satellite galaxies, validating a key SHAM assumption. 3 Chaves-Montero et al. (2016) also found that SHAM using Vpeak was capable of reproducing the clustering of stellar mass selected galaxies in EAGLE with very high accuracy. However, these authors used a very high resolution dark-matter simulation (which would imply artificial disruption of subhaloes was negligible), and the accuracy on small scales varied depending on the assumptions regarding star formation in satellites and their tidal disruption. We illustrate the impact of mass resolution for SHAM in Fig. 2.We show the projected correlation function, wp(rp), and the monopole, quadrupole, and hexadecapole of the redshift space correlation function (ξ=0(r), ξ=2(r), ξ=4(r), respectively) for stellar-mass selected galaxies at z=0. We use an integration depth of 40 h−1Mpc for wp, and model the redshift space distortion effect by modifying the position of the galaxies along the z-axis. Solid lines show the results for TNG galaxies, whereas coloured lines show the SHAM predictions as computed using the TNG300-1-Dark, TNG300-2- Dark, & TNG300-3-Dark simulations. Note we display results for three different number densities. We see that SHAM in the TNG300-1-Dark agrees remarkably well with the TNG measurements for the densest galaxy sample, but it 3Note that, although not shown here, we found a similar scatter between Mpeak and stellar mass but with a significantly larger difference between centrals as satellite galaxies, which further motivates the use of Vpeak instead of Mpeak (see also the discussion in Campbell et al. 2018). MNRAS 508, 175–189 (2021) Downloaded from https://academic.oup.com/mnras/article/508/1/175/6368853 by guest on 04 November 2025 178 S. Contreras, R. E. Angulo and M. Zennaro Figure 2. The projected correlation function (wp(rp), leftmost panel), the monopole (ξ=0(r), second panel), quadrupole (ξ=2(r), third panel), and hexadecapole (ξ=4(r), rightmost panel) of the galaxies of the TNG300 simulation (black solid line) and the subhaloes of the TNG300-1-Dark, TNG300-2-Dark, & TNG300- 3-Dark dark matter only simulation (cyan, magenta, and yellow dashed lines, respectively). The galaxies/subhaloes samples have a number density of n=0.01, 0.00316, and 0.001 h3Mpc−3selected by stellar mass (for the TNG300) and Vpeak (for the dark matter only simulations). For visualization purposes, the measurements for the two lowest number densities have been shifted in the y-axis with the lowest number density at the top – wp(rp)andξ=0(r) have been added 0.7 and 1.4, and ξ=2(r) and ξ=4(r) 20 and 40. overpredicts the clustering for lower number densities. Conversely, SHAM in the lowest resolution simulation, TNG300-3-Dark, agrees very well with the sparsest sample but instead heavily underpredicts the clustering of the densest sample. Note that these problems are relatively confined to small scales in real space, but can affect even large scales in the case of the quadrupole and hexadecapole. All this is simply a consequence that in standard SHAM, the abundance of satellites depends on the numerical resolution of the parent dark matter simulation which might result in an excess or lack of satellite galaxies. For instance, galaxies in the TNG300 can be affected by stripping of the stellar mass, while subhaloes never decrease their Vpeak, producing an overabundance of satellite galaxies (Smith et al. 2016). Overall, this highlights an important limitation of SHAM and motivates our development of an extension to SHAM that addresses these problems. We use the subhalo abundance matching technique (SHAM) with Vpeak as a starting point of our mocks. First, we match the value of Vpeak of the dark matter-only simulation to the stellar mass function of the TNG300 simulation. Since most of our analysis is done using cuts in number density, the results are mostly independent of the stellar mass function used. We then (i) add a free parameter controlling the scatter (σlogM) between the stellar mass of galaxies and their Vpeak, (ii) explicitly model the level of assembly bias (following the procedure presented in Contreras, Angulo & Zennaro 2021), and (iii) regulate the disruption of satellites and the star formation rate with recipes inspired by semi-analytic galaxy formation models and by the semi-empirical model EMERGE (Moster et al. 2018). In the next sections, we will discuss our prescriptions and show the impact that the free parameters of our model have on the predicted redshift space clustering of galaxies. 3 THE SATELLITE TREATMENT IN SHAM In this section, we introduced two SHAM modifications related to the treatment of satellite galaxies. In Section 3.1, we describe how we follow subhaloes after they can no longer be identified by our simulations (also known as ‘orphan suhaloes’), and in Section 3.2 we explain our implementation to account for tidal disruption and stripping of stars. 3.1 Orphan implementation Due to the finite mass resolution of dark matter simulations, subhaloes can no longer be identified after they have lost a given fraction of their mass due to tidal stripping. To improve convergence among simulations with different resolutions, it is necessary to include ‘orphan’ subhaloes/galaxies. These are structures with a known progenitor which we expect to still exist in the halo (i.e. that should not have merged with the central structure or completely destroyed via tidal forces). Orphan substructures are normally used in semi-analytical models of galaxy formation (a.k.a. SAMs) to reproduce galaxy clustering (De Lucia & Blaizot 2007; Guo et al. 2011). While some SHAM prescriptions have included orphans in the past (e.g. Moster, Naab & White 2013), most do not include this feature, because it adds an additional level of complexity to the model and to the underlying dark matter simulation. While uncommon, Guo & White (2014) found that not using an orphan prescription in the SHAM will cause some differences in the galaxy clustering prediction, even when using high-resolution simulations. In this work, we include a flexible amount of orphan galaxies by tracking the most bound particle of subhaloes with no known descendant. We assume a subhalo to be merged with their central structure only when the time since accretion exceeds the dynamical friction time-scale of the subhalo. To compute the dynamical friction time, we use a modified version of the expression given by (Binney & Tremaine 1987, equation 7.26) tdyn.fric.=1.17 tmerger d2 host Vhost Mhost/1013 h−1M1/2 Gln(Mhost/Msub +1) Msub ,(1) where tmerger is a free parameter that effectively regulates the number of orphan galaxies; dhost is the distance of the subhalo to the centre of its host halo; vhost is the virial velocity of the host halo; Mhost is the virial mass of the host halo, and Msub is the subhalo mass. The only difference with the original equation (besides the free parameter) is the term M1/2 host, which we find improves our predictions for the number of satellite galaxies in high halo masses. MNRAS 508, 175–189 (2021) Downloaded from https://academic.oup.com/mnras/article/508/1/175/6368853 by guest on 04 November 2025 A flexible SHAM model 179 Figure 3. The impact of the free parameters in our extended SHAM for stellar-mass selected galaxies with a number density of n=0.01 h3Mpc−3at z=0. The first column displays the projected correlation function, wp(rp); whereas the second, third, and fourth columns, respectively, display the monopole ξl=0(r), quadrupole ξ=2(r), and hexadecapole ξ=4(r) of the redshift space correlation function. In each row, we systematically vary a different SHAM parameter, σlogM, tmerger,andfs, while keeping all the other parameters fixed. Red to blue colours indicate changes from low to high parameter values, as indicated by the legend. 3.2 Tidal disruption One of the most basic assumptions in SHAM is that the relation between stellar mass and subhalo property used (Vpeak in our case) is constant through time and the same for central and satellites. While the Vpeak–stellar mass relation is indeed very similar for centrals and satellites (c.f. Fig. 1), it is not identical. In Smith et al. (2016), the authors showed that satellite galaxies can decrease their stellar mass to the point of disappearing into the intracluster medium. This effect should remove some satellite galaxies of a given sample (most of them located near the centre of a halo). To account for this effect, we follow Moster et al. (2018)and assume that all the galaxies in our mock hosted by subhaloes with current mass below a fraction fsof their maximum subhalo mass during its evolution, Mpeak, are disrupted, Msub <fsMpeak.(2) In Chaves-Montero et al. (2016) the authors found that the clustering of SHAM carried out over a dark matter simulation with the same resolution as EAGLE (a 67.7 h−1Mpc high-resolution hydrodynamic simulation) yielded a larger clustering compared to stellar-mass selected EAGLE galaxies. Similar results were found for the TNG300 by Contreras et al. (2021) when comparing SHAM mocks using Vpeak. This is because Vpeak never decreases, contrary to the stellar mass of galaxies that can do so via tidal stripping. By including a disruption parameter, SHAM can better reproduce the galaxy clustering even in high-resolution simulations (see Appendix A. for more details). 3.3 Parameter variations We now explore the impact that these recipes and free parameters have in the predicted redshift space clustering of SHAM galaxies. Fig. 3shows the multipoles and the projected correlation function of galaxies at z=0 with a number density of 0.01 h3Mpc−3in our TNG-300-mimic simulation. Coloured lines in each row show SHAM predictions varying either the scatter between stellar mass and Vpeak,σlogM; the merger time-scale, tmerger; or tidal disruption, fs; while keeping all other parameters fixed. For comparison, we also display the results for TNG galaxies as black lines. First, we see that the σlogM parameter has a small but systematic effect in the clustering. As its value increases, smaller subhaloes will be included in the selected sample which will progressively reduce the mean halo mass and large-scale bias. Note that the scatter we empirically measure from Fig. 1is σlogM ∼0.125. We find the effect of this parameter to be small at all number densities. Only when setting the rest of the parameters to zero, we find an increase in the effect of σlogM at the sparsest number density, but even in this case, its effect is less significant than that of any other parameter of the model. Next we consider tmerger which has a significantly stronger impact. Large values of tmerger essentially imply that every subhalo ever accreted will still be present in the simulation – this is thus an estimate for a simulation with infinite mass resolution. Consequently, there will be a large population of satellite galaxies in high-mass haloes, increasing its contribution and the clustering amplitude. In contrast, lower values of tmerger implies orphan subhaloes will be destroyed immediately, thus lowering the clustering and approaching the case of no orphan subhaloes. Note that this parameter brackets the two limiting cases of no artificial disruption in N-body simulations and of every subhalo destroyed being due to numerical reasons. From all the parameters introduced in this model, tmerger is the only one that can increase the galaxy clustering at small scales and is key to reproduce the galaxy clustering on low-resolution simulations. Finally, as expected, larger (smaller) values of the tidal parameter fsmeans fewer (more) satellite galaxies and thus, a lower (higher) clustering. Whereas this effect appears negligible compared to the effect of orphan subhaloes and tmerger, this parameter becomes more relevant when dealing with high-resolution simulations. We remind the reader all the changes show before are for a fixed number density, meaning that the effect on the galaxy clustering is not only due to the galaxies we are directly removing/adding, but also those that are included/excluded from the sample to keep the total number of galaxy constant. Also, please notice that we provide this as an illustration of the parameters of the model and that in practice we expect some degree of degeneracy between all of them. MNRAS 508, 175–189 (2021) Downloaded from https://academic.oup.com/mnras/article/508/1/175/6368853 by guest on 04 November 2025 180 S. Contreras, R. E. Angulo and M. Zennaro 4 ASSEMBLY BIAS The concept of halo assembly bias (also known as secondary bias; Mao, Zentner & Wechsler 2018) was introduced by Sheth & Tormen (2004) and Gao, Springel & White (2005) and is defined as the difference in the clustering of haloes of the same mass but of different secondary property (e.g. halo age, concentration, spin). Halo assembly bias can also propagate to galaxy clustering when haloes of the same mass but different secondary properties have a different halo occupation (effect known as occupancy variation; Artale et al. 2018;Zehavietal.2018). This difference in the galaxy clustering is known as galaxy assembly bias (Croton, Gao & White 2007) and can cause differences in the correlation function of up to a 20 per cent for stellar mass-selected samples, depending on the redshift and number density of the galaxy sample in SAMs and hydrodynamic simulations (Contreras et al. 2019,2021). In Chaves-Montero et al. (2016), the authors found that the galaxy assembly bias signal predicted by SHAM can be up to 50 per cent lower than that measured in EAGLE. In Contreras et al. (2021) we extended this work and found that hydrodynamic simulations, semi-analytical models, and SHAM predict different levels of galaxy assembly bias and that its evolution with redshift and number density is also different. Since it is unclear that the level of assembly bias in SHAM is a robust prediction, here we include an extension so that any amount of galaxy assembly bias signal can be added. While the standard SHAM assumes that the scatter between the stellar mass and Vpeak (or any other subhalo property used) is completely random (i.e. it does not depend on any secondary property of the subhalo), our method adds a correlation with the individual large-scale bias of subhaloes (measured with the method developed by Paranjape, Hahn & Sheth 2018). This is done for central and satellite galaxies independently, forcing conservation of the satellite fraction at all number densities (see Section 4 of Contreras et al. 2021 for more details). To imprint a given level of galaxy assembly bias, the method uses two free parameters (one for central galaxies, Ac, and one for satellite galaxies, As) that control the correlation between stellar mass and individual bias of each subhalo at a given constant Vpeak. Using two independent parameters is necessary since the level of assembly bias may be different between centrals and satellites. Please notice that the value of the assembly bias parameters will be different from those in other models (e.g. Zentner et al. 2019). This is because our method is based on SHAM and not HODs, as is the case for most of the other methods available. Also, our model uses the individual large-scale bias of subhaloes and not the concentration of the host haloes (which can only imprint a smaller amount of assembly bias (Hadzhiyska et al. 2020; Xu, Zehavi & Contreras 2021).) These free parameters can have a value between 1 (perfect correlation between stellar mass and bias-per-object, which will mean a maximum assembly bias signal), and −1 (perfect anticorrelation between stellar mass and bias-per-object, which will mean a minimum assembly bias signal), with 0 being an uncorrelated signal, meaning the sample will have the same assembly bias signal as the standard SHAM. Although not shown here, we kindly refer to fig. 8 of Contreras et al. (2021) where we showed that the difference in the galaxy clustering on large scales for a SHAM with maximum and minimum assembly bias is around a factor three. 5 STAR FORMATION RATE Most of the empirical models available to create mocks are focused on galaxy samples selected by stellar mass-like properties (e.g. HOD, SHAM, SCAM), and comparatively much less effort has been made to model star-forming galaxy samples. This is because most of these models relate a given (sub)halo property to a galaxy property, and while these relationships are quite simple for stellarmass selected samples (i.e. the stellar mass is a good proxy for how massive a (sub)halo is) they are much more complicated for SFRs. More sophisticated models, such as semi-analytical models of galaxy formation, show no good agreement in how classic (sub) halo properties correlate with SFR-like galaxy properties (Contreras et al. 2013,2015). To account for this non-trivial relation between SFR-like galaxy properties and (sub) halo properties, empirical models have become more complicated when trying to reproduce these galaxies. Whereas a standard HOD (Zheng et al. 2005) can reproduce the galaxy clustering of a stellar mass selected sample with reasonable accuracy using only five free parameters, SFR-selected galaxy samples need HODs with 8 to 9 parameters for a reasonable description of galaxy clustering (Avila et al. 2020;Gonzalez-Perezetal.2020). Other empirical methods, such as age matching or conditional abundance matching (Hearin & Watson 2013; Hearin et al. 2014; Kulier & Ostriker 2015;Favoleetal.2021), have been able to reproduce the galaxy clustering using the additional correlation between SFR and secondary halo properties, such as halo age, but with some limitations. In this section, we implement a different approach to create starforming mock galaxies inspired by the fact that newly accreted mass provides gas to fuel the formation of new stars in galaxies and by the modelling carried out by EMERGE. For central galaxies, we assume that the SFR is proportional to the dark matter mass accretion rate, Mh t , multiplied by an efficiency function. We assume that the star formation efficiency is maximal at a halo mass M1(typically, ∼1012h−1M, near to the halo mass above which AGN feedback is important) and then it decays as a power law for both, more and less massive host haloes. The SFR of the central galaxies can be expressed as: SRFcentrals ∝Mh t 1 (Mh/M1)−β+(Mh/M1)γ.(3) Inspired by the standard abundance matching, we also test adding an additional scatter between the host halo mass and the SFR, finding no improvement in our predictions. For satellite galaxies, we assume the star formation rate is constant over a given time-scale and equal to that at the time when it reached its peak mass. For this, we stored the value of Mh twhen Msh =Mpeak while running the N-body simulations. We consider a given satellite galaxy quenched when tinfall, the time since accretion by a larger halo, is larger than a given quenching time-scale: tinfall >tquench ≡τ0×tdyn ×Mhost halo/1012h−1MτS,(4) with tdyn the dynamical time of the halo (0.1/H0), Mhost halo the host halo mass of the subhalo, and τ0&τSfree parameters. This last equation is slightly different from that used by EMERGE. There, the authors assume a direct relation with the stellar mass of galaxies, whereas here we assume a relationship with their host halo mass. Despite the large scatter between the host halo mass and the star formation rate of galaxies, we find that our specific implementation can better describe the SFR of TNG galaxies. Also, not using the stellar mass in the SFR calculation gives us more freedom when creating mocks with only SFR and not stellar mass. To avoid including additional free parameters in the model, we will omit any normalization parameter. This means that the value MNRAS 508, 175–189 (2021) Downloaded from https://academic.oup.com/mnras/article/508/1/175/6368853 by guest on 04 November 2025 A flexible SHAM model 181 Figure 4. Similar to Fig. 3but for the SFR parameters: the slopes in the star formation efficiency (β&γ), the halo mass where the peak of the star formation efficiency occurs (M1), the parameters that control deabundance of satellite galaxies τ0, and its dependence with the host halo mass τS. obtained for the SFR cannot be compared with other values in the literature, but that a galaxy sample with a fixed number density should include the same galaxy population. If the actual SFR value for each galaxy is needed, an abundance matching without scatter can be done to the full sample forcing it to mimic a specific SFR function. We show the impact of varying each of the SFR parameters in Fig. 4. Similar to Fig. 3, it shows the multipoles and the projected correlation function of galaxies at z=0 with a number density of 0.01 h3Mpc−3in our TNG-300-mimic simulation. Coloured lines in each row show the slopes in the star formation efficiency (β&γ), the halo mass of the peak of the star formation efficiency (M1), the parameters that control deabundance of satellite galaxies τ0, and its dependence with it host halo mass τS. The dependence of the galaxy clustering of βand γis weak, especially for β. Lower values of βand higher values of γwould mean selecting galaxies in more massive haloes. Since these haloes are more biased, the clustering of these galaxies should also be higher, as shown in the plot. The same occurs with M1,where a higher value of M1produces a more clustered galaxy sample, consistent with selecting more massive haloes. The dependence of the galaxy clustering on τ0and τSis strong, with higher values of these parameters implies higher amplitude galaxy clustering. This is expected since higher values of τ0include more satellite galaxies, which has a strong effect on the clustering of galaxies, especially at small scales, and higher values of τSwill preferentially select galaxies in more massive haloes (that have a higher bias). In the next section, we will show the results regarding galaxy clustering for stellar mass and SFR selection criteria. 6 REDSHIFT-SPACE GALAXY CLUSTERING PREDICTIONS In this section, we will compute the clustering predictions of our mocks and compare them to those in the TNG300. Since the simulation used to build the mocks has the same initial conditions and volume as the TNG300, the cosmic variance should contribute little to the differences in clustering. To find the best set of values for the free parameters in our extended SHAM, we minimize the chi-squared (χ2) between the predicted projected correlation function (wp), monopole (ξ=0), quadrupole (ξ=2), and hexadecapole (ξ=4) of the two-point correlation functions. We employ these statistics from 0.3 to 20 h−1Mpc. We compute uncertainties using 1000 jackknife samples, following the procedure of Zehavi et al. (2002) & Norberg et al. (2009). For simplicity, we assume no correlation between different scales (i.e. we only use the elements of the diagonal of the covariance matrix). To minimize χ2, we use a particle swarm optimization (also known as PSO; Kennedy & Eberhart 1995). In this technique, a group of particles (known as the swarm) move with random velocities inside the hyperspace we would like to explore. In each step, we compute the value of χ2of the particular parameter combination. The velocities of the particles are then updated so they can point towards the location where each of the particles found its minimum χ2and to the position with the minimum χ2found by all the particle as a swarm. The amount of velocity that remains constant after each step, how strong is the deviation to the best local and global value, are the parameters that control how much ‘exploration’ or ‘exploitation’ the PSO will have. We chose a more ‘explorative’ setup using 16 particles with MNRAS 508, 175–189 (2021) Downloaded from https://academic.oup.com/mnras/article/508/1/175/6368853 by guest on 04 November 2025 182 S. Contreras, R. E. Angulo and M. Zennaro Figure 5. The projected correlation function, wp(rp); the redshift space monopole, ξ=0; quadrupole, ξ=2; and the hexadecapole, ξ=4, for mock galaxies with a number density of n=0.01, 0.00316, & 0.001 h3Mpc−3selected according to their stellar mass in the TNG300 (black solid line), a basic SHAM (green dashed line), and our extended SHAM (red dashed line). As in Fig. 2, the two lowest number density samples have been displaced in the y-axis for clarity, with the lowest number density at the top. The vertical dotted line indicates the smallest scale, 0.3h−1Mpc, included in the fit of the free parameters of the extended SHAM model. Table 1. The best-fitting parameters of the extended SHAM for galaxies selected by stellar mass from the TNG300 (top part), the SDSS (middle part), and for galaxies selected by SFR of the TNG300 (bottom part). The fit is done for three different number densities and property cuts (as labelled) for the TNG300 and all samples simultaneously for the SDSS. Stellar mass TNG300 n/h3Mpc−3cut/h−1MσlogM tmerger fsfk,cfk,s 0.01 6.5 ×1090.134 0.636 0.034 0.377 −0.328 0.00316 2.2 ×1010 0.193 0.406 0.003 0.227 -0.127 0.001 4.5 ×1010 0.059 0.765 0.081 0.328 −0.032 Stellar mass SDSS σlogM tmerger fsfk,cfk,s 0.043 1.112 0.211 0.272 0.088 SFR TNG300 n/h3Mpc−3cut/Myr−1βγlog(M1)τ0τS 0.01 0.47 2.490 5.127 12.770 4.931 −0.363 0.00316 1.49 2.531 4.687 12.051, 3.869 -0.265 0.001 3.03 8.166 8.907 12.063 4.346 −0.307 400 steps each. The PSO typically converges in less than half of these steps. 6.1 Stellar mass predictions In Fig. 5, we show the clustering prediction of the TNG300 (black solid line) and our extended SHAM (red dashed lines) for three number densities, n=0.01 h3Mpc−3,n=0.00316 h3Mpc−3,&n= 0.001 h3Mpc−3based on the stellar mass of the galaxies. The shaded black region represents 1σuncertainty estimated from the jackknifes samples. For comparison, we show the clustering prediction of a standard SHAM with a scatter equally to that measured for TNG galaxies (σlog M=0.125). The best-fitting parameters for our extended SHAM are shown in the upper section of Table 1. The two highest number densities have similar values for their parameters, with some small differences compared to the lowest number density. Note, however, that since these correspond to physically different galaxies, there is no a priori reason to expect identical values for the SHAM parameters. Note that the clustering of a galaxy population created with these parameters may vary with cosmology or resolution, mostly because of the treatment of orphan galaxies. We can see that our extended SHAM model can successfully reproduce the clustering of the TNG300 hydrodynamic simulation for all number densities. This is remarkable considering that the TNG300 was carried out with a full physics model and with a resolution 64 times higher than our DMO counterpart. We also notice that the extended SHAM is more accurate compared to the standard SHAM, even for the highest number density. This is in agreement with previous results in the literature that indicate that orphan galaxies are necessary for SHAM to yield reliable results, even in low number densities (e.g. Guo & White 2014). We further notice that the clustering is well reproduced on all scales, even below 0.3 h−1Mpc, the minimum scale used in the fitting of the SHAM parameters (indicated as a vertical dotted line in the figure), suggesting the robustness of our model. We also notice that the clustering agrees better for the densest sample. Originally, we expected the opposite since at these number densities the clustering depends more on the additional parameters we include (such as the fraction of orphans). On the other hand, the statistical uncertainties of the clustering measurements are smaller for the denser samples, which allows the free parameter of the model to be better constrained and more accurate comparison. To study the precision of the mock in more detail, we show in Fig. 6 the halo occupation distribution (HOD), the galaxy assembly bias (GAB), and the kNN-CDF for the TNG300, our extended SHAM, and the standard SHAM for the same three galaxy samples previously considered. We emphasize that our model was not calibrated to reproduce these statistics and thus its performance is an independent test of the realism of the model. We find that our extended SHAM model is almost indistinguishable from the TNG300 regarding the mean number of galaxies (N) that populate the haloes of mass ‘M’ in the TNG300, for all number densities we explored. We find small differences at the transition between zero to one galaxy per halo, this is nevertheless not surprising since galaxy clustering poorly constrains σlog M,4the parameter that control this transition (see Zehavi et al. 2011). Note that to increase the accuracy of our comparison, we follow Contreras 4Notice that, while the σlog Mparameter in HOD and SHAM are related, they are two different concepts. MNRAS 508, 175–189 (2021) Downloaded from https://academic.oup.com/mnras/article/508/1/175/6368853 by guest on 04 November 2025 A flexible SHAM model 183 Figure 6. Top left-hand panel: The halo occupation distribution (HOD) for galaxies from the TNG300 (blue solid line), a basic SHAM technique (orange dotted line), and our extended SHAM model (green dashed line) for number densities of 0.01, 0.00316, and 0.001 h3Mpc−3basedontheir stellar mass, with the highest number density at the right of the plot. Top right-hand panel: The galaxy assembly bias signal, measured as the ratio of the 2PCF of the selected galaxy sample to that of its shuffled counterpart, following Croton et al. (2007). As in the top panel, solid black, dotted green, and dashed red lines show the predictions for the TNG300, basic SHAM and our extended SHAM. Bottom: The upper panel shows the k-nearest neighbour Cumulative Distribution Functions (kNN-CDF) for the TNG300 (solid lines) basic SHAM (dotted lines) and our extended SHAM model (dashed lines) for N=0, 1, 2, & 3 (red, green, blue, and magenta line, respectively). The lower part of this row shows the ratio of the kNN-CDF of the standard and extended SHAM to one of the TNG300. Overall, the agreement between galaxies in the TNG simulation and our extended SHAM is remarkable. et al. (2015) and first match the halo masses in both simulations, TNG300 and TNG300-mimic. The magnitude of GAB is estimated by computing the ratio between our galaxy sample and that where the position of galaxies was randomly shuffled among haloes of similar mass, following the procedure of Croton et al. (2007). This shuffling method eliminates any dependence of the galaxy population on any halo property other than mass (i.e. any assembly bias). We remind the reader that, although model has two distinct assembly bias parameters for central and satellite galaxies, we fit them simultaneously using the full galaxy clustering measurement. In the top righthand panel of Fig. 6we can see that the assembly bias level in our SHAM catalogues agrees remarkably well with that in the TNG300 simulations, with differences being below the per cent level. This confirms that our assembly bias modelling is actually correctly capturing the environmental dependences absent in SHAM. In a future work, we plan to use apply our model to directly constrain the level of assembly bias in observational clustering measurements. Finally, we test the kNN-CDF of the three number densities for k=0,1,2,& 3. The kNN-CDF (Banerjee & Abel 2021) measures the empirical cumulative distribution function of distances from a set of volume-filling Poisson distributed random points to the k-nearest data points. As discussed by (Banerjee & Abel 2021), this statistic mostly captures non-Gaussian features and higher-order correlations in a given spatial dataset. The performance of our extended SHAM shows to be highly accurate in all cases, with differences with respect to TNG samples typically below 5 per cent. In contrast, the standard SHAM displays notorious and systematic differences. It is surprising that, even at the lowest number density (where it is normally assumed that the standard SHAM has a good performance and where the clustering and HOD shows a good agreement), the predictions of the kNN-CDF are systematically below the TNG data. This shows the importance of the extensions of SHAM studied here, even in MNRAS 508, 175–189 (2021) Downloaded from https://academic.oup.com/mnras/article/508/1/175/6368853 by guest on 04 November 2025