Modelling the tightest relation between galaxy properties and dark matter halo properties from hydrodynamical simulations of galaxy formation

TL;DR

This study identifies the stellar mass at peak circular velocity as the best predictor of galaxy properties, demonstrating that simple models like SHAM can effectively predict galaxy clustering despite inherent scatter.

Contribution

It reveals that stellar mass at peak circular velocity provides the tightest galaxy-halo relation across all galaxy types, including satellites, and shows the minimal impact of scatter on clustering predictions.

Findings

Stellar mass at peak circular velocity yields the tightest galaxy-halo correlation.

Evolution correction reduces scatter in the galaxy-halo relation.

Simple SHAM models can accurately predict galaxy clustering despite scatter.

Abstract

We investigate how a property of a galaxy correlates most tightly with a property of its host dark matter halo, using state-of-the-art hydrodynamical simulations of galaxy formation EAGLE, Illustris, and IllustrisTNG. Unlike most of the previous work, our analyses focus on all types of galaxies, including both central and satellite galaxies. We find that the stellar mass of a galaxy at the epoch of the peak circular velocity with an evolution correction gives the tightest such correlation to the peak circular velocity $V_{\rm peak}$ of the galaxy's underling dark matter halo. The evolution of galaxy stellar mass reduces rather than increases scatter in such a relation. We also find that one major source of scatter comes from star stripping due to the strong interactions between galaxies. Even though, we show that the size of scatter predicted by hydrodynamical simulations has a negligible impact on the clustering of dense $V_{\rm peak}$-selected subhalo from simulations, which suggests that even the simplest subhalo abundance matching (SHAM), without scatter and any additional free parameter, can provide a robust prediction of galaxy clustering that can agree impressively well with the observations from the SDSS main galaxy survey.

Equations

$${\log }_{10}(M_{\mathrm{star}}[M_{\odot }{h}^{-2}]){|}_{z=0}={\log }_{10}(M_{\mathrm{star}}[M_{\odot }{h}^{-2}])(z)+\alpha \times z,$$

$$\xi (r_{\sigma },r_{\pi })=\frac{DD(r_{\sigma },r_{\pi })-2DR(r_{\sigma },r_{\pi })+RR(r_{\sigma },r_{\pi })}{RR(r_{\sigma },r_{\pi })},$$

$${\xi }_l(s)=\frac{2l+1}{2}{\int }_{-1}^{1}d\mu \xi (s,\mu )P_l(\mu ),$$

Full text

††pubyear: 2019††pagerange: Modelling the tightest relation between galaxy properties and dark matter halo properties from hydrodynamical simulations of galaxy formation–Modelling the tightest relation between galaxy properties and dark matter halo properties from hydrodynamical simulations of galaxy formation

1 Introduction

In modern analysis of the large-scale structure of galaxies, galaxies are usually treated as discrete points. The over-density field of galaxies is a Poisson sample of the underlying dark matter field Peebles (1980) and the statistics of galaxies does not depend on properties of galaxies. However, this assumption is indeed over-simplified Casas-Miranda et al. (2002). Galaxy formation is far more complicated than just a discrete point process. Galaxies are not formed in an isolated environment but rather interact frequently with one another. Not only can such interactions change the probability distribution of galaxies but also the properties of galaxy themselves. As a result, galaxies are biased tracers of the dark matter field. The way they trace the dark matter field depends both on their properties as well as their assembly histories. For instance, it is well known that red galaxies are more clustered than the blue ones; more luminous galaxies are more clustered than the fainter ones as well Zehavi et al. (2011). Even for galaxies with similar properties, such as those selected in the BOSS CMASS sample Dawson et al. (2012), where the galaxies have a fairly uniform mass distribution that peaks at around $\log(M/M_{\odot})\sim 11.3$, galaxy clustering is still dependent on the stellar-mass assembly history Montero-Dorta et al. (2017).

