Galaxy assembly bias of central galaxies in the Illustris simulation

TL;DR

This study investigates how central galaxy properties relate to halo assembly history in the Illustris simulation, revealing that stellar mass correlates more strongly with peak halo circular velocity than with halo mass, and proposing a simple model to incorporate assembly bias into halo models.

Contribution

It demonstrates that galaxy stellar mass and star formation properties are primarily linked to halo peak circular velocity and formation time, and introduces a model to include assembly bias in clustering analyses.

Findings

Stellar mass correlates more tightly with peak halo circular velocity than with halo mass.

Galaxy properties related to star formation history are more strongly connected to halo formation time and concentration.

A simple model links galaxy and halo bias factors based on property correlations, aiding in incorporating assembly bias into models.

Abstract

Galaxy assembly bias, the correlation between galaxy properties and halo properties at fixed halo mass, could be an important ingredient in halo-based modelling of galaxy clustering. We investigate the central galaxy assembly bias by studying the relation between various galaxy and halo properties in the Illustris hydrodynamic galaxy formation simulation. Galaxy stellar mass $\Mstar$ is found to have a tighter correlation with peak maximum halo circular velocity $\Vp$ than with halo mass $\Mh$. Once the correlation with $\Vp$ is accounted for, $\Mstar$ has nearly no dependence on any other halo assembly variables. The correlations between galaxy properties related to star formation history and halo assembly properties also show a cleaner form as a function of $\Vp$ than as a function of $\Mh$, with the main correlation being with halo formation time and to a less extent halo concentration. Based on the galaxy-halo relation, we present a simple model to relate the bias factors of a central galaxy sample and the corresponding halo sample, both selected based on assembly-related properties. It is found that they are connected by the correlation coefficient of the galaxy and halo properties used to define the two samples, which provides a reasonable description for the samples in the simulation and suggests a simple prescription to incorporate galaxy assembly bias into the halo model. By applying the model to the local galaxy clustering measurements in Lin et al. (2016), we infer that the correlation between star formation history or specific star formation rate and halo formation time is consistent with being weak.

Equations

$$V_{\mathrm{peak}}=170{\left(\frac{M_{\mathrm{h}}}{10^{12}{h}^{-1}{\mathrm{M}}_{\odot }}\right)}^{1/3}\mathrm{km}{\mathrm{s}}^{-1}$$

$$\rho =\frac{⟨xy⟩-⟨x⟩⟨y⟩}{{\sigma }_x{\sigma }_y},$$

$$b_G={\iint }_{G}dxdyp(x,y)b(x,y),$$

$$⟨x{⟩}_{x>t{\sigma }_x}=\frac{{\int }_{t{\sigma }_x}^{+\infty }dx{\int }_{-\infty }^{+\infty }dyxp(x,y)}{{\int }_{t{\sigma }_x}^{+\infty }dx{\int }_{-\infty }^{+\infty }dyp(x,y)}=\frac{{\sigma }_x\exp (-t^2/2)}{{\int }_t^{+\infty }dv\exp (-v^2/2)}$$

$$⟨x{⟩}_{y>t{\sigma }_y}=\frac{{\int }_{t{\sigma }_y}^{+\infty }dy{\int }_{-\infty }^{+\infty }dxxp(x,y)}{{\int }_{t{\sigma }_y}^{+\infty }dy{\int }_{-\infty }^{+\infty }dxp(x,y)}=\frac{\rho {\sigma }_x\exp (-t^2/2)}{{\int }_t^{+\infty }dv\exp (-v^2/2)},$$

$$f={\int }_t^{+\infty }dv\frac{1}{\sqrt{2\pi }}\exp (-v^2/2).$$

$$b_h=⟨b{⟩}_{x>t{\sigma }_x}=k⟨x{⟩}_{x>t{\sigma }_x}+b_c$$

$$b_g=⟨b{⟩}_{y>t{\sigma }_y}=k⟨x{⟩}_{y>t{\sigma }_y}+b_c.$$

$${\delta }_g^b=\rho {\delta }_h^b.$$

$$b_g=b_c+\frac{\partial b}{\partial x}\left[x_c+{\sigma }_x\frac{\rho }{\sqrt{2\pi }}\exp (-t^2/2)/f\right],$$

$$b_h=⟨b{⟩}_{x>t{\sigma }_x}=k_x⟨x{⟩}_{x>t{\sigma }_x}+k_y⟨y{⟩}_{x>t{\sigma }_x}+b_c$$

$$b_g=⟨b{⟩}_{y>t{\sigma }_y}=k_x⟨x{⟩}_{y>t{\sigma }_y}+k_y⟨y{⟩}_{y>t{\sigma }_y}+b_c,$$

$$⟨y{⟩}_{x>t{\sigma }_x}=\frac{{\int }_{t{\sigma }_x}^{+\infty }dx{\int }_{-\infty }^{+\infty }dyyp(x,y)}{{\int }_{t{\sigma }_x}^{+\infty }dx{\int }_{-\infty }^{+\infty }dyp(x,y)}=\frac{\rho {\sigma }_y\exp (-t^2/2)}{{\int }_t^{+\infty }dv\exp (-v^2/2)}$$

$$⟨y{⟩}_{y>t{\sigma }_y}=\frac{{\int }_{t{\sigma }_y}^{+\infty }dy{\int }_{-\infty }^{+\infty }dxyp(x,y)}{{\int }_{t{\sigma }_y}^{+\infty }dy{\int }_{-\infty }^{+\infty }dxp(x,y)}=\frac{{\sigma }_y\exp (-t^2/2)}{{\int }_t^{+\infty }dv\exp (-v^2/2)}.$$

$${\delta }_g^b=\frac{\rho +\gamma }{1+\rho \gamma }{\delta }_h^b,$$

$$b_g=b_c+\frac{\partial b}{\partial x}\left[x_c+\Delta x_c+{\sigma }_x\frac{{\rho }_{\mathrm{eff}}}{\sqrt{2\pi }}\exp (-t^2/2)/f\right],$$

Full text

††pubyear: 2018††pagerange: Galaxy assembly bias of central galaxies in the Illustris simulation–B

1 Introduction

It has been well established that galaxies form in dark matter haloes (White & Rees, 1978). As the first step to study galaxy formation and clustering, halo formation and clustering, which is dominated by gravity, have been extensively studied with analytic models (e.g. Press & Schechter, 1974; Bardden et al., 1986; Mo & White, 1996; Sheth & Tormen, 1999) and cosmological $N$-body simulations (e.g. Springel, 2005; Prada et al., 2012). It has been found that halo clustering depends not only on halo mass but also on halo assembly history and environment (e.g. Gao et al., 2005; Gao & White, 2007; Paranjape et al., 2018; Xu & Zheng, 2018; Han et al., 2019; Ramakrishnan et al., 2019), or halo definition (e.g. Villarreal et al., 2017; Mansfield & Kravtsov, 2019). This is called halo assembly bias, whose nature is still under investigation (e.g. Dalal et al., 2008; Castorina & Sheth, 2013; Shi & Sheth, 2018).

If galaxy properties are affected by halo formation and assembly history, halo assembly bias would translate to galaxy assembly bias. Operationally, galaxy assembly bias can be defined as that at fixed halo mass, the statistical galaxy content shows dependence on other halo variables or galaxy properties show correlations with halo assembly history (which would cause a certain degree of halo assembly bias to be inherited by galaxies and to show up in galaxy clustering). The widely adopted halo model (e.g. Cooray & Sheth, 2002) of interpreting galaxy clustering, such as the halo occupation distribution (e.g., Berlind & Weinberg, 2002; Zheng et al., 2005) or conditional luminosity function (e.g., Yang et al., 2003), makes the implicit assumption of no galaxy assembly bias. Such methods have been successfully applied to galaxy clustering (e.g. Zehavi et al., 2005; Zheng et al., 2007; Xu et al., 2018). However, if assembly bias is significant, neglecting it in the model would lead to incorrect inference of galaxy-halo connections and introduce possible systematics in cosmological constraints (e.g. Zentner et al. 2014, 2019; but see also McEwen & Weinberg 2016; McCarthy, Zheng & Guo 2019). Conversely, observationally inferred galaxy assembly bias would help understand galaxy formation.

