---
abstract: 'We examine the evolution of assembly bias using a semi-analytical model of galaxy formation implemented in the Millennium-WMAP7 N-body simulation. We consider fixed number density galaxy samples ranked by stellar mass or star formation rate. We investigate how the clustering of haloes and their galaxy content depend on halo formation time and concentration, and how these relationships evolve with redshift. At $z=0$ the dependences of halo clustering on halo concentration and formation time are similar. At higher redshift, halo assembly bias weakens for haloes selected by age, and reverses and increases for haloes selected by concentration, consistent with previous studies. The variation of the halo occupation with concentration and formation time is also similar at $z=0$ and changes at higher redshifts. Here, the occupancy variation with halo age stays mostly constant with redshift but decreases for concentration. Finally, we look at the evolution of assembly bias reflected in the galaxy distribution by examining the galaxy correlation functions relative to those of shuffled galaxy samples which remove the occupancy variation. This correlation functions ratio monotonically decreases with larger redshift and for lower number density samples, going below unity in some cases, leading to reduced galaxy clustering. While the halo occupation functions themselves vary, the assembly bias trends are similar whether selecting galaxies by stellar mass or star formation rate. Our results provide further insight into the origin and evolution of assembly bias. Our extensive occupation function measurements and fits are publicly available and can be used to create realistic mock catalogues.'
author:
- |
    S. Contreras$^{1,2,3}$, I. Zehavi$^{4}$, N. Padilla$^{2,5}$, C. M. Baugh$^{6}$, E. Jiménez$^{2}$, I. Lacerna$^{7,8}$\
    $^{1}$Centro de Estudios de Física del Cosmos de Aragón (CEFCA), Plaza San Juan 1, Planta-2, Teruel, E-44001, Spain\
    $^{2}$Institúto Astrofísica, Pontifica Universidad Católica de Chile, Santiago, Chile\
    $^{3}$Donostia International Physics Center (DIPC), Manuel Lardizabal pasealekua 4, 20018 Donostia, Basque Country, Spain\
    $^{4}$Department of Physics, Case Western Reserve University, Cleveland, OH 44106, USA\
    $^{5}$Centro de Astro-Ingeniería, Pontificia Universidad Católica de Chile, Santiago, Chile\
    $^{6}$Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham, DH1 3LE, UK\
    $^{7}$Instituto de Astronomía, Universidad Católica del Norte, Av. Angamos 0610, Antofagasta, Chile\
    $^{8}$Instituto Milenio de Astrofísica, Av. Vicuña Mackenna 4860, Macul, Santiago, Chile
bibliography:
- 'Biblio.bib'
---

\[firstpage\]

cosmology: theory — galaxies: evolution — galaxies: formation — galaxies: haloes — galaxies: statistics — large-scale structure of universe

Introduction {#Sec:Intro}
============

Cosmic structure evolves hierarchically in the cold dark matter model. Density fluctuations grow by gravitational instability and form dark matter haloes, which evolve via accretion and mergers with other haloes [@Press:1974]. [@WhiteRees:1978] formulated the basis of modern galaxy formation theory starting from this concept, postulating that galaxies form inside dark matter haloes via the cooling of gas, star formation and mergers of galaxies. This framework is the basis of semi-analytic models of galaxy formation (SAMs; see, e.g., @Baugh:2006 [@Benson:2010] for reviews). These models use the merger histories of dark matter haloes as the starting point to model galaxy formation. The first SAMs used merger trees constructed using Monte-Carlo approaches based on the extended Press-Schechter theory (e.g., @Lacey:1993 [@Cole:1994; @Kauffmann:1993]), while modern SAMs use merger trees extracted from high-resolution dark matter simulations (e.g., @Kauffmann:1999 [@Bower:2006; @DeLucia:2007; @Lagos:2008; @Benson:2012; @Jiang:2014; @Croton:2016; @Lagos:2018; @Stevens:2018]). This opens up the prospect of studying environmental influences on the formation histories and properties of dark matter haloes and the impact on the galaxies they host.

The framework that led to SAMs also underpins the development of the halo occupation distribution (HOD) approach as an empirical description of galaxy clustering (e.g., @Peacock:2000 [@Berlind:2002; @Cooray:2002; @Zheng:2005]). The HOD formalism characterizes the relationship between galaxies and dark matter haloes in terms of the probability distribution that a halo of virial mass $M_h$ contains $N$ galaxies of a given type, together with the spatial and velocity distribution of galaxies inside haloes. An assumed cosmology and a specified shape of the HOD then allows us to predict any galaxy clustering statistic. The HOD approach is a powerful way to interpret observed galaxy clustering measurements, essentially transforming correlation function measurements to the relationship connecting galaxies with haloes (e.g., @Zehavi:2011 [@Coupon:2012] and references therein). It is also a useful method to characterize the predictions of galaxy formation models in a concise form that allows us to quantify the galaxy-halo relation (e.g., @Zheng:2005 [@C13; @C17]). Another important application of the HOD approach is to facilitate the generation of realistic galaxy mock catalogues by populating dark matter haloes from an N-body simulation with galaxies that reproduce a particular target clustering measurement. This method has become increasingly popular due to the growing demand for such catalogues for planning for and interpreting the results from large galaxy surveys and due to its good performance and low computational cost (e.g., @Manera:2015 [@Zheng:2016]). In the standard HOD framework mass is the only halo property that plays a role. This foundation of the HOD method has its origins in the Press-Schechter formalism and the uncorrelated nature of the random walks used to describe halo assembly in excursion set theory. This leads to the prediction that the halo environment is correlated with halo mass but not with how the halo is assembled [@Bond:1991; @Lemson:1999; @White:1999]. This is, however, not the case for haloes in $N$-body simulations where halo populations of the same mass but with a different ‘secondary property’ display different clustering, an effect that is now generally termed [*(halo) assembly bias*]{}. This was convincingly demonstrated in the Millennium N-body simulation of [@Springel:2005] by [@Gao:2005] who showed the age-dependence of the clustering of haloes of the same mass (see also @Sheth:2004); this dependence of halo clustering on secondary properties besides mass was later extended to, e.g., concentration, spin, substructure (e.g., @Wechsler:2006 [@Gao:2007; @Jing:2007; @Lacerna:2012; @Xu:2017; @Villarreal:2017; @Mao:2018]).

[@Croton:2007] used a SAM applied to the Millennium Simulation to show that halo assembly bias also impacts the clustering of galaxies, an effect that is now commonly referred to as [*galaxy assembly bias*]{}, namely halo assembly bias as reflected in the galaxy distribution (see also @Zhu:2006 [@Zu:2008; @Lacerna:2011; @Chaves-Montero:2016]). This can potentially have important implications for interpreting galaxy clustering using the HOD framework (e.g., @Zentner:2014). Detecting galaxy assembly bias has proven challenging and controversial. Despite some studies which claim to have uncovered the existence of assembly bias in the observable Universe (e.g., @Berlind:2006 [@Yang:2006; @Cooper:2010; @Wang:2013b; @Lacerna:2014; @Hearin:2015; @Miyatake:2016; @Saito:2016]) others argue that the impact of assembly is small [@Abbas:2006; @Blanton:2007; @Tinker:2008; @Tinker:2011; @Lin:2016; @Zu:2016a; @Dvornik:2017] or that the assembly bias signal could be a result of different systematics (e.g. @Campbell:2015b [@Zu:2016b; @Zu:2017; @Busch:2017; @Sin:2017; @Tinker:2017a; @Lacerna:2018]).

This is the latest in a series of papers examining the spatial distribution of galaxies predicted by SAMs. [@C13] examined the clustering and HOD predicted by SAMS from different groups and found that the models give robust clustering predictions when the galaxies are selected by properties that scale with the halo mass (such as stellar mass). [@C15] studied how predicted galaxy properties (such as stellar mass, cold gas mass, star formation rate, and black hole mass) correlate with their host halo mass in different SAMs. [@C17] examined how the predicted HOD form evolves with redshift in SAMs. We proposed a parametric form for the evolution of the HOD fitting parameters that can be used when constructing mock galaxy catalogs or for consistently fitting clustering measurements at different epochs.

Finally, in [@Zehavi:2017] (hereafter Z18) we use SAMs to investigate how the galaxy content of dark matter haloes is influenced by the large-scale environment and halo age at $z=0$, for galaxy samples selected by their stellar mass, finding distinct variations of the halo occupation functions. We show that haloes which form early have more massive central galaxies, and thus start hosting them at lower halo mass, and fewer satellite galaxies, compared to late-forming haloes. We also find similar results in hydrodynamical simulations [@Artale:2018]. These [*occupancy variations*]{}, namely the dependence of the HOD on halo properties other than mass, are intimately related to assembly bias, as it is their effect combined with [*halo*]{} assembly bias that gives rise to [*galaxy*]{} assembly bias.

Here, we build on our previous studies and investigate the evolution of assembly bias and specifically the occupancy variations in SAMs. We extend the analysis of Z18 in a number of ways: 1) we study a wide range of redshifts between $z=0$ and $z=3$; 2) we explicitly examine separately the different manifestations of assembly bias, namely halo assembly bias, occupancy variation, and galaxy assembly bias; 3) we consider galaxy samples constructed using two properties, stellar mass and star formation rate (SFR); and 4) we select haloes using two secondary parameters, halo formation time and concentration. We use the [@Guo:2013a] SAM which is a recent galaxy formation model from the Munich group implemented in a Millennium class N-body simulation with a WMAP-7 cosmology.