Since in the realistic case galaxies are not good tracers of the dark matter field, the key question then becomes how to accurately quantify the relation between them. The first approach for this is to do modeling. In $\Lambda$CDM, there has been significant progresses in modeling such biases over the past decade Scoccimarro (2004); Taruya et al. (2010); de la Torre et al. (2017); de la Torre & Guzzo (2012); Bianchi et al. (2012); Bianchi et al. (2015). The galaxy bias model now can yield reasonable accuracy in the quasi-linear regime. However, despite the progress, this approach is still unsatisfactory in several aspects. First, the accuracy and utility of these models have to be tested against mock galaxy surveys, as the galaxy bias models are essentially phenomenological and empirical. However, the mock catalogues are usually built for a particular survey or specific to a particular selection of galaxies. It is unclear how accurate these bias models can be applied to samples with a wide range of properties and complex biases. Moreover, it also remains difficult to assess the extent to which the underlining simplified assumptions and intrinsic limitations in the mock catalogues affect the calibration of these models in the first place. In addition, aside from the accuracy and utility issue, another concern is that current bias models are tested against mock galaxy catalogues based on simulations only in $\Lambda$CDM. It is unknown whether or not these models can still work in a modified gravity model. This is an important test because, unlike in $\Lambda$CDM, the relationship between galaxies and the dark matter field is much more complicated in modified gravity (MG) models. Galaxies in MG models are not even directly related to the dark matter field, as the formation, clustering, and motion of galaxies are dictated by a different potential $\phi_{+}=(\psi+\phi)/2$ rather than the lensing potential $\phi_{-}=(\psi-\phi)/2$ that is directly related to the true dark matter field. The relationship between $\phi_{+}=(\psi+\phi)/2$ and $\phi_{-}=(\psi-\phi)/2$ can be very complex He et al. (2015) in a modified gravity model.

The second approach makes use of N-body simulations, coupled to some phenomenological frameworks, such as the halo occupation distribution (HOD) Jing et al. (1998); Peacock & Smith (2000); Berlind & Weinberg (2002); Zheng et al. (2005, 2007); Leauthaud et al. (2011); Zehavi et al. (2011); Guo et al. (2015, 2016) or the conditional luminosity function (CLF) Yang et al. (2003, 2009), to link galaxies to dark matter halos. The basic idea of this approach is that galaxies reside in dark matter halos, namely, the densest regions of the underlying dark matter field. And the probability of the distribution of galaxies is only dependent on the masses of dark matter halos. Although these assumptions are obviously over-simplified which neglects some important effects such as the assembly bias Gao et al. (2005); Gao & White (2007), the HOD/CLF modeling turns out to be very successful in reproducing galaxy clustering even at very small scales for galaxy samples with a wide variety of different properties and complex biases Zehavi et al. (2011); Zheng et al. (2007); Guo et al. (2015, 2016). However, despite this, when applied to testing the underlining cosmological models, these frameworks have obvious limitations: they are too flexible and, in general, lack of strong physical motivations. In a modified gravity model, for example, even naively using the same HOD parameterisations as in $\Lambda$CDM, the framework can still be tuned to reproduce desired galaxy clustering Hernández-Aguayo et al. (2018), which is clearly unfeasible given that the dark matter halo properties and the processes of galaxy formation in a modified gravity model ought to be very different from those in $\Lambda$CDM.