The existence and strength of galaxy assembly bias are still a matter far from settled, either in theory or in observation. Galaxy assembly bias has been investigated in hydrodynamic (e.g. Berlind et al., 2003; Mehta, 2014; Chaves-Montero et al., 2016; Artale et al., 2018; Martizzi et al., 2019; Bose et al., 2019) or semi-analytic galaxy formation models (e.g. Croton et al., 2007; Contreras et al., 2019; Zehavi et al., 2018, 2019), focusing on the effect on galaxy occupation function and galaxy clustering. The results seem to depend on the implementation details of star formation and feedback. Studying galaxy assembly bias from observation has the difficulty of determining halo mass, and the results are not conclusive (e.g. Yang et al., 2006; Berlind et al., 2006; Lin et al., 2016; Zu et al., 2017; Guo et al., 2017). Given the potential importance of galaxy assembly bias in modelling galaxy clustering, in this paper we study the correlation between various central galaxy properties and halo properties in the Illustris hydrodynamic simulation (Vogelsberger et al., 2014a) at the halo level, aiming at providing useful insights in describing galaxy assembly bias.

The structure of the paper is as follows. In section  2, we introduce the simulation and the galaxy and halo catalogues. Then in section 3, we investigate the relation between galaxy and halo properties, with primary galaxy-halo properties in section  3.1 and general galaxy-halo properties in section  3.2. In section  3.3, we present a simple model to connect galaxy assembly bias with halo assembly bias. Finally, we summarise and discuss the results in section  4.

2 Simulation and galaxy-halo catalogue

In this work, we use galaxies and haloes from the state-of-the-art hydrodynamic galaxy formation simulation Illustris111http://www.illustris-project.org (Vogelsberger et al., 2014a; Nelson et al., 2015) to study galaxy assembly bias, which is able to produce different type of galaxies seen in observation (Vogelsberger et al., 2014b; Genel et al., 2014). In particular, we use the Illustris-2 simulation, which has a box size of 75$h^{-1}$Mpc on a side, and contains $910^{3}$ dark matter particles of mass $5\times 10^{7}{\rm M}_{\odot}$ and the same number of baryon particles of mass $1\times 10^{7}{\rm M}_{\odot}$. The mass resolution is sufficient for our purpose of studying (central) galaxies in haloes of more massive than a few times $10^{10}h^{-1}{\rm M}_{\odot}$. The simulation adopts a spatially-flat cosmology with the following parameters: $\Omega_{\rm m}=0.27$, $\Omega_{\rm b}=0.0456$, $h=0.70$, $n_{\rm s}=0.963$, and $\sigma_{8}=0.809$.

The haloes in the Illustris database are identified with the friends-of-friends (FoF) algorithm. As this algorithm is notorious for having the probability of bridging two separate halos into one halo, we apply the phase-space halo finder Rockstar (Behroozi et al., 2013a) to identify haloes, with dark matter particles extracted from 30 snapshots of the simulation (from z=3.94 to z=0). We then build the merger tree using the consistent-tree algorithm (Behroozi et al., 2013b).

We focus our study on the relation between central galaxies and dark matter haloes and the assembly effect at $z=0$. Galaxies in the Illustris-2 simulation are assigned to the Rockstar haloes. For each galaxy, if its distance $d_{gh}$ to the centre of a halo is smaller than the virial radius $R_{\rm vir}$ of the halo, it is assigned to this halo. In rare cases, a galaxy may be assigned to more than one haloes, as halos are not spherical. We then put the galaxy into the halo that corresponds to the lowest ratio of d_{gh}/$$R_{\rm vir}. With galaxies assigned to haloes, for each halo we define the most massive galaxy inside 0.2$R_{\rm vir}$ from the halo centre to be the central galaxy. If there is no galaxy inside 0.2$R_{\rm vir}$, we simply choose the most massive galaxy inside $R_{\rm vir}$ as the central galaxy.

With the catalogue of central galaxies and the associated host haloes, we study the relation between galaxy properties and halo properties and the assembly effects. The halo properties we focus on are:

• (1)

$M_{\rm h}$, halo mass enclosed in a volume with mean density of 200 times the background density of the universe;

• (2)

$V_{\rm peak}$, peak maximum circular velocity of the halo over its accretion history;

• (3)

$c$, halo concentration parameter, defined as the ratio of halo virial radius to scale radius;

• (4)

$a_{\rm M/2}$, cosmic scale factor when the halo obtains half of its current ($z=0$) total mass;

• (5)

$\dot{M}_{\rm h}$, halo mass accretion rate near $z=0$ (averaged between z=0 and z=0.197, about 2.4 Gyr, one dynamical time), in units of $h^{-1}{\rm M}_{\odot}\,{\rm yr}^{-1}$;

• (6)

$\dot{M}_{\rm h}/M_{\rm h}$, specific halo accretion rate, in units of ${\rm Gyr}^{-1}$.

The central galaxy properties we consider include:

• (1)

$M_{*}$, stellar mass (sum of masses of star particles within twice the stellar half mass radius);

• (2)

SFR, star formation rate within twice the stellar half mass radius;

• (3)

sSFR, specific star formation rate, the ratio of SFR to $M_{*}$;

• (4)

$g-r$, galaxy colour defined by the $g$-band and $r$-band luminosity.

3 Results

We aim at presenting the relation between central galaxies and haloes to learn about the correlation between halo formation and assembly and galaxy properties. In section 3.1 We first study how galaxy stellar mass depends on the primary halo properties ($M_{\rm h}$ and $V_{\rm peak}$). Then we investigate how various halo and galaxy properties are correlated in section 3.2. Finally, in section 3.3 we use a simplified model to describe the connection between galaxy and halo assembly bias factor.

3.1 Relationship between stellar mass and halo properties

Galaxy stellar mass is a primary property inferred from observation. The relation between this primary galaxy property and certain primary halo property (e.g., $M_{\rm h}$ and $V_{\rm peak}$) can be established based on subhalo abundance matching or modelling the stellar mass dependent clustering, which encodes information about galaxy formation. Here we show the relation predicted by the Illustris simulation and study how tight stellar mass correlates with $M_{\rm h}$ and $V_{\rm peak}$.

The top-left panel of Fig. 1 shows $M_{*}$ as function of $M_{\rm h}$, colour-coded with values of $V_{\rm peak}$. Galaxy stellar mass $M_{*}$ increases steeply with $M_{\rm h}$ at $\log[M_{*}/(h^{-1}{\rm M}_{\odot})]<12$ and then slowly at $\log[M_{*}/(h^{-1}{\rm M}_{\odot})]>12$, a trend similar to that inferred from observation (e.g. Behroozi et al., 2010; Leauthaud et al., 2012; Zu & Mandelbaum, 2015). The scatter in $M_{*}$ at fixed halo mass in the $M_{*}$–$M_{\rm h}$ relation decreases with increasing $M_{\rm h}$ (solid curve in Fig. 2), varying from about 0.3 dex at the low-mass end to about 0.17 dex at the high-mass end. The scatter at the high-mass end is consistent with the value $\sim 0.16$ dex inferred from galaxy clustering modelling (e.g. Tinker et al., 2017). One source of the scatter can be the halo formation history (Tinker, 2017), which may affect the growth history of stellar mass (either from star formation or galaxy merging; e.g. Gu et al. 2016).