[@Wechsler:2006] and [@Gao:2007] study the evolution of halo assembly bias in large $N$-body simulations using a mark-correlation statistic and the large-scale bias of the mass-halo cross-correlation, respectively. [@Hearin:2016] examine the redshift dependence of assembly bias in the context of an extension of the HOD framework that incorporates assembly bias (the so-called decorated HOD), finding that the impact of assembly bias on galaxy clustering weakens at higher redshift for samples with fixed stellar mass. We aim to comprehensively investigate the evolution of galaxy assembly bias using a physical galaxy formation model. We focus here on galaxy assembly bias as reflected in the halo occupation and galaxy clustering. To our knowledge this is the first work that explicitly examines the evolution of the occupancy variation, and as a consequence, of galaxy assembly bias. Our aim is to investigate the origin and evolution of assembly bias. This will enable the development of more sophisticated tests to search for assembly bias in the observable Universe. Our results will also help shape the design of new mock galaxy catalogues, which are necessary for the next generation of galaxy surveys.

The outline of this paper is as follows: in Section 2 we introduce the SAM used and describe the different galaxy and halo samples employed in this work. Section 3 shows our results regarding the evolution of halo assembly bias, while Section 4 presents our main results regarding the evolution of the occupancy variation. In Section 5 we study the impact of assembly bias on galaxy clustering and the evolution of galaxy assembly bias. Finally, in Section 6 we summarise our results and present our conclusions. We describe our publicly available occupancy variation measurements and parametric fits in the appendix. Throughout the paper masses are measured in $h^{-1}\, {\rm M_{\odot}}$, the SFR is measured in $\rm M_{\odot}/yr$ and distances are measured in $h^{-1}\, {\rm Mpc}$ and are in comoving units.