In this work, rather than adhering to the doctrine that galaxies are tracers of the dark matter field, we regard galaxies as tracers of the dark matter halos: galaxies illuminate the dark matter halos and their properties are tracers of the properties of the dark matter halos. This idea is in line with the philosophy of subhalo abundance matching Kravtsov et al. (2004); Vale & Ostriker (2004); Conroy et al. (2006); Moster et al. (2010); Guo et al. (2010); Behroozi et al. (2013); Reddick et al. (2014). In this paper, we will further explore how a property of galaxy physically correlates to a property of its host dark matter halo. In order to achieve this, we use hydro-dynamical simulations of galaxy formation. In contrast to the sim-analytical galaxy formation models that are based on the DMO ones Cole et al. (2000); Guo et al. (2011), hydrodynamical simulations can trace the complicated interactions between baryons and dark matter in a self-consistent way based on the first principles of gravity and hydrodynamics. This is of paramount importance as the motion and clustering of galaxies, along with their host dark matter halos, are primarily dictated by gravity and hydrodynamics. In order to further strength our results and prevent potential biases due to the choice of a particular simulation, in this work we adopt three different simulations EAGLE Schaye et al. (2015); Crain et al. (2015); McAlpine et al. (2016), Illustris Nelson et al. (2015); Vogelsberger et al. (2014) and IllustrisTNG  Nelson et al. (2018a); Pillepich et al. (2018); Springel et al. (2018); Nelson et al. (2018b); Naiman et al. (2018); Marinacci et al. (2018). These simulations are different in many aspects such as the numerical methods used, the model of subgrid baryonic physics and, most importantly, the properties of simulated galaxies.

In addition, in contrast to most of the previous analysis that only focuses on central galaxies (e.g.Refs. Matthee et al. (2017); Matthee & Schaye (2019)), our analysis includes all types of galaxies, especially including a significant fraction of satellite galaxies. Unlike the central galaxies, satellite subhalos can not be easily matched from the DMO simulations to the full baryonic physics hydrodynamical ones by identifying their dark matter particles. It can be as high as 30$\%$ low-mass satellites in the DMO simulations that can not find their counterparts in the hydrodynamical ones Chaves-Montero et al. (2016). We, therefore, present a new approach to analyse the simulation data, which is different from the method used in  Chaves-Montero et al. (2016).

This paper is organised as follows: in Section 2, we introduce the hydrodynamcial simulations and the galaxy catalogue used. In Section 3, we discuss the physical relationship between a property of galaxy and a property of its host dark matter halo before major disruptions; in Section 4, we discuss the impact of disruptions on the galaxy property-halo property relation; in Section 5, we discuss how to efficiently model galaxy clustering using subhalos from simulations. in Section 6, we summarise and conclude this work.

2 Hydrodynamical simulations

The first simulation we adopt is the EAGLE simulation Schaye et al. (2015); Crain et al. (2015); McAlpine et al. (2016). Among the EAGLE suite, we use the largest simulation with a box size of $67.77{\rm Mpc}/h$ along one side. The mass resolution of gas particle in this simulation is $1.81\times 10^{6}M_{\odot}$ and a DM particle is $9.70\times 10^{6}M_{\odot}$. We only use well resolved galaxies with $M_{\rm star}>8\times 10^{7}M_{\odot}$. For the Illustris Nelson et al. (2015); Vogelsberger et al. (2014) simulation, we use the highest resolution run with a box size of $75{\rm Mpc}/h$ along one side. A gas particle in this simulation is $1.3\times 10^{6}M_{\odot}$ and a DM particle is $6.3\times 10^{6}M_{\odot}$. We focus on galaxies with $M_{\rm star}>8\times 10^{7}M_{\odot}$. In the IllustrisTNG suite Nelson et al. (2018a); Pillepich et al. (2018); Springel et al. (2018); Nelson et al. (2018b); Naiman et al. (2018); Marinacci et al. (2018), we use the simulation with a box size of $75{\rm Mpc}/h$ along one side, the same as Illustris. A gas particle in this simulation is $1.4\times 10^{6}M_{\odot}$ and a DM particle is $7.5\times 10^{6}M_{\odot}$. Again, we focus only on galaxies with $M_{\rm star}>8\times 10^{7}M_{\odot}$. Further, the stellar (or gas) mass used in this work is the total mass of stars (or gases) that are bounded to a galaxy.

3 Pre-disruption