The colour code in $V_{\rm peak}$ in the top-left panel enables us to see how the scatter in the $M_{*}$–$M_{\rm h}$ relation may be connected to halo assembly. On average, $V_{\rm peak}$ and $M_{\rm h}$ are correlated, and the mean relation is found to be well described by

[TABLE]

in the Illustris simulation. However, there is scatter on top of the mean relation, and at fixed $M_{\rm h}$, the distribution of $V_{\rm peak}$ reflects that in the assembly history. As can be seen in the top-left panel, the assembly of haloes encoded in $V_{\rm peak}$ does contribute to the scatter in the $M_{*}$–$M_{\rm h}$ relation – at fixed $M_{\rm h}$, galaxies residing in haloes of higher $V_{\rm peak}$ tend to have higher stellar mass, especially at the low mass end.

To see how well the scatter in $M_{*}$ can be attributed to the scatter in $V_{\rm peak}$, in the top-right panel of Fig. 1, we plot the $M_{*}$–$V_{\rm peak}$ relation. It follows a similar trend seen in the $M_{*}$–$M_{\rm h}$ relation, steeper (shallower) dependence of $M_{*}$ on $V_{\rm peak}$ at the low (high) $V_{\rm peak}$ end, which is expected given the correlation between $V_{\rm peak}$ and $M_{\rm h}$. The $M_{*}$–$V_{\rm peak}$ relation appears to be tighter than the $M_{*}$–$M_{\rm h}$ relation, in the sense that at fixed $V_{\rm peak}$ the scatter in $M_{*}$ is lower than that at the corresponding $M_{\rm h}$ (see Matthee et al. 2017 for a similar result in terms of $z=0$ maximum halo circular velocity with the EAGLE simulation). The scatter varies from $\sim$0.28 dex at low $V_{\rm peak}$ to $\sim$0.13 dex at high $V_{\rm peak}$ (dashed curve in Fig. 2). In the $M_{*}$–$V_{\rm peak}$ plot (top-right panel), the contours are colour-coded by $M_{\rm h}$. Unlike the $M_{\rm h}$ case in the top-left panel, we find that $M_{*}$ does not show a clear dependence on $M_{\rm h}$ – at fixed $V_{\rm peak}$ (i.e. by taking a vertical cut in the plot), $M_{*}$ in haloes of different $M_{\rm h}$ appears to follow similar mean and scatter.

In the bottom-left panel, we show the effect of the other halo assembly property $a_{\rm M/2}$ on the $M_{*}$–$M_{\rm h}$ relation. Each curve show the mean $M_{*}$–$M_{\rm h}$ relation for haloes in one $a_{\rm M/2}$ bin. The scatter around the mean relation is illustrated with the shaded region, only shown for the latest forming haloes to avoid crowdedness. There is a clear and substantial dependence of the $M_{*}$–$M_{\rm h}$ relation on the assembly property $a_{\rm M/2}$– at fixed $M_{\rm h}$, haloes forming earlier (smaller $a_{\rm M/2}$) tend to host galaxies of higher $M_{*}$. Such a trend is consistent with previous work based on the EAGLE simulation (Matthee et al., 2017) and semi-analytic galaxy formation model (Zehavi et al., 2018).

Switching to the $M_{*}$–$V_{\rm peak}$ relation (bottom-right panel), we find that at fixed $V_{\rm peak}$ there is no dependence of $M_{*}$ on the assembly property $a_{\rm M/2}$ , with the curves of different $a_{\rm M/2}$ bins all falling on top of each other. As further shown in section 3.2, $V_{\rm peak}$ is able to absorb the effect on $M_{*}$ from any other halo assembly variable. This is consistent with the results using $z=0$ maximum halo circular velocity (Matthee et al., 2017). We will discuss this remarkable result in section 4.

The other commonly adopted halo circular velocity quantity is $V_{\rm max}$, the present-day maximum circular velocity of a halo. In Appendix A, we show that the $M_{*}$–$V_{\rm max}$ relation is similar to the $M_{*}$–$V_{\rm peak}$ relation, but with the former showing a little bit residual dependence on halo formation time. In Appendix A, by replacing $V_{\rm peak}$ used here from the full physics simulation with the one from the dark-matter-only simulation, we also show that baryon effect on the halo potential plays almost no role in driving $V_{\rm peak}$ to track $M_{*}$.

The tighter correlation between $M_{*}$–$V_{\rm peak}$, in comparison to $M_{*}$–$M_{\rm h}$, suggests that in galaxy clustering modelling galaxy assembly bias effect can be partially account for by switching from a $M_{\rm h}$–based model to a $V_{\rm peak}$-based model. That is, the galaxy-halo relation is parameterised as a function of $V_{\rm peak}$. This is in line with the finding by Chaves-Montero et al. (2016) using the EAGLE simulation, while their velocity quantity most strongly correlating with $M_{*}$ is slightly different, ${\rm V_{relax}}$, the highest value of the maximum circular velocity of a subhalo with a relaxation criterion imposed. Chaves-Montero et al. (2016) show that ${\rm V_{relax}}$ is able to capture the majority (but not all) of the galaxy assembly bias effect on galaxy clustering for stellar-mass-based samples. The tighter correlation between $M_{*}$–$V_{\rm peak}$ may also be the reason that when using subhalo abundance matching or its variants to model galaxy clustering, halo circular velocity-based models usually have good performances (e.g. Reddick et al., 2013; Guo et al., 2016).

To investigate how well assembly bias effect in Illustris central galaxy clustering can be absorbed by $V_{\rm peak}$, we perform a test by constructing shuffled galaxy catalogues (e.g. Croton et al., 2007; Zu, Zheng, Zhu & Jing, 2008) and compare the two-point correlation functions of the original and shuffled central galaxy samples, selected by stellar mass thresholds. By dividing haloes into small $V_{\rm peak}$ bins, we randomly shuffle the values of $M_{*}$ among haloes in each $V_{\rm peak}$ bin, which eliminates the dependence of $M_{*}$ on halo properties other than $V_{\rm peak}$. If $V_{\rm peak}$ can account for all assembly effects on central galaxy $M_{*}$, we would expect the $V_{\rm peak}$-shuffled sample to have the same two-point correlation function as the original, unshuffled sample. For comparison, we also similarly construct $V_{\rm max}$-shuffled and $M_{\rm h}$-shuffled samples. In the left panel of Fig. 3, the result is shown for the $\log[M_{*}/(h^{-1}{\rm M}_{\odot})]$>9 sample. The $M_{\rm h}$-shuffled sample show the largest deviation in the two-point correlation function from the original sample, with the assembly effect on galaxy clustering at a level of 10 per cent. The level of deviation is small (below $\sim$2 per cent on most scales above 1$h^{-1}{\rm Mpc}$) for the $V_{\rm peak}$-shuffled sample, suggesting that $V_{\rm peak}$ is able to absorb the majority of the assembly bias on stellar mass. The result with the $V_{\rm max}$-shuffled sample is in between, and $V_{\rm max}$ is not as efficient as $V_{\rm peak}$ in absorbing the assembly bias effect on stellar mass (consistent with that seen in Appendix A). For the $\log[M_{*}/(h^{-1}{\rm M}_{\odot})]$>10 sample (right panel), the assembly bias effect is also at the 10 per cent level, and both $V_{\rm peak}$ and $V_{\rm max}$ are able to roughly absorb more than half of it, with $V_{\rm peak}$ behaving marginally better than $V_{\rm max}$. We note that the overall trend is similar to that in Zehavi et al. (2019) based on galaxies in a semi-analytic galaxy formation model (top panel of their fig.8 for central galaxies). As a whole, we see that $V_{\rm peak}$ is able to absorb a large fraction of the assembly effect on central galaxy stellar mass, more efficient for lower mass haloes, suggesting that $V_{\rm peak}$ may serve as a good halo variable to replace $M_{\rm h}$ in modelling clustering of galaxies selected by stellar mass.