Theoretical background and sample definition {#Sec:GC}
============================================

In this section we describe the dark matter simulation and the semi-analytic model used in this paper. We also present the different galaxy and halo samples we utilize. Finally, we describe the techniques used to characterise the galaxy and halo samples.

The semi-analytic model {#Sec:Sam}
-----------------------

Semi-analytical modelling (SAMs) is one of the main tools used to study galaxy formation (see @Baugh:2006 [@Lacey:2016] for reviews). These models aim to follow the main physical processes involved in the formation and evolution of galaxies. Some of the processes modelled by the SAM are (i) the collapse and merging of dark matter haloes (ii) shock heating and radiative cooling of gas (iii) star formation (iv) supernovae, AGN, and photoionization feedback (v) chemical enrichment of gas and stars (vi) disc instabilities and (vii) galaxy mergers.

The SAM used here is that of Guo et al. (2013; hereafter G13). This model is a version of [L-GALAXIES]{} the SAM code developed by the Munich group [@DeLucia:2004; @Croton:2006; @DeLucia:2007; @Bertone:2007; @Guo:2011; @Henriques:2013; @Henriques:2015]. For an extended description of this model and its performance we refer the reader to Guo et al. (2013, see also @Guo:2016b and @C17). The outputs are publicly available from the Millennium Archive[^1]. G13 is the latest publicly available SAM of the Munich group that makes use of the Millennium-WMAP7 dark matter simulation. We will explore other SAMs in future work, but do not expect our conclusions to change.

N-body Simulation {#Sec:Mill}
-----------------

The G13 model is implemented in the Millennium-WMAP7 N-body simulation [@Guo:2013a; @Gonzalez:2014; @Lacey:2016]. This simulation has similar specifications to the original Millennium simulation of [@Springel:2005] but uses a WMAP7 cosmology (ie. $\Omega_{\rm \Lambda0}$ = 0.728, $\Omega_{\rm m0}$= $ \Omega_{\rm dm0}$+$\Omega_{\rm b0}$ = 0.272, $\Omega_{\rm b0}$ = 0.0455, $\sigma_8$ = 0.81, $n_{\rm s}$ = 0.967, $h$ = 0.704.). The simulation uses $2160^3$ particles in a periodic box of comoving volume $(\,500\, h^{-1}\, {\rm~Mpc})^3$ corresponding to a particle mass of $ 9.31\times10^8 h^{-1} \, {\rm M_{ \odot}}$ and a softening value of $5\, h^{-1}\, {\rm~kpc}$. There are 61 simulation snapshots output between $z=50$ and $z=0$.

Halo merger trees are constructed from the simulation outputs. Haloes are identified using the friends-of-friends ([FoF]{}) group finding algorithm [@Davis:1985] at each snapshot of the simulation, using a minimum of 20 particles per halo (equivalent to a mass of $ 1.86\times10^{10} h^{-1} \, {\rm M_{ \odot}}$). [SUBFIND]{} is then run on these groups to identify subhaloes [@Springel:2001]. Merger trees are constructed by linking a subhalo in one snapshot to a single descendant subhalo in the subsequent output, i.e., a subhalo merger tree. The semi-analytical code uses these merger trees as the starting point to build its galaxy catalogue. Here, the mass of a dark matter halo, $\rm M_{h}$, is defined as the mass within the radius where the halo overdensity is 200 times the critical density of the simulation (referred to as “m\_crit200” in the public database).

The galaxy and halo samples {#Sec:Samples}
---------------------------

### Classifying samples by galaxy properties

For the main part of our analysis we use samples defined by galaxy number density. To do this we rank the model galaxies either by stellar mass or SFR and include all galaxies above a particular value of the stellar mass or SFR threshold that provides the desired number density. We construct galaxy samples for three different number densities, $n = 0.01, 0.00316$ and $ 0.001 
\,h^{3}\, {\rm Mpc}^{-3}$, and for a wide range of redshifts: $z=0, 0.5, 1, 1.5,
2, 2.5$ and $ 3$. The samples are chosen to be evenly spaced in logarithmic number density with differences of half a decade in log abundance.

The cumulative comoving number density of galaxies ranked by stellar mass is often used to link galaxy populations across cosmic time (e.g., @Padilla:2010 [@Leja:2013; @Mundy:2015; @Torrey:2015; @C17]). This type of selection is preferred over using a constant stellar mass cut to select galaxies at different epochs since it mitigates the need to assume a specific evolution model for the stellar mass, is insensitive to systematic shifts in the calculation of stellar masses and can be readily applied to observations. It also facilitates the comparison with galaxy samples selected using different properties (here, e.g., with galaxies selected by their SFR). [@C13] also showed that the HOD predictions for samples defined in this way are robust among different SAMs at a fixed redshift.

Fig. \[Fig:AMF\] shows the cumulative stellar mass function (top panel) and SFR function (bottom panel) for all redshifts studied here. The horizontal dashed lines show the different number density cuts we consider. The galaxies selected in each case are those to the right of the intersection with their associated dashed line. The top panel exhibits the expected growth of the galaxy stellar mass with time, while the bottom panel shows that there are fewer star forming galaxies at low redshifts than at high redshift.

![The cumulative stellar mass function (top panel) and the cumulative SFR function (bottom panel) predicted by the G13 SAM for different redshifts (as marked). The dashed horizontal lines indicate the number densities of the samples used in this work.[]{data-label="Fig:AMF"}](Fig1a.pdf "fig:"){width="50.00000%"} ![The cumulative stellar mass function (top panel) and the cumulative SFR function (bottom panel) predicted by the G13 SAM for different redshifts (as marked). The dashed horizontal lines indicate the number densities of the samples used in this work.[]{data-label="Fig:AMF"}](Fig1b.pdf "fig:"){width="48.00000%"}

### Classification by halo properties

To investigate assembly bias we define subsets of the fixed number density galaxy samples by selecting haloes using two different intrinsic or secondary properties: formation time (age) and concentration.

We define the formation time of a halo as the redshift when its main progenitor reaches half of the halo’s present-day mass for the first time. This definition is commonly used in the study of assembly bias (e.g., @Gao:2004 [@Gao:2005; @Croton:2007], Z18). We note that the formation time of a halo is calculated at each redshift independently. We calculate the formation time using the merger trees available in the database and linearly interpolate the mass of the haloes between snapshots.

The other halo property we consider is the concentration. The halo concentration characterizes the density profile. It is canonically defined as $c_{\rm vir} = r_{\rm vir}/r_{\rm s}$, where $r_{\rm vir}$ is the virial radius of the halo and $r_{\rm s}$ is the inner “transitional” radius appearing in the @NFW:1996 profile, at which the density profile changes slope. It is often alternatively defined via the rotation curve of the halo, as the ratio between $V_{\rm max}$ and $V_{\rm vir}$, where $V_{\rm max}$ is the peak value of the rotation curve, $V^2_{\rm c} = G\ M(r)/r$, and $V_{\rm vir}$ the virial velocity of the halo, $V_{\rm vir} \equiv V_c(r_{\rm vir})$ [@Bullock:2001; @Gao:2007]. We utilize the latter definition here, which is directly calculable from simulation data and does not require any model fitting.

In order to explore the variation in clustering and halo occupation with halo age and concentration, following Z18, we rank the haloes by these properties and identify the 20 per cent oldest and youngest haloes (based on their formation time) and (separately) the 20 per cent of haloes that are most or least concentrated. These divisions are made in 0.1 dex bins of halo mass, so as to factor out the influence of the changing halo mass function on these properties; the 20 per cent extremes of the distribution set up in this way effectively have the same mass function as the overall sample. We also tested using binnings of 0.05 and 0.2 dex in halo mass finding no difference in our main results.

The HOD and the correlation function {#Sec:HOD_CF}
------------------------------------

To study the impact of assembly bias on galaxies we measure the halo occupation functions and the correlation functions for the various halo and galaxy samples.

The HOD formalism describes the “bias” relation between galaxies and mass at the level of individual haloes allowing us to characterize the galaxy-halo connection. The key ingredient is the halo occupation function, $\langle N(M_h) \rangle$, which represents, for a given galaxy sample, the average number of galaxies per halo as a function of halo mass (loosely referred to here also as the HOD). The commonly assumed shape for the halo occupation function is motivated by predictions of physical models such as SAMs and hydrodynamic simulations [@Berlind:2003; @Zheng:2005]. When inferring the HOD it is often useful to consider separately the contribution from central galaxies and that of the additional satellite galaxies populating the halo [@Kravtsov:2004; @Zheng:2005]. For stellar mass (or luminosity) threshold galaxy samples, the expected form of the central galaxies occupation function is a smoothed step function and roughly a power-law for the satellites. For samples defined by SFR or color, the shape of the halo occupation function is more complex to account for the paucity of blue/star forming centrals in massive haloes (e.g., @Zehavi:2005 [@Geach:2012; @C13; @Gonzalez:2018]). We emphasize that the HODs presented in this work are all calculated directly from the SAMs, rather than inferred from the clustering, as is commonly done with observational data.

The correlation function (CF) is the most fundamental measure of the spatial distribution of haloes and galaxies. It is defined as the excess probability of finding a pair of objects at a given separation compared to a random distribution. Following [@Gao:2007] and Z18, whenever we measure the CF for the full galaxy sample we calculate the auto CF (the correlation of a given sample of objects with respect to the same sample). In contrast, when we measure the CF of a subsample of galaxies (e.g., the ones associated with the 20 per cent earliest-formed haloes) we measure the cross CF between this sample and the full galaxy sample. As explained in Z18 (see specifically their Appendix B) using the cross CF increases the signal-to-noise of the measurements and facilitates the interpretation of the results compared with the use of the auto CF of the subsamples.

The evolution of halo assembly bias {#Sec:HAB}
===================================

There are two basic ingredients necessary for galaxy assembly bias: (i) halo assembly bias, namely the dependence of halo clustering on halo properties other than mass, and (ii) the variation in the galaxy content of haloes with these properties, which we refer to as the occupancy variation (see Z18). Galaxy assembly bias requires both effects to be present. In this paper, we study how halo assembly bias, the occupancy variation and the resulting galaxy assembly bias evolve with time. This will provide further insight into the nature and origin of assembly bias and may guide attempts to detect it in observational galaxy samples. We show the evolution of halo assembly bias in this Section. The evolution of the occupancy variation is discussed in Section \[Sec:OV\], and the evolution of galaxy assembly bias is presented in Section \[Sec:GAB\].

First we look at the evolution of halo assembly bias in the dark matter-only $N$-body simulation without reference to the SAM galaxies. We begin with a visual inspection of the distribution of haloes in the simulation. Fig. \[Fig:CW\] shows haloes in a slice of the Millennium-WMAP7 simulation at $z=0$ and $z=3$, distinguishing between those with early and late formation times and also those with high and low concentrations. Starting with the halo age dependence at $z=0$ (Fig. \[Fig:CW\], top-left double panels), we see that while both early-formed and late-formed haloes trace the same cosmic web, the early-formed haloes present a sharper view of the web and appear somewhat more clustered. The view of the cosmic web when highlighting the extremes of halo concentration at $z=0$ (Fig. \[Fig:CW\], top-right panels) is reminiscent of that using halo formation time, though the clustering differences are slightly less apparent in this case. The bottom half of Fig. \[Fig:CW\] shows the distribution of haloes chosen similarly, but now at $z=3$. As expected, the haloes overall appear less clustered than at $z=0$. The differences between the early-formed and late-formed haloes (bottom-left panels) are much smaller in this case, and interestingly for the concentration (bottom-right panels), it appears that the haloes with low concentration are in fact now more clustered than the ones with high concentration. We quantify all of the trends discussed above shortly below using the CF.

These results are in agreement with those of [@Gao:2007], who found that the halo assembly bias signal does not depend on redshift when the halo samples are selected using a fixed cut in peak height ($\nu = \delta_c/\sigma(M_h)$, where $\sigma(M_h)$ is the root mean square linear overdensity within a sphere with mean mass $M_h$, and $\delta_c(z)$ is the linear overdensity threshold for collapse at redshift z). For an increase in the peak height from $\nu = 0.5$ to $\nu = 1.55$ (the peak height values of the minimum halo mass of our simulation at $z=0$ and $z=3$, respectively), [@Gao:2007] found that the difference in the clustering signal between early and late-forming haloes decreases (with early formed haloes being more clustered at low $\nu$). For haloes selected by their concentration, they showed that at low peak height, high concentration haloes are more correlated than late forming, low concentration haloes, and for high peak height, low concentration haloes are more correlated than high concentration haloes. These results are equivalent to the redshift evolution trend found for a fixed halo mass cut. The same can be concluded if a fixed cut in the nonlinear mass for collapse is used [@Wechsler:2006].

![image](Fig2a.pdf){width="47.00000%"} ![image](Fig2b.pdf){width="47.00000%"} ![image](Fig2c.pdf){width="47.00000%"} ![image](Fig2d.pdf){width="47.00000%"}

To better understand how the age and concentration of haloes correlate with one another at different redshifts, Fig. \[Fig:CP\] shows the joint distribution of halo concentration and formation time at $z=0$ (left) and $z=3$ (right). We show both the distribution of the full set of haloes (contours) and the median concentration as a function of halo age for three narrow bins of halo mass (lines and errorbars). The jags in the contours in the $z=3$ panel are artificial being caused by the limited time resolution of the Millennium-WMAP7 simulation outputs at high redshift. At $z=0$ there is a clear trend of concentration increasing with formation redshift (as shown by the solid lines). On the other hand, at $z=3$ there is little variation of halo concentration with formation redshift which suggests that the assembly bias effect with concentration and halo age might be different. [@Xu:2018] also look at the correlation of halo bias with different secondary halo properties (including formation time and concentration). They show that this correlation changes dramatically with halo mass. Therefore, it is unsurprising that the evolution of the clustering signal with formation time and concentration of the halo is different (see also @Mao:2018 who studied the dependence of clustering on several halo secondary properties).

The concentration of haloes with a small number of particles (fewer than 200 particles or $M_h < 10^{11.26} h^{-1} {\rm M_{ \odot}}$, ie. the red lines in Fig. \[Fig:CP\]) could be underestimated due to resolution effects [@Trenti:2010; @Paranjape:2017]. Most of the following analysis uses halo masses above this threshold and ranks samples in halo concentration at fixed mass, rather than using its actual value; we expect that our results should be unaffected by this possible source of systematics.

![image](Fig3a.pdf){width="47.00000%"} ![image](Fig3b.pdf){width="47.00000%"}

We now look into the CF of halo samples with different concentrations and formation times for a fixed number density after rank ordering the haloes in decreasing mass at $z=0$ and $z=3$ (Fig. \[Fig:CF\_Ev\]). We use a halo number density of $n = 0.00618 \, h^{3}\, {\rm Mpc}^{-3}$, which is comparable to the number density of central galaxies in the $n = 0.01 \, h^{3}\, {\rm Mpc}^{-3}$ galaxy sample (when selecting galaxies by their stellar mass). The equivalent halo mass (peak height) cut for these samples are $10^{11.75}\ h^{-1}M_{\odot}$ (0.76) for $z=0$ and $10^{11.24}\ h^{-1}M_{\odot}$ (2.01) for $z=3$. In each subplot of Fig. \[Fig:CF\_Ev\], the black line in the top panel denotes the auto CF of the full halo sample, while the red (blue) lines correspond to the cross CF of the full sample with the 20 per cent oldest (youngest) haloes in the top row of the figure, or the 20 per cent highest (lowest) concentration haloes in the bottom row (see § \[Sec:HOD\_CF\]).

![image](Fig4a.pdf){width="33.00000%"} ![image](Fig4b.pdf){width="33.00000%"} ![image](Fig4c.pdf){width="33.00000%"} ![image](Fig4d.pdf){width="33.00000%"} ![image](Fig4e.pdf){width="33.00000%"} ![image](Fig4f.pdf){width="33.00000%"}

We find that at $z=0$, haloes with early formation times and high concentrations are more clustered than ones with late formation times and low concentrations. This is the well-studied behavior of halo assembly bias (e.g., @Gao:2005 [@Wechsler:2006; @Gao:2007]). The halo assembly bias effect, as reflected by the clustering differences, is slightly stronger for the case of halo formation time than for halo concentration and extends to smaller separations. We note that as we are measuring here halo (instead of galaxy) clustering, the scales involved are all in the so-called 2-halo regime. At higher redshift (e.g., $z=3$) there is no difference in the clustering measured for haloes at the extremes of the formation time distribution and low concentration haloes are more clustered than high concentration haloes, reversing the trend seen at the present day.

We reach the same conclusions as already inferred from Fig. \[Fig:CW\], namely that the halo assembly bias signal (i.e. the difference between the red and blue lines) decreases with increasing redshift for halo samples selected by age. For concentration, the evolution of halo assembly bias is stronger in the sense that the clustering differences reverse at high redshift. These trends are in agreement with the evolution of the halo assembly bias signal found in the original Millennium simulation by [@Gao:2007] where they found the same trends when the peak height increases from 0.76 to 2.01, that is the increase of peak height from our samples between $z=0$ and $z=3$ (see @Wechsler:2006 for a comparison using the nonlinear mass for collapse.)

To understand the origin of this difference between using age and concentration as the secondary parameter we show in the right panel of Fig. \[Fig:CF\_Ev\] the CF of the main progenitors of the $z=0$ haloes selected at $z=1, 2$ and $3$. We call this sample the ‘tracking sample’. For this sample only, the secondary property halo labels (i.e. in terms of the extremes of concentration or formation time) refer to the $z=0$ descendants. We find different trends for different scales. At large scales ($>10\ h^{-1}{\rm Mpc}$), the tracking sample shows the same clustering trend as their descendants at $z=0$ but with a higher amplitude. The behaviour of the tracking sample cannot be easily related to the [@Wechsler:2006] or [@Gao:2007] results, since there is no fixed mass cut (or peak height or nonlinear mass for collapse) in the tracking sample. We interpret this non-evolution in the halo clustering as a negligible change in the comoving position and abundance of the haloes in this redshift range. This means that the evolution of halo assembly bias is not caused by a change in the clustering of haloes with extreme values of formation time or concentration. Instead we attribute the evolution of the assembly bias signal at a fixed halo number density to a shift in the ranking of haloes according to their mass and secondary property. This means that, for example, haloes with the highest concentrations at $z=0$ are not necessarily the ones with the highest concentrations at $z=3$.

We demonstrate this shift in the ranking of the haloes in Table \[Table:SampleEv\]. Here we show that fewer than $40$ per cent of the progenitors of $z=0$ haloes were part of the original sample at $z=1$. At $z=3$ this number decreases to $\sim 20$ per cent. This shift also explains the different evolution of halo samples selected by age and concentration. Table \[Table:C\_Age\] shows that while at $z=0$ there is a $\sim 40$ per cent overlap between members of the early (late) formation time halo sample and the high (low) concentration halo sample this number decreases to $\sim 10$ per cent at $z=3$.

                   z=0   z=1   z=2   z=3
  --------------- ----- ----- ----- -----
   Early forming   100   29    23    22
   Late forming    100   18    18    17
    High conc.     100   37    23    19
     Low conc.     100   22    14    11

  : The progenitors of a 100 representative haloes from the $n = 0.00618\ h^3\,  {\rm Mpc}^{-3}$ sample at $z=0$, drawn from the early and late formation time, high and low concentration samples, that were part of those samples at $z=0,1,2$ and $3$. []{data-label="Table:SampleEv"}

                                z=0   z=1   z=2   z=3
  ---------------------------- ----- ----- ----- -----
   Early forming - high conc.   40    29    19    10
    Late forming - low conc.    41    27    17    13

  : The haloes in common in a representative sample of 100 haloes for a fixed number density of $n = 0.00618\ h^3\, {\rm Mpc}^{-3}$ in the early formation time - high concentration samples and in the late formation time -low concentration samples. The numbers are calculated as the ratio between the number of haloes in the interception of the two samples ($\rm A \cap B $) and the number of haloes in the union of both samples ($A \cup B$).[]{data-label="Table:C_Age"}

Different trends are seen at intermediate and small separations in the right panels of Fig. \[Fig:CF\_Ev\]. The progenitors of early formation time and high concentration haloes are more correlated on small scales and less correlated on intermediate scales (compared to haloes with late formation times and low concentrations). These are not the focus of our work presented here, and we provide just some heuristic considerations. As early-formed haloes grow faster at higher redshifts it may be expected that they exhibit stronger clustering on small scales at $z=3$ (since they accrete mass from nearby structures). The stronger clustering on intermediate scales for haloes with late formation times may be explained in terms of these haloes accreting more mass at lower redshifts and the structures that will merge with these haloes being in their vicinity but not immediate proximity.

Since at $z=0$ there is a $\sim 40$ per cent overlap between halo samples selected by age and concentration (Table \[Table:C\_Age\]) we can assume this explanation is also valid for the main progenitors of the haloes selected by concentration.

One might be concerned that the agreement between the correlation functions of early and late formation time haloes in the top-middle panel of Fig. \[Fig:CF\_Ev\] could be an artifact of the limited time resolution of the Millennium-WMAP7 simulation at high redshifts. To check this we also calculated the correlation functions using the P-Millennium simulation [@Baugh:2019], a dark matter only simulation with over four times as many snaphots as the Millennium-WMAP7 run and with a better mass resolution. We find the same trends as those presented in this work, confirming that our results are not a product of the finite time resolution of the dark matter simulation used.

The occupancy variation evolution {#Sec:OV}
=================================

In this section we show the evolution of the occupancy variation in the SAM, i.e., how the dependence of the HOD on a secondary halo property varies with time. This may provide us more insight into the nature and origin of this phenomena.

Z18 showed that in SAMs, when selecting galaxies at $z=0$ by their stellar mass, the predicted HOD depends on halo formation time as well as halo mass. They found that haloes with early formation times tend to start being populated by central galaxies (the main galaxy of a dark matter halo) at lower masses than those with late formation times, but they have a lower number of satellites. [@Artale:2018] show that this is also the case in hydrodynamic simulations. We find that the above results also hold for other redshifts. This is shown, for example, in Fig. \[Fig:HOD\_single\], where we plot the HOD at $z=1$ for $n = 0.01 \, h^3\, {\rm Mpc}^{-3}$ for galaxies ranked by their stellar mass. The occupancy variation for both central and satellite galaxies is clearly evident.

![(Top) The halo occupation distributions at $z=1$ for a galaxy sample with $n=0.01\ h^3\, {\rm Mpc}^{-3}$. The lines show the contribution of all galaxies (solid), centrals galaxies only (dotted) and satellite galaxies only (dashed). The red and blue lines show the halo occupation distributions for the $20$ per cent earliest and latest forming haloes, respectively. The shaded region represents the jackknife errors using 10 subsamples, and is effectively negligible over most of the range. (Bottom) The ratio between the different halo subsets (with the oldest haloes in red and youngest ones in blue) and the full sample HOD for central galaxies (dotted line) and satellites (dashed line). The lines are plotted only when there are at least 20 haloes per halo mass bin of 0.1 dex width.[]{data-label="Fig:HOD_single"}](Fig5.pdf){width="49.00000%"}

Fig. \[Fig:HOD\_1\], shows the HOD for $z=0$ (top) and $z=3$ (bottom) for the same sample selection. In the left panels of Fig. \[Fig:HOD\_1\] the lines represent the contribution from the $20$ per cent of haloes with the earliest (red) and latest (blue) formation times, while the right panels show the contribution from the haloes with the $20$ per cent highest and lowest concentrations. We remind the reader that these halo subsamples are constructed by selecting the haloes in narrow bins of halo mass. At $z=0$, the predictions for the high (low) concentration samples are similar to those with early (late) formation times. This is consistent with what we found in Section \[Sec:HAB\], that the behaviour of these samples in terms of clustering is similar at $z=0$, but now extended to the halo occupation with galaxies. This similarity is no longer present at $z=3$. At this redshift, the occupancy variation for haloes selected by their age is qualitatively similar to that at $z=0$ (and at $z=1$; see Fig. \[Fig:HOD\_single\]). For haloes selected by concentration, the occupancy variation decreases somewhat for the satellite galaxies and it almost disappears for the central galaxies as we go to $z=3$. These trends also hold for other number density samples.

![(Left) Same as Fig. \[Fig:HOD\_single\] but for $z=0$ (top) and $z=3$ (bottom). (Right) Same as the left panels, but for haloes selected by concentration instead of formation time.[]{data-label="Fig:HOD_1"}](Fig6a.pdf "fig:"){width="24.00000%"} ![(Left) Same as Fig. \[Fig:HOD\_single\] but for $z=0$ (top) and $z=3$ (bottom). (Right) Same as the left panels, but for haloes selected by concentration instead of formation time.[]{data-label="Fig:HOD_1"}](Fig6b.pdf "fig:"){width="24.00000%"} ![(Left) Same as Fig. \[Fig:HOD\_single\] but for $z=0$ (top) and $z=3$ (bottom). (Right) Same as the left panels, but for haloes selected by concentration instead of formation time.[]{data-label="Fig:HOD_1"}](Fig6c.pdf "fig:"){width="24.00000%"} ![(Left) Same as Fig. \[Fig:HOD\_single\] but for $z=0$ (top) and $z=3$ (bottom). (Right) Same as the left panels, but for haloes selected by concentration instead of formation time.[]{data-label="Fig:HOD_1"}](Fig6d.pdf "fig:"){width="24.00000%"}

We repeat the above analysis for galaxies selected by SFR in Fig. \[Fig:HOD\_2\]. The overall shape of the HOD at $z=0$ for SFR-selected galaxies is different than for galaxies selected by stellar mass due to the tendency of high mass halos to host non star forming (red) central galaxies, as discussed in Section \[Sec:HOD\_CF\]. Interestingly, the “dip” feature diminishes as one goes toward higher redshifts, possibly due to having less time for quenching mechanisms to occur. By $z=3$ the overall shape of the HOD, and in particular the central galaxies contribution, is very similar for the SFR-selected galaxy samples and the stellar mass selected samples (see also @Orsi:2008). We have verified that the transition in the shape of the HOD between $z=0$ and $z=3$ is smooth with increasing redshift. A large set of HOD measurements for different redshifts and number densities is being released with this paper (see Appendix A for more details).

![Same as Fig. \[Fig:HOD\_1\] but for galaxies selected by their SFR instead of stellar mass.[]{data-label="Fig:HOD_2"}](Fig7a.pdf "fig:"){width="24.00000%"} ![Same as Fig. \[Fig:HOD\_1\] but for galaxies selected by their SFR instead of stellar mass.[]{data-label="Fig:HOD_2"}](Fig7b.pdf "fig:"){width="24.00000%"} ![Same as Fig. \[Fig:HOD\_1\] but for galaxies selected by their SFR instead of stellar mass.[]{data-label="Fig:HOD_2"}](Fig7c.pdf "fig:"){width="24.00000%"} ![Same as Fig. \[Fig:HOD\_1\] but for galaxies selected by their SFR instead of stellar mass.[]{data-label="Fig:HOD_2"}](Fig7d.pdf "fig:"){width="24.00000%"}

At $z=0$, for the SFR-selected galaxy samples, early forming and high concentration haloes have a lower number of satellite galaxies compared to haloes with late formation times or low concentrations (same as for galaxies selected by stellar mass). For the central galaxies, at low halo masses, early forming and high concentration haloes have a larger number of central galaxies, while for higher halo masses they have a lower number of central galaxies compared to haloes with late formation times or low concentrations. The latter trend perhaps arises since the central galaxies in the early-formed high-mass halos have more time to be impacted by star formation quenching. At $z=3$, the HODs for galaxies selected by SFR display the same trends as those for galaxies selected by stellar mass. The occupancy variation (i.e. the difference between the red and blue lines) stays roughly constant for haloes selected by age. The occupancy variation with halo concentration decreases with redshift for the HOD of the satellites and nearly diminishes for that of the central galaxies.

The full redshift evolution of the occupancy variation is captured in Fig. \[Fig:HOD\_ratio\_1\] and Fig. \[Fig:HOD\_ratio\_2\], where we show the ratios of the HODs (as in the bottom subpanels of Fig. \[Fig:HOD\_single\]-\[Fig:HOD\_2\]) for halos selected by age and concentration at all redshifts explored, for galaxies selected by stellar mass and SFR, respectively. Here we corroborate that, for galaxies selected by stellar mass, the magnitude of the central galaxies occupancy variation is constant with redshift for haloes selected by age and it significantly decreases with increasing redshift (nearly diminishing by $z=3$) for haloes selected by concentration. The occupancy variation for the satellites part of the HOD progressively decreases for either age or concentration. The overall shift of the ratios toward lower halo mass with increasing redshift reflects the expected redshift evolution of the HOD (as studied for example by @C17; see their Fig. 5). For SFR selected galaxies (Fig. \[Fig:HOD\_ratio\_2\]), the occupancy variations decrease for both age and concentration, with a more pronounced trend for the latter.

![image](Fig8a.pdf){width="45.00000%"} ![image](Fig8b.pdf){width="45.00000%"}

![image](Fig9a.pdf){width="45.00000%"} ![image](Fig9b.pdf){width="45.00000%"}

The different evolution of the occupancy variation with age and concentration indicates a different origin for these two effects. Even though they appear similar at $z=0$, they evolve differently. We will further investigate their nature and origin in future work (Zehavi et al., in prep.). It is interesting to note that the evolution of the occupancy variation shows different trends compared to the evolution of the halo assembly bias found in Section \[Sec:HAB\], where the signal decreased for the halo samples selected by age but not for those selected by concentration as in the occupancy variation case. Both effects will influence the evolution of the galaxy assembly bias signal, as we will now show.

The evolution of galaxy assembly bias {#Sec:GAB}
=====================================

In this section we show the effect of assembly bias on the galaxy correlation function at different redshifts. As we did in Section \[Sec:HAB\], we measure the auto CF for the full galaxy sample as well as the cross CF of the full sample with the given subsample (e.g. early/late formed haloes).

![The auto-correlation and cross-correlation functions at $z=1$ for a stellar mass selected galaxy sample with number density of $0.01\ h^{3}\, {\rm Mpc}^{-3}$. Dashed lines correspond to the shuffled galaxy sample (i.e., without occupancy variation) and solid lines represent the original sample. The red (blue) lines show the cross-correlation function between the galaxies that populate the $20$ per cent earliest (latest) forming haloes and the full galaxy sample. The black lines show the auto-correlation function for the full sample. The top panel shows the different correlation functions, the middle panel shows the cross-correlation functions (red and blue lines) divided by the auto-correlation function (black), and the bottom panel shows the ratio between different correlation functions from the model and the ones from the shuffled sample. The shaded region and errorbars represent the jackknife errors calculated using 10 subsamples.[]{data-label="Fig:single_CF"}](Fig10.pdf){width="48.00000%"}

To study the impact of assembly bias on the CF we shuffle galaxies among haloes of the same mass, following the approach of [@Croton:2007] and Z18. This consists of taking all haloes in a given bin of halo mass (0.1 dex wide in our case; we also tested using a bin width of 0.05 and 0.2 dex and found no major difference in our results) and randomly reassigning the galaxy population between these haloes. Central galaxies are located at the position of the central galaxies they replace (except if there is no galaxy in a halo; in which case the new galaxy is located at the potential minimum of the halo). The satellite galaxies are moved together with their original central galaxy and retain the same relative positions to it in their new halo. The shuffling removes any potential connection to the assembly history of the halos, and effectively transforms the HOD of any halo subsample (e.g., for a range of halo formation times or concentrations) to be the same as the total HOD (making, e.g., the red and blue lines of Fig. \[Fig:HOD\_single\] be the same as the black line). This new galaxy sample will have, by construction, no occupancy variation.

The CFs for a stellar mass selected sample with number density $n = 0.01\ h^{3}\, {\rm Mpc}^{-3}$ at $z=1$ is shown in top panel of Fig. \[Fig:single\_CF\]. The auto CF of the full sample is shown in black, and the red and blue lines are the cross CF for the $20$ per cent earliest and latest forming haloes, respectively. The dashed lines show the CF for the shuffled samples. The shaded region and errorbars represent the jackknife errors calculated using 10 realisations for the real and the shuffled samples. The middle panel shows the ratios between the cross CF of the subsamples and the auto CF of the full sample, for both original (solid) and shuffled (dashed) galaxy samples. The bottom panel shows the ratios between the different CFs measured for the original (unshuffled) galaxy samples and the corresponding shuffled ones.

A value above unity for the black line in the bottom panel of Fig. \[Fig:single\_CF\] indicates that the original sample has a larger CF than that measured for the shuffled sample. These differences are the manifestation of galaxy assembly bias [@Croton:2007]. As explained in Z18, this arises from the combined effect of the occupancy variation and halo assembly bias. The central galaxies occupancy variation indicates a preferential occupancy of early-formed halos. These halos are more clustered, thus leading to a stronger clustering signal on large scales. The significant clustering differences on small scales come about from the satellites occupancy variation, where the increased number of satellites in late-forming halos gives rise to a stronger clustering in the 1-halo regime.

It is interesting to note that the clustering of galaxies in the late-forming haloes is stronger on large scales than that for galaxies in early-formed haloes, as can be seen in the middle panel of Fig. \[Fig:single\_CF\]. This is the opposite to the results found by Z18 at $z=0$, implying that trend evolves with redshift. This again arises from the inter-related effects of halo assembly bias and the occupancy variation. The dashed lines in the middle panel correspond to the shuffled galaxy samples, and reflect the same halo assembly bias trends seen in Fig. \[Fig:CF\_Ev\], modulated by satellite galaxies. The central occupancy variation acts to slightly increase this ratio for the galaxies in late-forming haloes and decrease it for the galaxies in early-formed halos (see corresponding discussion in Z18, specifically their §5.3), thus likely resulting in the reversed clustering trend seen.

The evolution of this ratio is individually presented Fig. \[Fig:CF\_Ratio1\]. The left panels show the ratio between the CF measured for galaxies in haloes with the 20 per cent earliest/latest formation times and the CF of the full sample, for $z=0,1,2$ and $3$. The panels on right show the same but for haloes selected by concentration. Solid lines show the original SAM galaxies and dashed lines show the shuffled sample results (i.e., with no occupancy variation). For the shuffled samples, we can see that the difference between the CFs of galaxies in the earliest/latest-forming haloes decreases with increasing redshift while the difference for haloes with the highest/lowest concentrations is reversed, with increasingly stronger clustering found for the haloes with a lower concentration. These trends are consistent with the evolution of the halo assembly bias signal shown in Section \[Sec:HAB\], where we found a decrease (flip) of the of the difference in clustering on large scales for haloes selected by age (concentration). At high redshift, the differences in the clustering of the SAM galaxies in early and late formed haloes come mostly from the occupancy variation. This is opposite to the situation of galaxies that live in high and low concentrated haloes, where the differences in their clustering come mostly from halo assembly bias.

![image](Fig11a.pdf){width="49.00000%"} ![image](Fig11b.pdf){width="49.00000%"}

For the galaxy population predicted by the SAM, galaxies in haloes with late formation times or low concentrations show stronger large-scale clustering at higher redshifts. For the concentration case, the clustering signal becomes identical to that measured for the shuffled galaxies. This is expected since, as shown in Section \[Sec:OV\], the occupancy variation of haloes with concentration decreases strongly with increasing redshift. For galaxies in haloes selected by formation time, the clustering of the galaxies in the latest forming haloes is stronger than the shuffled galaxies, while for the galaxies in the earliest-forming haloes it is lower than for the shuffled case. This is again consistent with what we found in § \[Sec:OV\], with the occupancy variation with age persisting to higher redshifts, with (at each redshift) late (early) formation effectively shifting the occupation toward higher (lower) halo masses, thus changing the clustering.

On small scales, the dashed lines in Fig. \[Fig:CF\_Ratio1\] are identical for both age and concentration, since there is no halo assembly bias in that regime (1-halo scales). The SAM galaxies in late-forming haloes are more correlated than those in early-forming haloes on small scales, at all redshifts. This is due to the increased number of satellites in early versus late forming halos, which persists at all redshifts (as seen in Fig. \[Fig:HOD\_ratio\_1\]). Galaxies selected by halo concentration exhibit a similar behavior – galaxies in low concentration haloes are more clustered than those in high concentration ones – at small-to-intermediate scales. On very small scales (below $\sim 0.1
h^{-1}\, {\rm Mpc}$), though, this trend flips. One might have expected the same small-scale behaviour with concentration due to the similar satellite occupancy variation. However, the concentration differences impact the clustering as well. E.g., for the low-concentration sample, even with more satellite galaxies, they are likely less concentrated (since they trace the dark matter distribution) and as a consequence, less clustered on very small scales.

In Fig. \[Fig:CF\_Ratio2\] we show, for completeness, the corresponding evolution of the CFs for galaxies selected by their SFR. We obtain the same trends found for galaxies selected by stellar mass. We also analysed other number density samples and found similar results for the evolution of the galaxy CFs.

![image](Fig12a.pdf){width="49.00000%"} ![image](Fig12b.pdf){width="49.00000%"}

Finally, we consider the evolution of the galaxy assembly bias signal. As previously mentioned, galaxy assembly bias is quantified in terms of the ratio between the CF of a galaxy sample and that of a shuffled sample, where the relation to halo assembly has been erased, as shown by the black line in the bottom panel of Fig. \[Fig:single\_CF\]. Fig. \[Fig:CF\_Ratio4\] presents our measurements for three different number densities ($n= 0.01$, $0.00316$ and $0.001 \, h^{-3}\, {\rm Mpc}^3$) over a range of redshifts, for galaxies selected by stellar mass and SFR. We find that this clustering ratio generally decreases for higher redshifts, and for lower number densities. Interestingly, this decrease can be large enough in some of these cases so that the original sample becomes less clustered than the corresponding shuffled sample, and the clustering [*difference*]{} changes sign and continues growing in magnitude. This typically occurs for lower number densities and at high redshifts. This ratio is overall lower for the galaxy samples selected by their SFR rather than stellar mass. Nonetheless, the trends with redshift and number density persist for these SFR selected samples,

![image](Fig13a.pdf){width="49.00000%"} ![image](Fig13b.pdf){width="49.00000%"}

Our results are in qualitative agreement with those found by @Hearin:2016 ([-@Hearin:2016]; their Fig. 8) over the limited redshift range they explore ($0<z<1$). However, we note that @Hearin:2016 compare samples with the same stellar mass thresholds at the different redshifts, not accounting for any stellar mass evolution. Effectively, this amounts to probing more massive galaxies (lower number densities) at higher redshifts, and thus it is impossible to separate the evolution they find from the expected number density dependence.

Again, the impact of assembly bias on galaxy clustering arises from the combined contributions of the occupancy variation and halo assembly bias. At relatively low redshift for the stellar mass selected samples, these typically combine to produce an increased clustering (see the discussion following Fig. \[Fig:CF\_Ratio1\] and in Z18). For example, using our results in the previous sections for $n=0.01 h^{3}\, {\rm Mpc}^{-3}$, we can explain the behaviour exhibited in the top-left panel of Fig. \[Fig:CF\_Ratio4\]. As we saw earlier for the halo age case, the level of halo assembly bias decreases while the level of the central galaxies occupancy variation remains similar, leading to a diminishing galaxy assembly bias effect. For lower number densities at high redshifts, halo assembly bias reverses sense (in a similar manner to concentration) such that early-formed halos become less clustered than the late-formed ones, and this gives rise to the reversed sense of galaxy assembly bias in those cases.

One could have a-priory envisioned a scenario in which the stochasticity involved in the galaxy formation processes would serve to weaken galaxy assembly bias over time. Alternatively, one might have expected the signature to grow with time (i.e., diminish as one goes to higher redshift), due to the hierarchical growth of structure. However, it seems that the evolution of assembly bias is far more intricate. The overall trend we find is that the CFs ratio monotonically increases with time (or decreases with increasing redshift; Fig. \[Fig:CF\_Ratio4\]). This leads to a change in the sign of the effect, i.e., whether the clustering of the galaxy sample is stronger or weaker due to assembly bias effects, as well as a shift in whether the magnitude of this clustering difference decreases or increases with time. This gets more complex as the amplitude of the clustering ratio varies with the specifics of the galaxy selection (e.g., stellar mass or SFR) and number density, thus it is non-trivial to predict which galaxy sample would show negligible or extreme assembly bias properties and at which redshift.

Summary & Conclusions
=====================

We use a state-of-the-art semi-analytic model of galaxy formation, the G13 SAM model, to study the origin and evolution of assembly bias in the galaxy distribution. We identify two separate contributions to this effect: halo assembly bias, which refers to the different clustering of haloes with different ‘secondary property’, and occupancy variation, the dependence of the number of galaxies in haloes of the same mass on a second property of the haloes. We isolate the evolution of these two effects for haloes selected by their concentration and formation redshift, two of the most common secondary properties used to measure assembly bias. The galaxy samples correspond to different number densities based on either ranked stellar mass or SFR. Our key results are shown in Figures \[Fig:CF\_Ev\], \[Fig:HOD\_ratio\_1\] and \[Fig:CF\_Ratio4\]. We now summarise our main findings:

-   At $z=0$ the concentration of dark matter haloes correlates with formation time. This correlation weakens at higher redshifts.

-   Haloes at $z=0$ with high concentrations or early formation times are more clustered than those with low concentrations or late formation times. At high redshift, there are no differences in the CF measured for haloes with different formation times, but low concentration haloes are more correlated than high concentration ones.

-   Haloes ranked to have an extreme concentration or formation time at a given redshift do not necessarily have the same ranking at other redshifts. We found that the main progenitors of $z=0$ haloes display clustering similar to that measured for their descendants. This means that the evolution of the halo assembly bias signal is not caused because a set of haloes (e.g., high concentration haloes) change their clustering over time, but because haloes change their ranking in terms of a secondary property.

-   At $z=0$, haloes with early formation times or high concentrations are populated by galaxies starting at lower halo masses (for a fixed cut in stellar mass) but they have fewer satellite galaxies for a fixed mass compared to haloes with late formation times or low concentrations.

-   For galaxies selected by SFR we generally find similar occupancy variation trends to those found for galaxies selected by stellar mass (though different shape of the HOD). Haloes with early formation times or high concentrations are first populated by galaxies at a lower mass and have fewer satellite galaxies at a given mass compared to haloes with late formation times or low concentrations. The one difference is that at higher halo masses, where the central galaxies occupation drops, there are less centrals in haloes with early formation times or high concentrations than for those with either late formation times or low concentrations.

-   The occupancy variation for central galaxies in haloes with different formation times stay roughly constant as a function of redshift for a fixed galaxy number density and for galaxies selected by either stellar mass or SFR. The corresponding satellite galaxies occupancy variation decreases somewhat with increasing redshift.

-   The occupancy variation for galaxies in haloes with different concentrations diminishes for the central galaxies and satellites with increasing redshift, for both stellar mass or SFR selected galaxy samples.

-   The evolution of the CF of galaxy samples without occupancy variation (i.e., the shuffled samples) reflects the same trends on large scales as the evolution of halo assembly bias for haloes selected by age or concentration; the CF differences for galaxies in haloes with early and late formation times decreases with look back time, while the CF of galaxies in low-concentration haloes increases relative to the CF of galaxies in high-concentration haloes when going to higher redshifts.

-   The CF of galaxies hosted by haloes with late formation times or low concentration increases relative to the CF of galaxies in haloes with early formation times or high concentrations, respectively, with increasing redshift.

-   The occupancy variation tends to increase the amplitude of the CF of galaxies that live in haloes with either late formation times or low concentrations, and decrease it for galaxies that live in haloes with early formation times or high concentrations.

-   Galaxy assembly bias as measured by the ratio between the CF of the model galaxies and that of the shuffled galaxies decreases with redshift, going below 1 in some cases. This CFs ratio is generally smaller for lower number densities and for SFR-selected samples.

The different evolution of halo assembly bias and the occupancy variation with age and concentration likely points to a different origin for the dependence on these two secondary parameters. This is further corroborated by their lack of correlation at high redshift. In general, we find similar trends in the evolution of assembly bias, for both the occupancy variation and galaxy assembly bias, for galaxies selected by SFR versus stellar mass. This is quite impressive considering that galaxy samples selected by stellar mass and by SFR exhibit quite different behaviours in the SAMs [@C13; @C15], and may be relevant for upcoming surveys.

The results shown here will help to inform theoretical models of assembly bias and the development of observational tests to detect its existence (or absence) in the Universe. They can also be used to construct improved mock galaxy catalogues incorporating assembly bias (as standard HOD mocks do not include this effect). For these purposes we are releasing all the HODs and occupancy variation measures obtained in this work, as well as parametrised fits for them (see Appendix A for more details).

Acknowledgements {#acknowledgements .unnumbered}
================

We thank the anonymous referee for insightful comments that helped improve the presentation of this manuscript. We thank Peder Norberg for useful discussions. This work was made possible by the efforts of Gerard Lemson and colleagues at the German Astronomical Virtual Observatory in setting up the Millennium Simulation database in Garching. SC, IZ, NP, and EJ acknowledge the hospitality of the ICC at Durham University. SC & NP acknowledge support from a STFC/Newton-CONICYT Fund award (ST/M007995/1 - DPI20140114) and Anillo ACT-1417. SC is also supported by the European Research Council through grant ERC-StG/716151. IZ acknowledges support by NSF grant AST-1612085 and by a CWRU Faculty Seed Grant. NP & EJ are further supported by “Centro de Astronomía y Tecnologías Afines” BASAL PFB-06 and by Fondecyt Regular 1150300. This project has received funding from the European Union’s Horizon 2020 Research and Innovation Programme under the Marie Skłodowska-Curie grant agreement No 734374. The calculations for this paper were performed on the ICC Cosmology Machine, which is part of the DiRAC-2 Facility jointly funded by STFC, the Large Facilities Capital Fund of BIS, and Durham University and on the Geryon computer at the Center for Astro-Engineering UC, part of the BASAL PFB-06, which received additional funding from QUIMAL 130008 and Fondequip AIC-57 for upgrades.

HOD catalogues
==============

To facilitate the creation of mock galaxy catalogues with occupancy variation, which can be used for the creation of mocks with galaxy assembly bias, we are making public the HODs calculated in this work. The HODs are calculated for the following number densities, $n = 0.0316, 0.01, 0.00316, 
0.001, 0.000316$ and $0.0001 \, h^3\, {\rm Mpc}^{-3}$, for galaxies ranked either by stellar mass or SFR. The following redshifts are used, $z=0,0.5,1,1.5,2,2.5$ and $ 3$, and the haloes are selected in 10 bins of ranked age and concentration. This yields more than 1,800 HODs in total. This material can be found at <https://github.com/hantke/-HOD_Extractor2>

Additionally, we provide in this same repository the HOD fitting parameters for all galaxy samples selected by stellar mass, for the commonly used 5-parameter model introduced by @Zheng:2005. The HOD parameters of the galaxies selected by SFR can not be well represented by this standard parametrization (see, e.g., @C13) and will be investigated further in future work. The 5-parameters model is given by $$\langle N_{\rm gal}(M_{\rm h})\rangle = \langle N_{\rm cen}(M_{\rm h})\rangle + \langle N_{\rm sat}(M_{\rm h})\rangle,
\label{Eq:Cen_HOD}$$ with $$\langle N_{\rm cen}(M_{\rm h})\rangle = \frac{1}{2}\left[ 1 + {\rm erf} \left( \frac{\log M_{\rm h} - \log M_{\rm min}}{\sigma_{\log M}}  \right) \right]
\label{Eq:Cen_HOD2}$$ and $$\langle N_{\rm sat}(M_{\rm h})\rangle = \left( \frac{M_{\rm h}-M_{\rm cut}}{M^*_1}\right)^\alpha,
\label{Eq:Sat_HOD}$$ where $M_{\rm h}$ is the host halo mass and $ {\rm erf}(x)$ is the error function, $${\rm erf}(x) = \frac{2}{\sqrt{\pi}} \int_{0}^{x} e^{-t^2} {\rm d}t.$$ $M_{\rm min}$ is the mass where, on average, half of the haloes are occupied by a central galaxy (i.e., $\langle N_{\rm cen}(M_{\rm min})\rangle = 0.5$); $\sigma_{\log M}$ characterises the width of the transition from zero to one central galaxy per halo, where $\sigma_{\log M}=0$ represents a vertical step-function transition; $\alpha$ is the slope of the satellite HOD power-law; $M_{\rm cut}$ is the minimum halo mass at which haloes can host a satellite galaxy and $M^*_1$ is the normalization. Note that we provide instead the value of a related parameter, $M_{1}$, the halo mass where on average there is one satellite per halo (i.e., $\langle N_{\rm sat}(M_{1})\rangle = 1$) and is equal to $ M_{\rm cut} + M^*_1$.

As an example, we provide in Table \[Table:HOD\] the HOD parameters for the galaxy samples with $n=0.01\ h^3\, {\rm Mpc}^{-3}$ at $z=0,1,2,3$ and for the $10$ per cent oldest haloes, the $10$ per cent youngest haloes, and the full halo sample. The full set of fitted parameters is provided in our public release website. Together with the HODs and their parameters (stored in two [HDF5]{} files), we are also releasing tools to read and plot the HODs to facilitate their analysis.

                                                    $z=0$                                   
  -------------------------- --------------- ------------------- ---------- --------------- ---------
                              $M_{\rm min}$   $\sigma_{\log M}$   $\alpha$   $M_{\rm cut}$   $M_{1}$
    10$ \% $ oldest haloes        11.44             0.10            1.15         11.99        12.82
          All haloes              11.62             0.21            0.99         11.83        12.57
   10$ \% $ youngest haloes       11.92             0.30            0.84         11.74        12.31
                                                    $z=1$                                   
                              $M_{\rm min}$   $\sigma_{\log M}$   $\alpha$   $M_{\rm cut}$   $M_{1}$
    10$ \% $ oldest haloes        11.36             0.16            1.07         11.97        12.57
          All haloes              11.59             0.28            0.92         11.81        12.34
   10$ \% $ youngest haloes       11.94             0.37            0.85         11.70        12.18
                                                    $z=2$                                   
                              $M_{\rm min}$   $\sigma_{\log M}$   $\alpha$   $M_{\rm cut}$   $M_{1}$
    10$ \% $ oldest haloes        11.15             0.11            1.08         11.66        12.33
          All haloes              11.41             0.30            0.92         11.60        12.13
   10$ \% $ youngest haloes       11.74             0.34            0.86         11.46        11.97
                                                    $z=3$                                   
                              $M_{\rm min}$   $\sigma_{\log M}$   $\alpha$   $M_{\rm cut}$   $M_{1}$
    10$ \% $ oldest haloes        10.89             0.13            1.01         11.47        12.07
          All haloes              11.17             0.32            0.93         11.34        11.91
   10$ \% $ youngest haloes       11.51             0.33            0.85         11.24        11.70

  : The HOD parameters described in Eqns. A1 to A4 for galaxy samples corresponding to a number density of $n=0.01\ h^{3}\, {\rm Mpc}^{-3}$ selected by stellar mass. From top to bottom we show the parameters for $z = 0,\ 1,\ 2\ \&\ 3$, for the $10$ per cent oldest haloes, for the full halo sample and for the $10$ per cent youngest haloes.[]{data-label="Table:HOD"}

\[lastpage\]

[^1]: <http://gavo.mpa-garching.mpg.de/Millennium/>