The first galaxy property we should look at is the total baryonic mass $M_{\rm star+\rm gas}$ of a galaxy, namely, the total masses of stars and cold gas in the galaxy. This is motivated by the recent discovery of the mass discrepancy relation McGaugh et al. (2016), which finds that there is a tight empirical relation between the radial dependence of the enclosed baryonic-to-dynamical mass ratio and the baryonic acceleration. Such a relation also indicates that the total baryonic mass of a galaxy should be tightly correlated with the circular velocity of its host dark matter halo. In this work, we investigate this issue using hydrodynamic simulations. In order to simplify our analysis, we adopt the following strategy: before discussing a galaxy’s property at the current time, we first examine the galaxy’s property at the epoch of the peak value of the maximum circular velocities during its merger history $v_{\rm peak}$. This is to use a galaxy’s property before major disruption, which can avoid any drastic changes in the properties of galaxies. We then discuss the impact of disruption on a galaxy’s property later on.

The upper panels of Fig. 1 show the total baryonic mass of a galaxy as a function of $v_{\rm peak}$ in the EAGLE (left), Illustris (middle) and IllustrisTNG (right) simulations at three representative redshifts $z=1.49,0.87,0$. Galaxies in Fig. 1 are selected when their maximum circular velocities peak at the three specified redshifts during their merger histories. Hence, the value of the total baryonic mass for each galaxy in Fig. 1 is at $v_{\rm peak}$. Further, note that the samples used here include both centrals and satellites. From the lower panels of Fig. 1, EAGLE and IllustrisTNG simulations show relatively big scatter in the total baryonic mass-$v_{\rm peak}$ relation, which is around $\sigma(\Delta\log_{10}M_{\rm star+\rm gas})\sim 0.3$. The scatter also depends on redshift. However, Illustris simulation shows a smaller scatter, which is around $\sigma(\Delta\log_{10}M_{\rm star+\rm gas})\sim 0.2$. That hydrodynamical simulations do not give a consistent result on the total baryonic mass-$v_{\rm peak}$ relation indicates that the behavior of the gas component in a galaxy is dependent on the details of baryonic physics implemented in the simulations.

Next, we turn to the stellar component of a galaxy. The upper panels of Fig. 2 show the stellar mass of a galaxy (the value at the epoch of $v_{\rm peak}$) as a function of $v_{\rm peak}$ for the EAGLE (left), Illustris (middle) and IllustrisTNG (right) simulations at three different redshifts $z=1.49,0.87,0$, the same as in the previous case. In contrast to the total baryonic mass, all three simulations show a very similar tight relation between stellar mass and $v_{\rm peak}$ at a fixed redshift. From the lower panels of Fig. 2, different simulations even at different redshifts also show a very similar size of scatter. However, the zero-points in the scaling relations are different. This is due to the fact that a galaxy’s stellar mass evolves with time. A longer time of star-formation can lead to more stellar mass formed in a galaxy and the zero-points of the scaling relation is, therefore, higher at low redshift and lower at high redshift. Clearly, the evolution of the stellar mass plays an important role here. If we neglect such evolution effect and simply take the stellar mass at the epoch of $v_{\rm peak}$, the scatter in the stellar mass-$v_{\rm peak}$ relation will be very large. Therefore, the evolution of galaxy stellar mass should be accounted for. However, before we proceed to present how to correct for this effect, we first look at the specific star forming rate (sSFR), namely, the star-forming rate per unit stellar mass of a galaxy.

Figure 3 shows the sSFR of a galaxy as a function of $v_{\rm peak}$ at three representative redshifts in the EAGLE simulation. Around the mean values (represented by solid lines and usually called the main sequence galaxies), the scatter of the sSFR is very big. However, the average value of the sSFR is nearly a constant at a fixed redshift, which is almost independent of $v_{\rm peak}$, except the most massive ones at redshift zero. Thus, one can model such evolution of stellar mass as

[TABLE]

where $\alpha$ is a constant, which is different between the EAGLE, Illustris and IllustrisTNG simulations ($\alpha_{\rm EAGLE}=0.324\,,\,\alpha_{\rm Illustris}=0.243,,\,\alpha_{\rm IllustrisTNG}=0.194$). This is because the stellar mass functions and the star forming rates in these simulations are different. The stellar mass-$v_{\rm peak}$ relation at different redshifts then can be corrected to redshift zero using the above equation.