In what follows, besides a further investigation of the dependence of stellar mass on halo assembly variables, we also extend the study to other galaxy properties.

3.2 Relationship between galaxy properties and halo properties

With the 7 halo properties ($M_{\rm h}$, $V_{\rm peak}$, $c$, $a_{\rm M/2}$, $\lambda$, $\dot{M}_{\rm h}$, and $\dot{M}_{\rm h}/M_{\rm h}$) and the 4 central galaxy properties ($M_{*}$, SFR, sSFR, and $g-r$ colour), we study the correlations between them. To aid the discussion, we also present the correlations among halo properties and those among galaxy properties.

3.2.1 At fixed $M_{\rm h}$

Fig. 4 shows the correlation between each pair of the halo and galaxy properties for central galaxies in haloes of a narrow mass bin, $\log[M_{\rm h}/(h^{-1}{\rm M}_{\odot})]=12.0\pm 0.1$. In each contour panel, the contours indicate the 68.3 and 95.4 per cent of the distribution of the pair of properties. The number labelled in each panel is the Pearson correlation coefficient $\rho$ of the two properties, indicating how strong the correlation is. It is calculated as

[TABLE]

where $x$ and $y$ denote the two properties, $\langle\rangle$ means average, and $\sigma_{x}$ and $\sigma_{y}$ are the standard deviations of $x$ and $y$. The panel at the top of each column shows the marginalised distribution of the property labelled at the $x$-axis of the column.

The panels with red contours (i.e. the top 6 rows and left 6 columns of contour panels) display the correlation between halo properties. Within the small but finite halo mass bin, all but one halo property shows almost no correlation with $M_{\rm h}$ (correlation coefficient close to zero). The exception is $V_{\rm peak}$, and the correlation is simply driven by the $V_{\rm peak}\propto M_{\rm h}^{1/3}$ mean relation. At fixed $M_{\rm h}$, any pair of halo properties show a substantial correlation (with $|\rho|$ above 0.2). The nearly perfect correlation ($\rho=0.983$) between $\dot{M}_{\rm h}$ and $\dot{M}_{\rm h}/M_{\rm h}$ is a consequence of fixed $M_{\rm h}$. Overall the correlation trend is that haloes of higher $V_{\rm peak}$ are more concentrated, form earlier, spin more slowly, and have lower accretion rate, which have been seen in previous work (e.g. Jeeson-Daniel et al., 2011; Han et al., 2019; Xu & Zheng, 2018).

The panels with black contours (i.e. the bottom 4 rows and the left 7 columns of contour panels) show the correlation between halo and galaxy properties. The correlation between $M_{*}$ and $M_{\rm h}$ shows up because of the finite size of the halo mass bin. At fixed $M_{\rm h}$, the central galaxy stellar mass $M_{*}$ correlates with all other halo properties – haloes of higher $V_{\rm peak}$, higher concentration, earlier formation, lower spin, and lower accretion rate tend to host more massive central galaxies. On the contrary, the SFR shows no strong correlation with any halo properties. The most significant one is with halo formation time ($\rho\sim 0.16$), with on average higher SFR in haloes of later formation. Given the substantial correlation between $M_{*}$ and halo properties and the weak or lack of correlation between SFR and halo properties, the sSFR ($\equiv{\rm SFR}/M_{*}$) is expected to correlate well with halo properties, but in an trend opposite to and weaker than that with $M_{*}$. This is indeed the case. The most significant correlation is with $V_{\rm peak}$ or $a_{\rm M/2}$ (both with $|\rho|\sim 0.4$). The correlation between sSFR and the average halo accretion rate over the past dynamic time is there but not strong ($\rho\sim 0.16$). The correlation between galaxy colour $g-r$ and halo properties essentially follows the case of sSFR (with a sign change in $\rho$; redder galaxies having lower sSFR).

The panels with blue contours (i.e. the right 3 columns of contour panels), the correlations between pairs of galaxy properties at fixed halo mass are shown. The SFR positively correlates with $M_{*}$ ($\rho\sim 0.54$), and the mean relation has a slope close to unity, ${\rm SFR}\propto M_{*}$. It resembles the star-forming main sequence (e.g. Brinchmann et al., 2004; Speagle et al., 2014; Santini et al., 2017). That is, even if we only consider central galaxies in haloes of fixed mass, the star-forming main sequence emerges. Note that this SFR–$M_{*}$ correlation is not driven by the correlation of SFR and $M_{*}$ with a common halo variable we consider here. In fact, from Fig. 4, it can be seen that their correlation with a common halo variable may lead to the opposite effect. For example, haloes of earlier formation tend to host central galaxies of higher $M_{*}$ and lower SFR, and naively this would imply an anti-correlation between SFR and $M_{*}$, opposite to what is found here. While it is possible that the halo-level star-forming main sequence is related to a halo variable not considered here, it is more likely that the sequence is driven by baryonic physics, which may have complicated dependence on or decouple from halo formation history. Unlike the SFR, the sSFR shows little dependence on $M_{*}$. However, the sSFR is tightly correlated with the SFR ($\rho\sim 0.88$). Given that sSFR$\equiv$SFR/$M_{*}$, the pattern in the mutual correlations among $M_{*}$, SFR, and sSFR can be achieved if the SFR–$M_{*}$ correlation coefficient is close to the ratio of the scatters in $\log M_{*}$ and $\log{\rm SFR}$222To see this, let $x=\log M_{*}$, $y=\log{\rm SFR}$, and $z=\log{\rm sSFR}=y-x$. We can derive the relation among the correlation coefficients, $\rho_{xz}/\rho_{yz}=(\rho_{xy}-\sigma_{x}/\sigma_{y})/(1-\rho_{xy}\sigma_{x}/\sigma_{y})$. For $\rho_{xz}$ to be near zero, we have $\rho_{xy}\sim\sigma_{x}/\sigma_{y}$. , which appears to be the case. Galaxy $g-r$ colour strongly correlates with sSFR and follows the same trends as sSFR in its correlations with $M_{*}$ and SFR.

Overall, at fixed halo mass, galaxy properties other than SFR show significant correlations with one or more halo properties, manifesting galaxy assembly bias at the level of haloes. The correlations among galaxy properties, however, may largely result from baryonic physics, given that the trend cannot be simply explained by their correlation with halo properties. In section 3.1, it is found that switching from $M_{\rm h}$ to $V_{\rm peak}$ can remove the dependence of $M_{*}$ on other halo properties. We now extend the investigation to other galaxy properties.

3.2.2 At fixed $V_{\rm peak}$

Fig. 5 is similar to Fig. 4, but the correlations are presented for haloes at a fixed $V_{\rm peak}$ bin, $\log[V_{\rm peak}/({\rm km\,s^{-1}})]=2.23\pm 0.03$. The correlations among halo properties (in panels with red contours) are similar to those in Fig. 4, and there are additional correlations between $M_{\rm h}$ and other halo properties.

The galaxy-halo correlations are shown in panels with black contours. The finite bin size in $V_{\rm peak}$ makes the $M_{*}$–$V_{\rm peak}$ correlation show up. Other than this (and the one with $M_{\rm h}$), $M_{*}$ does not correlate with any other halo properties at fixed $V_{\rm peak}$, reinforcing the result in section 3.1. It indicates that the correlations of $M_{*}$ with halo properties seen at fixed $M_{\rm h}$ (Fig. 4) can be attributed to the $M_{*}$–$V_{\rm peak}$ correlation and the correlation of $V_{\rm peak}$ with other halo properties. For the SFR, at fixed $V_{\rm peak}$, it correlates significantly with $M_{\rm h}$ and $a_{\rm M/2}$, higher SFR in haloes of higher mass and later formation. The correlations between SFR and other halo properties are weak. The correlations between sSFR (or colour) and halo properties closely follow the SFR case.