Figure 4 shows the evolution corrected stellar mass-$v_{\rm peak}$ relation at redshift zero. We only use the EAGLE simulation here for illustrative purposes. In Fig. 4, we also show the results at three representative redshifts as control. After taking into account the effect of stellar mass evolution, the evolution corrected stellar mass of galaxies are tightly correlated with $v_{\rm peak}$ (right panel in Fig. 4). The intrinsic scatter of this relation is very small. Note that here we have used the galaxy’s stellar mass that would be at redshift zero if there were no disruption rather than the true current stellar mass of the galaxy.

In addition, it is worth noting that although the scatter in the sSFR-$v_{\rm peak}$ relation is very large (see Fig. 3), the stellar mass and $v_{\rm peak}$ end up with a very tight correlation. This is indeed not a surprise because the sSFR has big fluctuations during the accretion history of stellar mass Matthee & Schaye (2019), which leads to the big scatter in the sSFR-$v_{\rm peak}$ relation. However, in the stellar mass-$v_{\rm peak}$ relation, this effect can be largely smoothed out as the stellar mass is the temporal cumulation of the sSFR.

4 Post-disruption

Unlike the pre-disruption stellar mass-$v_{\rm peak}$ relation, after $v_{\rm peak}$ great complexity comes in due to the interactions between galaxies. In order to illustrate this, in Fig. 5, we show the ratio of the stellar mass (pink points) and the maximum circular velocity $v_{\rm max}$ (blue points) of a galaxy at the present day to those values at the epoch of $v_{\rm peak}$ as a function of the stellar mass for the EAGLE (left), Illustris (middle) and IllustrisTNG (right) simulations. As shown by the blue points, the values of $\log_{10}(V_{\rm max}/V_{\rm peak})$ for a significant fraction of galaxies are much less than zero, which means that the disruption of the dark matter component of the galaxies is prevalent (in most cases $v_{\rm max}$ at the present day is much smaller than the value at the epoch of $v_{\rm peak}$). In contrast, the stellar component of most galaxies (pink points), on the other hand, grows rather than decreases significantly from the epoch of $V_{\rm peak}\,$ to the present day. This happens even for the cases that the dark matter and gas components undergo a dramatic loss. This indicates that even when a galaxy becomes a satellite merging into a large system, galaxies in most cases can still continue forming stars for a long time. Only in very extreme cases, the stellar components can suffer dramatic losses along with their dark matter components due to disruptions (pink points with $\log_{10}[M_{\rm star}(z=0)/M_{\rm star}(V_{\rm peak})]\ll 0$).

Figure 6 shows the post-disruption total baryonic mass-$V_{\rm peak}$ relation, compared with the post-disruption stellar mass-$V_{\rm peak}$ relation for different hydro-dynamical simulations. The left panels are for the EAGLE simulation, the middle ones are for the Illustris simulation and the right ones are for the IllustrisTNG simulations. The lower panels of Figure 6 show scatters for the above relations. After the epoch of $V_{\rm peak}$, both stars and gases suffer disruption. However, since there is only gravitational interaction between stars and there are significant non-gravitational interactions between gases(e.g. energetic feed-backs), gas is more prone to disruption than that of stars. The stellar mass-$V_{\rm peak}$ relation, therefore, is much tighter than that of the total baryonic mass-$V_{\rm peak}$ relation. Further, all three hydro-dynamical simulations yield a consistent result.

Figure 7 shows various stellar mass-$V_{\rm peak}$ relations for the EAGLE (left), Illustris (middle) and IllustrisTNG simulations. We show the evolution corrected relation (blue), in which stellar mass is taken at the epoch of $V_{\rm peak}$ but evolution corrected to redshift zero using Eq. (1). In the same plot, we also present the current stellar mass (post-disruption)-$V_{\rm peak}$ relation (pink) and the current stellar mass-$V_{\rm peak}$ relation but without disrupted galaxies $\log_{10}[M_{\rm star}(z=0)/M_{\rm star}(V_{\rm peak})]>0$ (brown). The lower panels show scatters in the various relations. The typical scatter on stellar mass in the post-disruption stellar mass-$V_{\rm peak}$ relation is around $\sigma(\Delta\log_{10}M_{\rm star})\sim 0.3$ in the EAGLE and IllustrisTNG simulations but slightly smaller $\sigma(\Delta\log_{10}M_{\rm star})\sim 0.25$ in the Illustris simulation. As in our analysis we include both centrals and satellites, the scatters here are slightly bigger than the values reported in the previous work that only use centrals (e.g. Refs. Matthee et al. (2017); Wechsler & Tinker (2018)). However, if we do not take into account the disrupted galaxies and only consider the galaxies that gain stellar mass after $V_{\rm peak}$, namely, $\log_{10}[M_{\rm star}(z=0)/M_{\rm star}(V_{\rm peak})]>0$, the scatter can be significantly reduced blow $\sigma(\Delta\log_{10}M_{\rm star})<0.2$ for galaxies with $V_{\rm peak}>100[\rm km/s]$ (brown). In this case, the scatter is close to the ideal case of the evolution corrected value using stellar mass at the epoch of $V_{\rm peak}$ (comparing the blue and brown lines in the lower panels of Fig. 7 ). This indicates that after $V_{\rm peak}$ when a galaxy is accreted by a larger system, stellar mass striping is a major source of scatter to the relation between the current stellar mass and $V_{\rm peak}$. Therefore, the tightest relation between a galaxy property and a halo property can only be achieved in terms of quantities before this disruption.

Nevertheless, it is worth pointing out that even in the case of post-disruption, the current stellar mass and $V_{\rm peak}$ still exhibit a very tight correlation. In Fig. 8, we compare such a relation (solid lines) with the evolution corrected value (dashed lines). Unlike Fig. 7, in Fig. 8 we show the scatter on $\log_{10}V_{\rm peak}$. This has more practical meaning as in observations we can only measure a galaxy’s current stellar mass (post-disruption). Given the observed galaxy catalogue, we usually add scatter in the $V_{\rm peak}$-ranked subhalo catalogues in $N$-body simulations and then compare them with observations. From Fig. 8, even considering the post-disruption galaxies, the scatter on $\log_{10}V_{\rm peak}$ is still very small $\sigma(\Delta\log_{10}V_{\rm peak})<0.08$, only except the most massive cases, which, however, account for only a small fraction of the total galaxies.

5 modelling galaxy clustering

The existence of such a tight correlation between galaxy stellar mass and $V_{\rm peak}$ has an important implication on how to efficiently model galaxy clustering from N-body simulations. As demonstrated in the previous sections, even for hydrodynamical simulations with very different subgrid baryonic physics and numerical methods used, they still produce similar predictions on such relation (see Fig. 8). Therefore, the predicted galaxy stellar mass-$V_{\rm peak}$ relation should be robust unless there are serious limitations in our current understanding of galaxy formation. Given such a tight correlation, galaxy clustering should be able to be modeled by using $V_{\rm peak}$-selected subhalos, which corresponds to stellar mass selected galaxies in the real observations. In this section, we will test this point with observations.

Before comparing simulations with observations, we first test the impact of scatter in the stellar mass-$V_{\rm peak}$ relation on the clustering of subhalos. In our analyses, we use the Small MultiDark Planck simulation (SMDPL) Klypin et al. (2016), which uses $3840^{3}$ dark matter particles in a box of $400{\rm Mpc}/h$ along one side. We focus on the redshift space since the redshift space clustering can be directly measured in the real observations. We add scatter in the halo catalogue as follows. We first take subhalos from the original catalogue and then replace their values of $V_{\rm peak}$ by drawing a random number around the logarithm of $V_{\rm peak}$ ( $\log_{10}V_{\rm peak}$) within 1$\sigma(\Delta\log_{10}V_{\rm peak})=0.08$ scatter. Then we re-rank all the subhalos in the catalogue and generate 200 realisations for such catalogues. The choice of the value of scatter is motivated by Fig. 8, from which $\sigma(\Delta\log_{10}V_{\rm peak})=0.08$ is large enough for most galaxies only except the most massive ones. The most massive galaxies, however, only account for a small fraction of the total number of galaxies.