In terms of the galaxy-galaxy correlations, the trends are similar to those seen at fixed halo mass, except that the sSFR and colour now show clear correlations with $M_{*}$.

As a whole, using $V_{\rm peak}$ as the halo variable largely removes the correlations between $M_{*}$ and other halo assembly properties, and the dependences of SFR, sSFR, and colour on halo assembly variables follow each other.

3.2.3 Dependence on $M_{\rm h}$ and $V_{\rm peak}$

The correlations shown in Fig. 4 and Fig. 5 are for haloes of $\log[M_{\rm h}/(h^{-1}{\rm M}_{\odot})]\sim 12.0$ and $\log[V_{\rm peak}/({\rm km\,s^{-1}})]\sim 2.23$. To obtain a full picture, in Fig. 6 we present the $M_{\rm h}$ and $V_{\rm peak}$ dependent Pearson correlation coefficients for the various pairs of galaxy and halo properties, by performing the calculation in different $M_{\rm h}$ and $V_{\rm peak}$ bins, respectively. The panels correspond to those in Fig. 4 and Fig. 5, and the correlation shown in a panel of a given row and column is between the property as labelled at the far left of the row and that at the bottom of the column. In each panel, the solid (dashed) curve is the dependence of $\rho$ on $M_{\rm h}$ ($V_{\rm peak}$), with zero correlation marked by the black dotted curve. Note that only $V_{\rm peak}$ is shown on the $x$-axis and the corresponding $M_{\rm h}$ can be obtained according to $M_{\rm h}\propto V_{\rm peak}^{3}$ from equation (1).

The panels with red curves show the correlations among halo properties. If we limit to halo properties other than $M_{\rm h}$ and $V_{\rm peak}$, we find that the correlation of any pair of the assembly variables only weakly depends on $M_{\rm h}$ or $V_{\rm peak}$ if any (manifested by the nearly flat curves) and that the correlation strength does not depend on whether we use $M_{\rm h}$ or $V_{\rm peak}$ bins (manifested by the highly overlapped solid and dashed curves).

For the galaxy-halo correlations (in panels with black curves), in terms of $M_{\rm h}$ dependence, the strongest correlation between galaxy and halo property is found between $M_{*}$ and $V_{\rm peak}$/$c$/$a_{\rm M/2}$ in low mass halos, with $|\rho|\sim 0.5$–0.6. It holds true in the full range of haloes considered here that using $V_{\rm peak}$ largely removes the correlation between $M_{*}$ and any other halo assembly variable (dashed curves around zero). The only exception is that $M_{*}$ appears to be slightly anti-correlated with $c$ in low-$V_{\rm peak}$ haloes. With $V_{\rm peak}$ the dependences of correlation on $V_{\rm peak}$ for SFR, sSFR, and colour closely track each other, which is not the case for those on $M_{\rm h}$. With $V_{\rm peak}$ as the primary halo variable, star formation related properties (SFR, sSFR, and colour) mainly show dependence on halo formation time and then halo concentration.

For galaxy properties (in panels with blue curves), the correlation between SFR and $M_{*}$ reaches a maximum in haloes of $\sim 10^{12}h^{-1}{\rm M}_{\odot}$. It weakens in haloes of higher $M_{\rm h}$ or $V_{\rm peak}$, probably because galaxies move away from the star-forming main sequence and passive evolution starts to dominate. However, the tight sSFR–SFR correlation persists over the full $M_{\rm h}$ or $V_{\rm peak}$ range. For colour, the correlation with SFR and sSFR is weak in low-$M_{\rm h}$ or low-$V_{\rm peak}$ haloes and becomes stronger in haloes of higher $M_{\rm h}$ or $V_{\rm peak}$.

The results indicate that in the Illustris simulation galaxy formation ties to $V_{\rm peak}$ more closely than $M_{\rm h}$. In comparison with the $M_{\rm h}$-based results, we find that in terms of $V_{\rm peak}$, galaxy properties show a cleaner trend in the correlation with other halo assembly variable, such as the lack of correlation for $M_{*}$ and the similar correlation pattern for SFR/sSFR/colour. It suggests that $V_{\rm peak}$–based halo model would be a good choice for capturing galaxy assembly bias effect (e.g. with $M_{*}$-based galaxy samples) and for studying galaxy assembly bias (e.g. with SFR/sSFR/colour-based samples). With $V_{\rm peak}$ as the primary halo variable in the model, halo formation time and concentration would be the main options for the secondary variable to describe the relation between halos and star formation related quantities, with the former preferable.

3.3 Assembly bias of central galaxies

With the set of galaxy and halo properties investigated in section 3.2, we do not find a galaxy property that 100 per cent correlates with a halo assembly property. For the $M_{\rm h}$ dependence, the strongest correlation has $|\rho|\sim$0.5–0.6, between $M_{*}$ and $V_{\rm peak}$/$c$/$a_{\rm M/2}$. It means that halo assembly bias cannot be fully inherited by galaxies and be fully translated to galaxy assembly bias. Galaxy assembly bias should be different from halo assembly bias. For example, for haloes of the same mass, we can split them into two halo samples of low and high concentrations and then split central galaxies into two galaxy samples with low and high $M_{*}$. There would be a difference in the clustering of the two halo samples, as well as in that of the two galaxy samples. Given that $M_{*}$ is not perfectly correlated with $c$, we expect that the difference in the galaxy samples is smaller than that in the halo samples. As connecting galaxy assembly bias to halo assembly bias at the halo level can be an important ingredient in incorporating assembly bias effect into clustering model, we develop a simplified model below to understand the connection between galaxy and halo assembly bias.

Let us consider haloes at fixed $M_{\rm h}$ (or $V_{\rm peak}$) and focus on one halo assembly variable $x$ (e.g. $a_{\rm M/2}$ or concentration) and one galaxy property $y$ (e.g. $M_{*}$ or SFR). Without losing generality, $y$ is assumed to be positively correlated with $x$. The joint distribution of $x$ and $y$ is illustrated by an ellipse in the top panel of Fig. 7. We can form a halo sample by selecting the fraction $f$ of haloes with the highest $x$ (to the right of the vertical dashed line, the region inside the dashed curve) and a galaxy sample by selecting the same fraction of central galaxies with the highest $y$ (above the red-purple dividing line, the region in red). Then, what is the relation between the bias factors of the galaxy and halo sample? To answer the question, in addition to $p(x,y)$, we also need the dependence of the halo bias factor $b(x,y)$ on $x$ and $y$. For a sample $G$ defined by cuts in galaxy property $y$ or more generally cuts in the $x$-$y$ plane, the corresponding bias factor from haloes of fixed $M_{\rm h}$ or $V_{\rm peak}$ can be expressed as

[TABLE]

where $G$ indicates the region defined by the cuts.

To proceed, we make the following assumptions – (1) The joint distribution of galaxy and halo properties follows a 2-dimensional (2D) Gaussian $p(x,y)$, characterised by the centre $(x_{c},y_{c})$, standard deviations $\sigma_{x}$ and $\sigma_{y}$, and the correlation $\rho$ between $x$ and $y$. We can take $(x_{c},y_{c})=(0,0)$ by shifting $x$ and $y$. A 2D Gaussian function is a reasonable approximation for the distributions seen in Fig.4 and Fig. 5, which also follows the Taylor expansion of the logarithmic of the distribution function to the second order. (2) Galaxy property $y$ has a strong dependence on the halo assembly property $x$, and only weakly on other halo assembly variables. That is, any galaxy assembly bias effect is a result of inheriting from the assembly bias effect in halo property $x$, not other halo assembly properties. With this assumption, $b(x,y)$ has no additional dependence other than $x$, simplifying to $b(x,y)=b(x)$. (3) At the fixed $M_{\rm h}$ (or $V_{\rm peak}$), halo bias factor is linear with respect to halo property $x$, $b(x)=kx+b_{c}$, the first order approximation from Taylor expansion. The value $b_{c}$ is the bias factor at $x=0$, which is also the average halo bias factor for haloes at mass $M_{\rm h}$ (or $V_{\rm peak}$). The slope $k$ is the first derivative of $b$ with respect to $x$, $k=\partial b/\partial x$. In Appendix B, we present and test an extension of the model by loosing the assumptions (2) and (3), and the tests confirm that the simplification here is reasonable.

For the top $f$ fraction of haloes with the highest $x$ and that of galaxies with the highest $y$, the distribution of halo property $x$ is shown in the bottom panel of Fig. 7 as the dashed curve and red curve, respectively. They are the projections of the region inside the dashed curve and that in red. Clearly the selected halo and galaxy samples differ in the mean halo property $x$, lower for the galaxy sample. The mean values can be calculated as

[TABLE]

and

[TABLE]

where $t$ is determined by having the correct fraction $f$,

[TABLE]

The bias factors for the halo and galaxy samples are then

[TABLE]

and

[TABLE]

We can characterise the assembly bias effect by the fractional difference between the bias factor of the selected halo/galaxy sample and the average halo bias factor at $M_{\rm h}$ (or $V_{\rm peak}$), $\delta_{h}^{b}=(b_{h}-b_{c})/b_{c}$ and $\delta_{g}^{b}=(b_{g}-b_{c})/b_{c}$. Based on equations (4)–(8), we have

[TABLE]

That is, the assembly bias effect of the galaxy sample is weaker than that of the halo sample, by a factor equal to the correlation coefficient of the galaxy and halo property. Only in the case that galaxy and halo properties are tightly correlated (with zero scatter; $|\rho|=1$) does halo assembly bias effect completely translate to galaxy assembly bias effect. The connection in equation (9) is also valid for samples defined by the property range bounded by two percentiles (i.e. bin samples rather than threshold samples considered here).

To test how well the simple model works, we choose a pair of halo and galaxy properties to construct the halo and galaxy samples. We then measure the two-point correlation functions (2PCFs) of the halo and galaxy samples in each $M_{\rm h}$ and $V_{\rm peak}$ bin. To reduce the uncertainty, the large scale bias factor of a given halo sample is derived from the ratio of the halo-matter two-point cross-correlation function and the matter auto-correlation function (e.g. Xu & Zheng, 2018), $b_{h}=\xi_{hm}(r)/\xi_{mm}(r)$, averaged over scales of 5-18 $h^{-1}{\rm Mpc}$. The bias factor for the galaxy sample is similarly derived. We consider samples based on halo formation time $a_{\rm M/2}$ and galaxy colour/SFR/sSFR, which show $M_{\rm h}$ ($V_{\rm peak}$) dependent correlation coefficient (Fig. 6). The results are shown in Fig. 8. In the top-left panel, the thick dotted red (blue) curves show the assembly bias quantity $\delta_{h}^{b}$ for the 50 per cent oldest (youngest) haloes as a function of $M_{\rm h}$. The solid red (blue) curves are $\delta_{g}^{b}$ for the 50 per cent reddest (bluest) central galaxies. Both quantities decrease with increasing halo mass, i.e. the assembly bias effect becomes weaker for more massive haloes. The thin dashed curves are the same as the dotted curves but modulated by the correlation coefficient between colour and $a_{\rm M/2}$, i.e. $\rho\delta_{h}^{b}$, the prediction for $\delta_{g}^{b}$ from the simple model. The model works well in reproducing the halo mass dependent galaxy assembly bias effect based on halo assembly bias and the galaxy-halo correlation. The middle-left and bottom-left panels are for galaxies selected according to SFR and sSFR. The right panels show the assembly bias effect as a function of $V_{\rm peak}$. In all the cases, $\delta_{g}^{b}$ can be well described by $\rho\delta_{h}^{b}$, which supports the effectiveness of the simple model.

The success of the model suggests that an easy recipe could be developed to incorporate galaxy assembly bias into the halo model. The contribution of central galaxies to the galaxy bias factor, in its full form in the simple model, is

[TABLE]

where $t=t(f)$ is from equation (6). As an example, let us use the halo mass $M_{\rm h}$ as the primary variable in the halo model and consider a galaxy property (colour) that correlates with halo formation time ($a_{\rm M/2}$ or \log$$a_{\rm M/2}). In the halo model, besides the average halo bias $b_{c}$, we also need to know how the halo bias changes with $a_{\rm M/2}$ ($\partial b/\partial x$), the mean value of $a_{\rm M/2}$ ($x_{c}$), and the scatter in $a_{\rm M/2}$ ($\sigma_{x}$), all as a function of $M_{\rm h}$. As usual, we can construct fitting formulae for those four quantities based on $N$-body simulations. The quantities $f$ and $\rho$ belong to the description of the galaxy-halo relation, which can be parameterised. For $f$, it is simply the occupation fraction for haloes at $M_{\rm h}$. For $\rho$, the results in Fig. 6 suggest that a quadratic form with $M_{\rm h}$ would suffice. To compute the galaxy bias factor for a galaxy sample, we only need to perform a 1D integral over halo mass. Certainly it is straightforward to include assembly bias effect by populating dark matter haloes in $N$-body simulations. However, the above proposal has its virtue for analytic calculations in theoretical investigations.

The simple model can also be applied to observation to infer the correlation between galaxy and halo properties. Lin et al. (2016) construct samples of central galaxies from the Sloan Digital Sky Survey data and study the assembly bias effect from the 2PCF measurements. Early and late galaxy samples are defined according to either star formation history (SFH) or sSFR. Weak lensing measurements are used to verify that the host haloes of the early and late galaxies are of similar halo mass (around $10^{12}h^{-1}{\rm M}_{\odot}$). Lin et al. (2016) compare the difference in the early and late galaxy clustering to that in the early and late formed haloes, and do not find evidence for galaxy assembly bias. For SFH-based galaxy samples, they find the ratio of the early to late galaxy bias factor to be $1.00\pm 0.12$ (see their fig.5). If we take the mean of the bias factors of the two galaxy samples as the average halo bias and the uncertainty comes from two similar error bars added in quadrature, the measurement gives $\delta_{g}^{b}=0$ with an uncertainty 0.085. For halos around $10^{12}h^{-1}{\rm M}_{\odot}$, $\delta_{h}^{b}\sim 0.25$, with early- and late-forming haloes (Fig. 6). The coefficient of the correlation between SFH and halo formation time is then inferred to be $\rho=\delta_{g}^{b}/\delta_{h}^{b}=0.00\pm 0.34$. For the sSFR-based samples, the ratio of the early to late galaxy bias factor is $1.07\pm 0.14$ (their fig.5). We infer $\delta_{g}^{b}=0.034\pm 0.099$, and with $\delta_{h}^{b}\sim 0.25$ the coefficient of the correlation between sSFR and halo formation time is constrained to be $\rho=0.14\pm 0.40$. For both cases, the correlation between galaxy and halo property is consistent with being small. It implies that galaxy SFH and sSFR at most only loosely track halo formation. Similar measurements with large samples can reduce the uncertainty in the inferred correlation, which would help test galaxy formation models (e.g. by comparing to those in Fig. 6).

4 Summary and Discussion

Properties in galaxies residing in haloes of the same mass may have a dependence on certain aspects of halo assembly or formation history, which is termed as galaxy assembly bias. Studying galaxy assembly bias and its relation to halo properties can help improve the halo model of galaxy clustering and yield insights in galaxy formation and evolution. Using the Illustris cosmological hydrodynamic galaxy formation simulation, we investigate the central galaxy assembly bias effect through studying the relation among a set of galaxy and halo properties.

The main results can be summarised as follows.

• (1)

Central galaxy stellar mass $M_{*}$ has a tighter relation with $V_{\rm peak}$ than with $M_{\rm h}$, manifested by the smaller scatter in $M_{*}$ at fixed $V_{\rm peak}$ than that at fixed $M_{\rm h}$. Once the assembly effect of $M_{*}$ on $V_{\rm peak}$ is included, $M_{*}$ shows nearly no correlation with any other halo assembly properties.

• (2)

The correlations between halo assembly properties and other galaxy properties also appear cleaner if studied at fixed $V_{\rm peak}$, which reveal that galaxy SFR, sSFR, and colour mainly correlate with halo formation time (and to a less extent with halo concentration).

• (3)

A simple model is presented to show the relation between galaxy and halo assembly bias, which is linked by the correlation coefficient of the galaxy and halo property in consideration.

The Illustrius simulation produces a relation between central galaxy stellar mass and halo mass ($M_{*}$–$M_{\rm h}$) similar to that inferred from observation. We find that the scatter in the relation is closely related to halo assembly properties. For example, at fixed $M_{\rm h}$, haloes of higher $V_{\rm peak}$ or earlier formation tend to host galaxies of higher $M_{*}$. If we choose $V_{\rm peak}$ to be the primary halo variable, the scatter in $M_{*}$ at fixed $V_{\rm peak}$ is reduced compared to that at fixed $M_{\rm h}$. Remarkably, once switched to $V_{\rm peak}$, $M_{*}$ appears to have nearly no dependence on other halo assembly properties, at least for those considered in our study (including halo concentration, formation time, spin, accretion rate, and specific accretion rate). The property $V_{\rm peak}$, which is an indication of the maximum potential depth over the assembly history of haloes, is able to capture almost all the assembly effect in galaxy stellar mass. The results are in broad agreement with the study using the EAGLE simulation in terms of the $z=0$ maximum halo circular velocity (Matthee et al., 2017). The reason for $V_{\rm peak}$ to be the fundamental property in determining $M_{*}$ is likely related to the accretion and response of baryons in the gravitational potential. As for the scatter, the simulation noise in Illustris does not contribute much (Genel et al., 2019). It could be related to chaotic or stochastic processes in star formation and feedback (Matthee et al., 2017; Genel et al., 2019). Further study is needed to understand the cause of the correlation between $M_{*}$ and $V_{\rm peak}$ and the origin of the scatter.

We present the correlation between each pair of galaxy and halo property in terms of the Pearson correlation coefficient. Besides $M_{*}$, the other galaxy properties (SFR, sSFR, and colour) show a more consistent and clear trend with $V_{\rm peak}$ than with $M_{\rm h}$. At fixed $V_{\rm peak}$, those other galaxy properties related to star formation are found to mainly correlate with halo formation time and concentration, with stronger correlation with the former. The relatively nice behaviour in the correlations with galaxy properties in the $V_{\rm peak}$-based investigation suggests that it would be advantageous to use $V_{\rm peak}$* as the primary variable* in halo model of galaxy clustering, in particular in modelling stellar-mass-based samples (while the combination of or interpolation between halo mass and circular velocity quantities remains as a possibility to investigate, in the spirit of Lehmann, Mao, Becker, Skillman & Wechsler 2017). To further model SFR-, sSFR-, or colour-selected samples of galaxies, our investigation suggests to introduce halo formation time as the secondary halo variable (and to a less extent, halo concentration). This is in line with the age-matching model of Hearin & Watson (2013), who assumes monotonic mapping between galaxy colour and some variant of halo formation time at fixed galaxy luminosity (or stellar mass). That is, there exists a perfect correlation between the galaxy and halo property. Our investigation shows, however, that the correlation coefficient between galaxy SFR/sSFR/colour and halo formation time should be included as one important ingredient in the model.

The Illustris simulation is able to reproduce the observed star-forming main sequence reasonably well (Sparre et al., 2015). Here we find that at fixed $M_{\rm h}$ or $V_{\rm peak}$ the relation between SFR and $M_{*}$ of central galaxies follows the star-forming main sequence. Interestingly, compared with the SFR–$M_{*}$ relation, SFR shows a tighter correlation with sSFR ($\rho\sim 0.9$) over the full halo $M_{\rm h}$ or $V_{\rm peak}$ range. It is necessary to test its validity with observations and study its origin by tracking SFR and stellar mass growth of individual galaxies in simulations. The correlation between SFR and sSFR or between any pair of galaxy properties ($M_{*}$, SFR, sSFR, and colour) considered here cannot be explained solely by their common dependence on one halo assembly variable. Baryonic processes in galaxy formation (like star formation and feedback) likely play a major role in shaping such correlations.

For the effect of assembly bias on galaxy clustering, at fixed $M_{\rm h}$ or $V_{\rm peak}$, we come up with a simple model to relate the bias factors of a galaxy sample and the corresponding halo sample, which are connected by the correlation coefficient of the galaxy and halo properties used to define the two samples. It gives a reasonable description for the samples constructed with the simulation. It suggests a simple prescription to incorporate galaxy assembly bias into the halo model. By applying the simple model to the galaxy clustering measurements in Lin et al. (2016), we infer that the correlation between SFH/sSFR and halo formation time is consistent with being weak ($\rho\sim 0$–0.14). The simple model can be further tested with other hydrodynamic simulations, like EAGLE (Schaye et al., 2015) and IllustrisTNG (Nelson et al., 2019), and semi-analytic models, which can also provide further insights on parameterising the correlations between galaxy and halo properties. While our study in this paper focuses on central galaxies, we plan to carry out similar investigations for satellite galaxies to complete the picture of galaxy assembly bias at the halo level to help improve the halo model.

Acknowledgements

We thank David Weinberg, Jeremy Tinker and Hélion du Mas des Bourboux for useful comments. The support and resources from the Center for High Performance Computing at the University of Utah are gratefully acknowledged. This research was supported by the Munich Institute for Astro- and Particle Physics (MIAPP) of the DFG cluster of excellence “Origin and Structure of the Universe”, and ZZ thanks the hospitality of MIAPP, where the revision of the manuscript was undertaken.

Appendix A The $M_{*}$–$V_{\rm max}$Relation and the Baryon Effect on the $M_{*}$–$V_{\rm peak}$Relation

In section 3.1, we show that central galaxy $M_{*}$ has a tight correlation with $V_{\rm peak}$ and there is no further dependence on $a_{\rm M/2}$ at fixed $V_{\rm peak}$. The other velocity quantity commonly considered is $V_{\rm max}$, the present-day maximum circular velocity of a halo. While both $V_{\rm peak}$ and $V_{\rm max}$ characterise the depth of the potential well of a halo, the computation of $V_{\rm max}$ in a simulation is more straightforward, without tracking the whole formation history of a halo. We find that the $M_{*}$–$V_{\rm max}$ relation is similar to the $M_{*}$–$V_{\rm peak}$ relation, as shown in the left panel of Fig. 9. However, the $M_{*}$–$V_{\rm max}$ relation has a residual dependence on $a_{\rm M/2}$, with on average higher stellar mass in older haloes. This suggests that $V_{\rm peak}$ better captures the effect of halo formation time on the final central galaxy stellar mass than $V_{\rm max}$.

In the right panel of Fig. 9, we show the standard deviation in $\log M_{*}$ as a function of $V_{\rm peak}$ and $V_{\rm max}$, respectively. We fit a straight line to derive the mean relation between \log$$V_{\rm peak} and \log$$V_{\rm max}, and the correspondence between $V_{\rm peak}$ and $V_{\rm max}$ in the figure comes from such a fit. The $V_{\rm max}$ dependence of the scatter in $M_{*}$ only slightly differs from the $V_{\rm peak}$ dependence. In particular, compared to the $V_{\rm peak}$ dependence of the scatter in $M_{*}$, that for the $V_{\rm max}$ dependence is smaller (higher) for haloes with $V_{\rm peak}$ below (above) $\sim$300 ${\rm km\,s^{-1}}$. Although $V_{\rm peak}$ explains the dependence of $M_{*}$ on formation time better than $V_{\rm max}$, the tighter $M_{*}$–$V_{\rm max}$ relation in lower $V_{\rm peak}$ haloes indicates that $M_{*}$ could also depend on variables other than halo formation time (e.g. halo environment).

In section 3.1, the quantity $V_{\rm peak}$ used in the analyses is from the full physics simulation of Illustris. One may wonder whether baryons are partly responsible for the tight $M_{*}$–$V_{\rm peak}$ relation, as gas cooling with subsequent star formation can deepen the potential and thus change $V_{\rm peak}$. We explore the baryon effect on the $M_{*}$–$V_{\rm peak}$ relation by replacing $V_{\rm peak}$ and $a_{\rm M/2}$ with those in the dark-matter-only (DMO) simulation, after matching Rockstar haloes found in the Illustris-2-Dark simulation to those in the Illustris-2 simulation. In the left panel of Fig. 10, we show the ratio of full physics (hydro) $V_{\rm peak}$ to the DMO $V_{\rm peak}$ as a function of DMO $V_{\rm peak}$. At each DMO $V_{\rm peak}$ bin, we further divide haloes into two halves, with high and low $M_{*}$, respectively. Haloes with more massive central galaxies at fixed DMO $V_{\rm peak}$ tend to have higher hydro-to-DMO $V_{\rm peak}$ ratio, reflecting the effect of baryons on deepening the potential. The difference between hydro $V_{\rm peak}$ and DMO $V_{\rm peak}$ is larger in haloes of lower DMO $V_{\rm peak}$. However, the overall difference is small, e.g. below $\sim$3 per cent. This difference plays a negligible role in the $M_{*}$–$V_{\rm peak}$ relation, as shown by comparing the right panel of Fig. 10 to the bottom-right panel of Fig. 1.

Appendix B A More General Form of Central Galaxy Assembly Bias

In section 3.3 we develop a simple model to describe central galaxy assembly bias based on halo assembly bias and galaxy-halo relation. In that model, bias is assumed to depend only on one halo property at fixed halo mass. The model works reasonably well based on the tests with Illustris simulation. Since halo assembly bias is found to have a multi-variate dependence (e.g. Mao, Zentner & Wechsler, 2018; Xu & Zheng, 2018; Han et al., 2019), it is useful to extend the simple model for more general applications. Instead of increasing the dimensionality of the model and describing halo assembly bias as a function of multiple halo assembly variables at fixed halo mass, we choose to model halo bias as a function of one halo assembly property ($x$) and one galaxy property ($y$). That is, the bias dependence on any additional halo assembly properties is absorbed into that on the galaxy property (as galaxy property can be related to those halo properties). Additionally, the bias factor is assumed to be linearly dependent on halo and galaxy property, $b(x,y)=k_{x}x+k_{y}y+b_{c}$. Therefore, the contours of constant bias now have the possibility of deviating from vertical in the top panel of Fig. 7. The coefficient $k_{y}$ captures the dependence of $y$ on other halo properties and the scatter in the $y$-$x$ relation.

With this model, equations (7) and (8) become

[TABLE]

and

[TABLE]

where the new terms $\langle y\rangle_{x>t\sigma_{x}}$ and $\langle y\rangle_{y>t\sigma_{y}}$ can be calculated as

[TABLE]

and

[TABLE]

Similar to what we do in section 3.3 and based on equations (13), (14), (11), (12), (4), and (5), we obtain the following relation between the fractional halo and galaxy assembly bias ($\delta_{h}^{b}$ and $\delta_{g}^{b}$),

[TABLE]

where $\gamma=(k_{y}\sigma_{y})/(k_{x}\sigma_{x})$ is an indicator of the relative importance of the bias dependence on galaxy and halo property. For $\gamma$=0 (e.g. with $k_{y}=0$), this model reduces to the simple one in section 3.3, with the bias gradient along the $x$ direction (bias contours being vertical). In general, $|\gamma|<|\rho|$ corresponds to bias contours close to vertical, while $|\gamma|>|\rho|$ close to horizontal. In the following we will inverstigate which case represents that in the simulation by calculating $\gamma$.

We examine the performance of equation (15) using the Illustris simulation and test the importance of $\gamma$. Same as in section 3.3, halo assembly property ($x$) we consider is the formation time $a_{\rm M/2}$, and galaxy properties ($y$) include colour, SFR, and sSFR. At a fixed $M_{\rm h}$ or $V_{\rm peak}$ bin, haloes are divided into four sub-samples, depending which lower/upper halves of $x$ and lower/upper halves of $y$ they belong to. We measure the bias factors of the sub-samples and derive $k_{x}$ and $k_{y}$ from fitting them with $b=k_{x}\langle x\rangle+k_{y}\langle y\rangle+b_{c}$. Together with the measured standard deviations $\sigma_{x}$ and $\sigma_{y}$, we infer $\gamma$ in equation (15.)

The results are shown in Fig. 11. It is similar to Fig. 8, except that a set of dash-dotted curves are added to show the predicted $\delta_{g}^{b}$ by equation (15) (and that larger bins are used for haloes of high $M_{\rm h}$/$V_{\rm peak}$ to reduce the noise in the bias measurements for sub-samples). Clearly the more general model of equation (15) leads to better agreements with the measured galaxy assembly bias, and in almost of the cases the predictions just fall on top of the measurements. In general, however, we find that $\gamma$ only plays a minor role in determining the fractional assembly bias, compared to $\rho$. For example, we find $\gamma=-0.13\pm 0.22$ with $a_{\rm M/2}$ and galaxy colour across all $M_{\rm h}$ bins, lower in amplitude than $\rho$ (about 0.3–0.5, as shown in Fig. 6). Similar results are also obtained with other galaxy properties that we used in Fig. 8. That is, as long as the halo assembly property we identify is closely connected to galaxy property, any additional dependence of bias (on other halo properties, through that on galaxy property) is weak. The result indicates that considering only the bias dependence of such a halo property is sufficient to describe the most useful part of halo bias in modelling galaxy bias and the simple model in section 3.3 would work reasonably well. It would be interesting to see how general this holds by testing with more galaxy and halo properties and by using different galaxy formation simulations.

If one would like to incorporate any additional bias dependence into the model, the counterpart of equation (10) for the contribution of central galaxies to the galaxy bias factor reads

[TABLE]

where $\Delta x_{c}=\gamma(\sigma_{x}/\sigma_{y})y_{c}$ and $\rho_{\rm eff}=\rho+\gamma$. Instead of parameterising $\rho$ (as a function of $M_{\rm h}$), we now only need to describe two quantities, the offset $\Delta x_{c}$ in halo assembly property and the effective correlation coefficient $\rho_{\rm eff}$. As with equation (10), the expression follows the spirit of halo model – $b_{c}$, $\partial b/\partial x$, $x_{c}$, and $\sigma_{x}$ are halo properties, while $\Delta x_{c}$, $\rho_{\rm eff}$, $t$, and $f$ describe the galaxy-halo relation, which includes the correlation and scatter between the galaxy and halo property and the galaxy occupation function.

To achieve more precise modelling, it is also possible to express halo bias as linear functions of more than two halo properties, instead of absorbing that into the dependence of galaxy property, and to include quadratic dependence of each halo property. While such extensions are straightforward in our framework, it requires multi-variate descriptions of halo and galaxy properties with a substantially increased number of parameters, making the model less tractable and less attractive. In practice, the simple model presented in section 3.3 would be a useful first step.