In order to measure the redshift-space two point correlation function $\xi(r_{\sigma},r_{\pi})$, we use the Landy and Szalay estimator Landy & Szalay (1993)

[TABLE]

where $DD$, $DR$ and $RR$ are the data-data, data-random, random-random pair counts and $r_{\sigma}$, $r_{\pi}$ are the separations of galaxy pairs perpendicular and parallel to the line-of-sight direction, respectively. $\xi(r_{\sigma},r_{\pi})$ can be expanded in terms of Legendre polynomials

[TABLE]

where $P_{l}(\mu)$ is the Legendre polynomial of order $l$, $s=\sqrt{r_{\sigma}^{2}+r_{\pi}^{2}}$ and $\mu=r_{\pi}/s$. We use logarithmic bins in $s$ ($d\log_{10}s=0.2$) and measure $s$ up to $s_{\rm max}=25{\rm Mpc}/h$. We numerically work out the integration in Eq. (3) with $d\mu=0.05$.

Figure 9 shows the impact of scatter on the predicted multipoles of redshift space distortions (RSD) of subhalos with a number density of $n=0.01[{\rm Mpc}/h]^{-3}$. The solid lines are for the results without scatter. The dashed lines are for the results with scatter. The left panel of Fig. 9 shows the RSD multipoles $\xi_{0,2,4}$ multiplied by the redshift space separation $s$, which makes the RSD multipoles more sensitive to changes at small scales. The right panel shows similar results but for RSD multipoles $\xi_{0,2,4}$ multiplied by $s^{2}$. The larger index of $s$ in this case can make the RSD multipoles more sensitive to changes at large scales. Overall, in both cases, the scatter has a very limited impact on the RSD multipoles of subhalos. This result indeed can be expected. Given the high number density of our samples, the subhalos that would be in and out of the catalogue due to scatter only account for a small fraction of the total samples. The clustering of subhalos, therefore, is stable against such scatter. Further, note that in the above analysis, we have adopted the parallel approximation to add the RSD effects for subhalos in simulations.

Figure 10 shows the predicted multipoles (monopole $\xi_{0}$, quadrupole $\xi_{2}$, and hexadecapole $\xi_{4}$) of RSDs (black solid lines) compared to the measurements from the SDSS data (symbols with error bars), which is volume-limited and complete in galaxy stellar mass with a number density of $n=0.01[{\rm Mpc}/h]^{-3}$ (see Ref. He et al. (2018) for details). The predictions from subhalos agree impressively well with observations even at relatively large scales $r\sim 20{\rm Mpc}/h$. Note that here we use a realistic mock catalogue by collating 8 replicas of the box and place the observer at the centre. This mock catalogue has the same survey mask as the real data. The subhalos are simply selected by ranking $V_{\rm peak}$ from the SMDPL simulation. So there is no free parameter in our mock catalogue. In order to demonstrate the robustness of our RSD measurements, we test three different estimators of stellar masses: a template-fit method originally adopted in the NYU catalogue with the SDSS model magnitudes (stars) Blanton et al. (2005), the same template-fit method but using SDSS Petrosian magnitudes (circles), and a single-color method (triangles) Yang et al. (2007).

6 Conclusions

In this work, we have investigated the galaxy property-halo property relation using hydrodynamical simulations of galaxy formation. Unlike the conventional phenomenological frameworks such as HOD or CLF, the advantage of hydrodynamical simulations is that it can provide a clear physical picture about how a property of a galaxy relates to a property of its host dark matter halo. Based on state-of-the-art hydrodynamical simulations, such as EAGLE, Illustris and IllustrisTNG, we study under what circumstance a galaxy’s property could correlate most tightly with a property of its host dark matter halo. In addition, in contrast to most of the previous work, our analyses include all types of galaxies, especially containing a significant fraction of satellites. Unlike the central galaxies, satellite subhalos can not be easily matched from the DMO simulations to the full baryonic physics hydrodynamical ones. Therefore, we adopt a novel approach that is different from the one used in Ref. Chaves-Montero et al. (2016) to analyse the simulation data. The main findings of our work are summarised as follows:

• •

Despite the presence of the observed mass discrepancy relation, hydrodynamical simulations do not show a consensus on the tight correlation between the total baryonic mass of a galaxy and the circular velocity of its host dark matter halo if considering all types of galaxies at the present time (post-disruption) (see Fig. 6). Within $\Lambda$CDM, this can be expected as the gas component of a galaxy is strongly affected by various feed-backs in the processes of galaxy formation in the first place. The gas component can also have non-gravitational interactions and, therefore, is more prone to disruptions than that of the dark matter and stellar components since the latter two only have gravitational interactions among them. Therefore, such a tight correlation can not be expected to exist in satellite galaxies .

• •

The stellar mass of a galaxy at the epoch of $V_{\rm peak}$ with an evolution correction correlates most tightly to $V_{\rm peak}$ of the galaxy (see Fig. 4).

• •

After accretion, star formation in most galaxies can last for a quite long time. Most galaxies, therefore, can still gain significant amounts of stellar mass (see Fig. 5). The stellar mass evolution reduces rather than increases the scatter in the current stellar mass-$V_{\rm peak}$ relation.

• •

Aside from the intrinsic scatter, star stripping is the main cause for the scatter in the current stellar mass (post-disruption)-$V_{\rm peak}$ relation (see Fig. 7).

• •

Even for the current stellar mass (post-disruption)-$V_{\rm peak}$ relation, hydrodynamical simulations still predict a very small scatter on $\log_{10}V_{\rm peak}$ (see Fig. 8), which means that the clustering of stellar mass selected galaxies in observations can be well modeled by $V_{\rm peak}$-selected subhalos in simulations.

• •

Since the scatter predicted by hydrodynamical simulations has a very limited impact on the clustering of dense $V_{\rm peak}$-selected subhalos (see Fig. 9), even simple subhalo abundance matching (without scatter and free parameter) can yield a robust prediction of galaxy clustering. We show that when compared with the SDSS main galaxy samples, the predictions of SHAM mock catalogue agree impressively well with the observations (see Fig. 10).

It is worth noting that, in the above conclusions, we have assumed that baryons have a limited impact on the positions and motions of subhalos. This assumption is based on the findings reported in Refs. Chaves-Montero et al. (2016); Hellwing et al. (2016); McAlpine et al. (2016). The authors there have compared the clustering of subhalos that can be well matched from DMO simulations to the hydrodynamical ones with the clustering of the simulated galaxies. They find that the effect of baryons is very small on scales greater than $>1{\rm Mpc}/h$. However, an important issue here is that Refs. Chaves-Montero et al. (2016); Hellwing et al. (2016) only focus on the matched subhalos. But in fact a significant fraction of subhalos in DMO runs can not find their counterparts in the hydrodynamical ones in the first place as dark matter substructures are prone to disruptions due to numerical errors. It is unclear yet how to gauge the effect of baryons for those un-matched subhalos. Moreover, as recently pointed out by Refs. van den Bosch & Ogiya (2018); van den Bosch (2017), even modern state-of-the-art DMO simulations might not reliably resolve the dark matter substructures due to the insufficient force resolution in simulations. So to exactly pin down the physical effect of baryons, rather than numerical artifacts, on the clustering and motions of dark matter substructures are indeed highly non-trivial. A detailed analysis is needed, which, however, is beyond the scope of this work and will be explored in our future work.

Acknowledgements

We acknowledge the Virgo Consortium for making their simulation data available. The eagle simula- tions were performed using the DiRAC-2 facility at Durham, managed by the ICC, and the PRACE facil- ity Curie based in France at TGCC, CEA, Bruy‘eres- le-Chˆatel. J.H.H. acknowledges support of Nanjing University and this work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure.