On the prospect of using the maximum circular velocity of halos to encapsulate assembly bias in the galaxy-halo connection

TL;DR

This study examines whether replacing halo mass with maximum circular velocity ($ ext{V}_{max}$) in galaxy-halo models captures assembly bias effects, finding limited success and suggesting more complex properties like $ ext{V}_{peak}$ may be needed.

Contribution

The paper evaluates the effectiveness of using $ ext{V}_{max}$ instead of halo mass to incorporate assembly bias into galaxy-halo connection models, revealing its limitations.

Findings

Using $ ext{V}_{max}$ reduces central galaxy occupancy variation.

Satellite occupancy variation increases with $ ext{V}_{max}$.

Galaxy assembly bias remains largely unaffected by switching to $ ext{V}_{max}$.

Abstract

We investigate a conceptual modification of the halo occupation distribution approach, using the halos' present-day maximal circular velocity, $\vmax$, as an alternative to halo mass. In particular, using a semi-analytic galaxy formation model applied to the Millennium WMAP7 simulation, we explore the extent that switching to $\vmax$ as the primary halo property incorporates the effects of assembly bias into the formalism. We consider fixed number density galaxy samples ranked by stellar mass and examine the variations in the halo occupation functions with either halo concentration or formation time. We find that using $\vmax$ results in a significant reduction in the occupancy variation of the central galaxies, particularly for concentration. The satellites occupancy variation on the other hand increases in all cases. We find effectively no change in the halo clustering dependence on concentration, for fixed bins of $\vmax$ compared to fixed halo mass. Most crucially, we calculate the impact of assembly bias on galaxy clustering by comparing the amplitude of clustering to that of a shuffled galaxy sample, finding that the level of galaxy assembly bias remains largely unchanged. Our results suggest that while using $\vmax$ as a proxy for halo mass diminishes some of occupancy variations exhibited in the galaxy-halo relation, it is not able to encapsulate the effects of assembly bias potentially present in galaxy clustering. The use of other more complex halo properties, such as $\vpeak$, the peak value of $\vmax$ over the assembly history, provides some improvement and warrants further investigation.

Full text

1. Introduction

In the standard cosmological framework, galaxies form, evolve and reside in dark matter halos. It is of fundamental importance to understand the relation between galaxies and dark matter halos. How do the galaxies populate the halos? How do the properties of galaxies depend on halo mass or other characteristics of the halos? What role does the environment have in the halo occupation? These questions lie at the core of our understanding of galaxy formation. They are also crucial if we are to take full advantage of the next generation of galaxy surveys aimed at measuring galaxy clustering with unparalleled accuracy. The cosmological constraints from these data will no longer be dominated by statistical errors but rather the accuracy of our theoretical models. Understanding how galaxies relate to the underlying dark matter is thus essential for optimally using galaxies as a cosmological probe.

The formation and evolution of the dark matter halos is dominated by gravity and can be predicted accurately using high-resolution cosmological numerical simulations and analytic models. The formation of the galaxies and their relation to the dark matter halos, however, is more complex and depends on the detailed physical processes leading to the varied observed galaxy properties. Studying in detail the galaxy-halo connection is therefore of paramount importance (see Wechsler & Tinker 2018 for a review).

A powerful approach to explore galaxy formation within dark matter halos is semi-analytic modeling (SAM; e.g., Cole et al. 1994, 2000; Benson et al. 2003; Baugh 2006; Croton et al. 2006; Bower et al. 2006; Somerville et al. 2008). In such models, halos identified from high resolution $N$-body simulations are “populated” with galaxies using analytical prescriptions for the evolution of baryons through cosmic time. These models have been successful in reproducing many measured properties including the galaxy luminosity and stellar mass functions (see, e.g., Croton et al. 2006; Bower et al. 2006; Fontanot et al. 2009; Guo et al. 2011, 2013; Henriques et al. 2015; Lagos et al. 2018; Baugh et al. 2019).

Alternatively, hydrodynamic simulations follow the physical processes which govern the behavior of baryons by solving the fluid equations, while also modeling some of the unresolved processes with sub-grid prescriptions (see, e.g., Somerville & Davé 2015). Cosmological hydrodynamical simulations have started to play a major role in the study of galaxy formation and evolution. Two recent efforts, the Illustris Project (Vogelsberger et al. 2014; Pillepich et al. 2018) and the EAGLE simulation (Schaye et al. 2015), have set the state-of-the-art in hydrodynamical calculations. These ambitious simulations are still however significantly smaller than large-scale structure dark-matter only simulations and harder to fine-tune to the observations.

A useful approach to empirically connect galaxies with dark matter halos is the Halo Occupation Distribution (HOD) framework (e.g., Jing et al. 1998; Benson et al. 2000; Peacock & Smith 2000; Seljak 2000; Scoccimarro et al. 2001; Berlind & Weinberg 2002; Berlind et al. 2003; Cooray & Sheth 2002; Yang et al. 2003; Kravtsov et al. 2004; Zheng et al. 2005). The HOD formalism characterizes the relationship between galaxies and halos in terms of the probability distribution, ${\rm P(N|{\rm M_{vir}})}$, that a halo of virial mass ${\rm M_{vir}}$ contains ${\rm N}$ galaxies of a given type, together with the spatial and velocity distributions of galaxies inside halos. The key ingredient is the halo occupation function, ${\rm\langle N({\rm M_{vir}})\rangle}$, representing the average number of galaxies as a function of halo mass. The advantage of this approach is that it does not rely on assumptions about the (poorly understood) physical processes that drive galaxy formation and can be directly constrained from the observations.

When considering the halo occupation function it is often useful to separate the contribution of central galaxies, namely the main galaxy at the center of the halo, and that of the additional satellite galaxies that populate the halo (Kravtsov et al. 2004; Zheng et al. 2005; Jimenez et al. 2019). Standard applications assume a cosmology and a parametrized form for the halo occupation functions motivated by predictions of SAMs and hydrodynamic simulations (e.g., Zheng et al. 2005). The HOD parameters are then constrained using galaxy clustering measurements from large surveys, the galaxies abundance and the predicted halo clustering. The approach has been demonstrated to be a powerful theoretical tool to study the galaxy-halo connection, transforming clustering measurements into a physical relation between galaxies and dark matter halos. It has been successful in explaining the shape of the galaxy correlation function, its dependence on galaxy properties and environmental dependence (e.g., Zehavi et al. 2004, 2005, 2011; Berlind et al. 2005; Abbas & Sheth 2006; Skibba et al. 2006; Tinker et al. 2008; Coupon et al. 2012). It has also become an increasingly popular method to create realistic mock catalogs, by populating halos in large $N$-body simulations (e.g., Manera et al. 2015; Zheng & Guo 2016; Smith et al. 2017; DeRose et al. 2019), important for planning and analysis of current and upcoming surveys.

A central assumption in the standard applications of the HOD framework is that the galaxy content of halos depends only on the host halo mass. The origins of this assumption is in the Press-Schechter formalism (Press & Schechter 1974; Lacey & Cole 1993) and the uncorrelated nature of random walks describing halo assembly, resulting in a correlation of the halo environment with its mass but not with its assembly history (Bond et al. 1991; White 1999; Lemson & Kauffmann 1999). This assumption has been challenged by explicitly demonstrating in large $N$-body simulations that the clustering of halos varies with halo formation time, concentration, substructure occupation and spin, at fixed halo mass (Sheth & Tormen 2004; Gao et al. 2005; Gao & White 2007; Wechsler et al. 2006; Jing et al. 2007; Wetzel et al. 2007; Lazeyras et al. 2017; Sato-Polito et al. 2018), an effect that has been generally referred to as “halo assembly bias”.

To what extent is the galaxy distribution impacted by the assembly bias of their host halos is an actively debated topic. If galaxy properties closely correlate with halo formation history this would lead to a dependence of the galaxy content on large-scale environment and a corresponding change in the amplitude of galaxy clustering on large scales. The latter is commonly referred to as “galaxy assembly bias” (a misnomer referring to the manifestation of halo assembly bias in the galaxy distribution; GAB hereafter). The predictions for GAB have been explored with simulated galaxies (e.g., Zhu et al. 2006; Croton et al. 2007; Zu et al. 2008; Zentner et al. 2014; Chaves-Montero et al. 2016; Romano-Diaz et al. 2017; Xu et al. 2019), while the observational evidence for it remains inconclusive and controversial. Several suggestive detections have been put forward (Cooper et al. 2010; Wang et al. 2013; Lacerna et al. 2014b; Watson et al. 2015; Hearin et al. 2015; Miyatake et al. 2016; Montero-Dorta et al. 2017) while other studies indicate the impact of assembly bias is small (Abbas & Sheth 2006; Blanton & Berlind 2007; Tinker et al. 2008; Lacerna et al. 2014a; Lin et al. 2016; Zu & Mandelbaum 2016; Walsh & Tinker 2019) and that previous claimed detections are due to systematics (e.g., Campbell et al. 2015; Zu et al. 2017; Sin et al. 2017; Tinker et al. 2017; Sunayama & More 2019).

If significant, such GAB can have direct implications for interpreting galaxy clustering using the HOD framework (Pujol & Gaztanaga 2014; Zentner et al. 2014). The existence of GAB essentially implies a dependence of the halo occupation functions on the secondary halo parameters in addition to mass. This “occupancy variation” (Zehavi et al. 2018; Z18 hereafter) is essentially the crucial link between halo assembly bias and GAB; both halo assembly bias and occupancy variation are required to produce GAB.

Z18 explore the predicted occupancy variation in SAMs applied to the Millennium simulation, focusing on the dependence of the galaxy content of halos on large-scale environment and on halo formation time. They find distinct occupancy variations with central galaxies in denser regions preferentially occupying lower-mass halos. A similar, but significantly stronger, trend is found with halo age, where early-formed halos are more likely to host central galaxies at lower halo mass. A reverse trend is seen for the occupation of satellites, with early-forming halos having fewer satellites. Z18 investigate the stellar mass-halo mass relation for central galaxies in the SAM, finding a clear dependence on secondary halo properties, e.g., at fixed halo mass, older halos tend to host more massive galaxies, which drive the centrals occupancy variation. They also examine the resulting GAB, arising from the combined effect of halo assembly bias and the occupancy variations. Similar results for the occupancy variation have also been obtained in the hydrodynamical simulations EAGLE, Illustris and IllustrisTNG (Artale et al. 2018; Bose et al. 2019).

Contreras et al. (2019, hereafter C19) extended the work of Z18 to study the cosmic evolution of assembly bias and occupancy variation in the Guo et al. (2013) SAM. They explore the redshift range between 0 and 3, considering galaxy samples selected by either stellar mass or star formation rate (SFR), and selecting halos by both halo formation time and halo concentration. At the present epoch, both halo concentration and halo formation time produce similar features when examining either halo clustering or the occupancy variation. These two halo properties however exhibit different evolutionary scenarios. The GAB signature monotonically decreases when going to higher redshift (and for lower number density samples), reversing its sense in some instances.

Building on the work of Z18 and C19, we set out here to explore a conceptual modification of the HOD approach which shifts from halo mass to an alternative proxy, aimed at encapsulating assembly bias in the galaxy-halo connection. We also draw upon abundance matching techniques (e.g., Conroy et al. 2006; Reddick et al. 2013; Chaves-Montero et al. 2016; Lehmann et al. 2017). Such methods typically associate dark matter (sub)halos with galaxies using a monotonic relation between a galaxy property (such as stellar mass or luminosity) and a halo property like its infall halo mass, its maximum circular velocity, ${\rm V_{max}}$, or the peak value of ${\rm V_{max}}$ across the assembly of the halo, ${\rm V_{peak}}$. These are expected to have a tighter relation with the galaxy properties compared to the virial mass of the halo, ${\rm M_{vir}}$.

Given the distinct trends with concentration (or halo age) in the galaxy-halo connection and their crucial role in producing the occupancy variations (Z18), we choose ${\rm V_{max}}$ as the alternate halo property. ${\rm V_{max}}$, the maximal value of the circular velocity inferred from the halo mass distribution, effectively characterizes the gravitation potential. Using the Guo et al. (2013) SAM applied to the Millennium WMAP7 simulation, we explore the extent to which switching from virial mass to the halo’s maximum circular velocity in the HOD is able to encapsulate the main features of assembly bias. We consider all three aspects of this phenomenon: occupancy variation, GAB, and halo assembly bias. While some aspects are certainly improved, we find (spoiler alert) that, regrettably, ${\rm V_{max}}$ is unable to “capture” the essence of assembly bias and improve upon the use of ${\rm M_{vir}}$ in measuring GAB. Still, we think it is educational and useful to go through this exercise.

The overall motivation and potential applications of this work are many. First, we can obtain further insight into assembly bias and on how best to capture its essential features. Beyond the increased theoretical understanding of this complex phenomenon, had this conceptual change to the HOD successfully encapsulated assembly bias, it would have provided varied practical applications. The modified HOD could then have been utilized as an efficient tool for creating mock catalogs from simulations, which would reproduce realistic galaxy clustering while incorporating assembly bias. Moreover, it could have been potentially invaluable for inferring the galaxy-halo connection from galaxy clustering free of the effects of assembly bias.

The outline of the rest of the paper is as follows. Section 2 describes the galaxy samples we use. In Section 3 we present the occupancy variation with halo concentration and age when using ${\rm V_{max}}$ instead of halo mass. In Section 4 we examine the relation between stellar mass and the different halo properties. Section 5 presents our results for galaxy assembly bias, and we conclude in Section 6. Appendix A investigates halo assembly bias with ${\rm V_{max}}$. Appendix B includes additional results for SFR-selected samples and samples at higher redshift. Finally, Appendix C shows some results when alternatively using ${\rm V_{peak}}$.

2. The Galaxy Samples

2.1. Simulation and Galaxy Formation Model

We use the Millennium WMAP7 N-body simulation (Guo et al. 2013), which is similar to the original Millennium simulation (Springel et al. 2005) but with cosmological parameters consistent with the seven-year WMAP data (Komatsu et al. 2011). The simulation has a comoving boxsize of $500\,h^{-1}\,{\rm Mpc}$ on a side and follows $2160^{3}$ particles with a mass resolution of $9.36\times 10^{8}h^{-1}{\rm M}_{\odot}$. Multiple outputs of the simulation are available at different snapshots from redshift $50$ to the present day. At each snapshot, halos are identified using a friends-of-friends algorithm (Davis et al. 1985), and SUBFIND (Springel et al. 2001) is then run on these to identify subhalos. Halo merger trees are constructed by following the evolution of the halos and subhalos with time.

Semi-analytic modeling is a fundamental methodology to model the evolution of baryons over cosmic time in a cosmological context (see, e.g., Baugh 2006; Lacey et al. 2016). SAMs aim to model the main physical processes involved in galaxy formation and evolution, grafted onto a halo merger tree. Some of these processes include gas cooling, star formation, feedback from supernovae and active galactic nucleii, chemical enrichment, dark matter halo mergers and galaxy mergers. The SAM used in our work is that of Guo et al. (2013), which is a version of the L-GALAXIES SAM code of the Munich group, based on Guo et al. (2011) and applied to the Millennium WMAP7 simulation. It is publicly available from the Millennium Archive††6  http://gavo.mpa-garching.mpg.de/Millennium/.

For the main part of our analysis we use galaxy samples with different number densities, ranked by the stellar mass of the galaxies. We use three number densities, $0.0316,\,0.01,$ and $0.00316\,h^{3}\,{\rm Mpc}^{-3}$, corresponding to stellar mass thresholds of $\sim 1.6\times 10^{9},\,1.2\times 10^{10}$ and $3.4\times 10^{10}h^{-1}{\rm M}_{\odot}$ , respectively. The samples are approximately evenly spaced in logarithmic number density, and follow the choices made in Z18. Figure 1 shows the cumulative stellar mass function and the three number density cuts used. The galaxies selected in each case are those to the right of the intersection with the corresponding horizontal line. Pertinent results for galaxy samples selected by ranked SFR are shown in Appendix B.

2.2. Halo Properties

We consider as the mass of the halo the virial mass, namely the mass enclosed within the virial radius of the halo corresponding to a density of 200 times the critical density of the universe. This choice, which we denote here as ${\rm M_{vir}}$, is largely a matter of convention but has been shown to roughly correspond to the boundary at which halos are in approximate dynamical equilibrium (e.g., Cole & Lacey 1996). We hereafter interchangeably refer to ${\rm M_{vir}}$ as the halo mass (but cf. Jiang et al. 2014 for an alternate mass definition). The circular velocity profiles of the halos obtained from the dark matter distribution, ${\rm V_{c}^{2}=GM(r)/r}$, have a well defined maximum, which serves as a natural halo size scale denoted as ${\rm V_{max}}$. This quantity is robustly defined for simulated halos and effectively characterizes the depth of the gravitational potential (e.g., Bullock et al. 2001; Conroy et al. 2006; Diemand & Moore 2011). We consider here in detail the use of ${\rm V_{max}}$ as a proxy for halo mass.

For the secondary parameters, we use the halo concentration and halo formation time, following C19. These are regarded as two fundamental parameters related to the halo assembly often used in assembly bias studies (e.g., Gao et al. 2005; Wechsler et al. 2006; Gao & White 2007; Mao et al. 2018; C19; Bose et al. 2019). The halo concentration characterizes the density profile. It is canonically defined as ${\rm C_{vir}=r_{vir}/r_{s}}$, where ${\rm r_{vir}}$ is the virial radius of the halo and $r_{\rm s}$ is the inner transitional radius appearing in the Navarro, Frenk & White (1996) profile, at which the density profile changes slope. It is often alternatively defined as the ratio between ${\rm V_{max}}$ and ${\rm V_{vir}}$, where ${\rm V_{vir}}$ is the virial velocity of the halo, ${\rm V_{vir}}\equiv{\rm V_{c}(r_{vir})}$. We use the latter definition here, which is directly calculable from simulation data and does not require any model fitting, and can be applied to the lowest-mass halos considered, which might have too few particles to fit a density profile (Bullock et al. 2001; Gao & White 2007; Diemand & Moore 2011).

The formation time of the halo is defined as the redshift at which the main progenitor of the halo first reaches half of its current mass, denoted as ${\rm z}_{0.5}$. We calculate it from the halo merger trees of the simulation, linearly interpolating the halo mass among the time snapshots available. This definition has been very commonly used in assembly bias studies (e.g., Gao et al. 2005; Gao & White 2007; Croton et al. 2007; Z18; Han et al. 2019). As we use the halo formation time simply to rank the halos as early-forming or late-forming we do not anticipate any dependence on the specific definition (cf., the study of Li et al. 2008).

Figure 2 shows the relation between ${\rm V_{max}}$ and ${\rm M_{vir}}$ for the halos in the simulation, color coded either by halo concentration or formation time. We plot only 1% of the halos, randomly chosen and ordered, to avoid overcrowding. As expected, there is a tight relation between the two, with the scatter arising distinctly from the differences in concentration. For fixed halo mass, ${\rm V_{max}}$ is directly related to the concentration (by definition), with more concentrated halos corresponding to larger ${\rm V_{max}}$. A similar relation is noted when color coding the halos by formation time, albeit with some scatter in the dependence on halo age, which arises from the scatter between concentration and halo age.

3. Occupancy Variations

We now proceed to examine the halo occupation functions and their variations with halo concentration and age. For each of these properties, we rank the halos by the secondary property in narrow (0.1 dex) bins of halo mass and identify the 20% extremes of the distribution. This factors out the halo mass dependence on these parameters and allows us to examine the occupancy variations for halos of the same mass. We also do the same for fine (0.04 dex) bins in ${\rm V_{max}}$. We have verified for both cases that our results are insensitive to the exact binning choice. Given the tight relation between ${\rm M_{vir}}$ and ${\rm V_{max}}$ there is significant overlap between the identification binned by either mass proxy.

Figure 3 presents the halo occupation functions for the $n=0.01\,h^{3}\,{\rm Mpc}^{-3}$ galaxy sample and their variation with concentration. The left-hand side shows the “standard” occupation functions as a function of halo mass. We see that more concentrated halos start hosting central galaxies at lower halo mass, while they host on average fewer satellites per halo. These are in full agreement with previous results explored in detail by Z18 and C19. The uncertainties on the occupation functions, estimated from jackknife resampling, are negligible over the range of halo masses plotted here, and hence they are not included. (They only start becoming noticeable for masses larger than $\sim 10^{14}h^{-1}{\rm M}_{\odot}$ where there are very few halos; See, e.g., Fig. 3 of Z18 or Fig. 5 of C19).

The right-hand panel of Fig. 3 shows our new results for the halo occupations now using ${\rm V_{max}}$ as our proxy for halo mass. We see that the general shape of the HOD remains the same. However, there are significant changes to the occupancy variation of central galaxies and satellites. We find that the ${\rm V_{max}}$ occupancy variation for central galaxies is very nearly diminished in this case, with all halos, most concentrated ones and least concentrated halos exhibiting similar occupation by galaxies. The ${\rm V_{max}}$ occupation functions for satellites, however, exhibit larger differences (namely, increased occupancy variation).

We also examine the differences in the occupancy variation for the other galaxy samples in Figure 4, finding some dependence on the number density, or rather stellar mass threshold, of the sample. For the lower number density (more massive galaxies; left-hand side) switching from ${\rm M_{vir}}$ to ${\rm V_{max}}$ only partially removes the central occupancy variation, while it increases further the satellites occupancy variation. While for the higher number density (less massive galaxies; right-hand side) the central occupancy variation is nearly overly compensated by the switch to ${\rm V_{max}}$, resulting in the least concentrated halos being occupied by central galaxies at slightly lower ${\rm V_{max}}$ than the most concentrated halos.

This change in behavior is not really unexpected and simply stems from the dependence of the (standard) occupancy variation with number density, as exhibited in the top panels of Figure 4 and studied in detail in Z18. The halo mass occupancy variation for both centrals and satellite galaxies increases with stellar mass (i.e., with decreased number density). Switching from ${\rm M_{vir}}$ to ${\rm V_{max}}$ amounts to a roughly constant shift in the mapping from ${\rm M_{vir}}$ to ${\rm V_{max}}$, for the 20% least/most concentrated halos, as can be seen in Fig. 2. This results in the slight over-compensation, nearly full compensation and partial compensation of the ${\rm V_{max}}$ occupancy variation, respectively, for the three cases with decreasing number density. Thus the nearly diminished ${\rm V_{max}}$ occupancy variation seen in Fig. 3 is somewhat coincidental.

To get a better sense of the amount of occupancy variations captured by switching to ${\rm V_{max}}$, we explore the variations also with regard to another fundamental halo parameter, the halo formation time. Figure 5 shows the occupation functions as a function of ${\rm M_{vir}}$ and ${\rm V_{max}}$, for the three number densities explored, now focusing on the 20% early-formed halos and 20% late-forming halos. When examining the standard HOD as a function of ${\rm M_{vir}}$, we find the well-studied variations with age (Z18; C19), with older halos more likely to host central galaxies at lower mass and to have fewer satellites. These variations are reduced when switching to ${\rm V_{max}}$ (bottom panels of Fig. 5), for all number densities, but only partially so. This is in accord with the emerging understanding that galaxy assembly bias is not fully governed by a single parameter (or simple combination thereof; Croton et al. 2007; Xu et al. 2019; Bose et al. 2019). The reduction noted here is likely due to the general correlation between halo formation time and concentration, and hence ${\rm V_{max}}$. This correlation of halo age with ${\rm V_{max}}$ is also seen in the right panel of Fig. 2.

4. Relation between Galaxy and Halo Properties

The study of Z18 first showed that the essential features of the centrals occupancy variation stem from the dependence of the galaxies’ stellar mass on secondary halo properties at fixed halo mass. They examined the stellar mass-halo mass relation (see their Fig. 9), finding that at fixed halo mass, more massive centrals tend to reside in older or denser halos. This directly produces the changes in occupation we observe, and the level of scatter in these secondary trends controls the strength of the occupancy variation. Similar diagnostics were also utilized in other works such as Matthee et al. (2017); Tojeiro et al. (2017); Artale et al. (2018); Xu et al. (2019).

We revisit this relation in the context of our work. The left panel of Figure 6 examines the relation between stellar mass and host halo mass for the central galaxies in the SAM, and the galaxies are color coded by the halo concentration. In an analogous fashion, we find here that at fixed halo mass, the more concentrated halos tend to host more massive central galaxies. The three horizontal lines mark the stellar mass thresholds used to define our samples. By examining the population of centrals above each threshold, we can visually see how the occupation variations come about, including the galaxies in the most-concentrated halos at lower halo mass and vice versa. The extent of the spread in the horizontal direction also directly correlates with the magnitude of the variations, producing larger variations for the higher stellar mass thresholds (lower number densities), as seen in Fig. 3 and 4. We note, however, that the color-coding by concentration here is done globally, while when defining the samples for the occupation functions, the 20% extremes are defined as a function of halo mass, so one should be careful when making the comparison.

The right panel of Fig. 6 shows the relation between the stellar mass of the central galaxies and ${\rm V_{max}}$ of their host halos, once again color-coded by concentration. The stellar mass-${\rm V_{max}}$ relation has overall slightly less scatter, as expected (see Chaves-Montero et al. 2016; Matthee et al. 2017). In this case, we see a less-obvious trend with concentration. Examining the three cases separately (the three horizontal lines) we can intuitively understand the result of each. For the middle sample, there is a very mild trend with concentration. The top sample (lowest number density) has in fact a larger horizontal scatter and trend with concentration (after accounting for the different concentration markings noted above), causing the remaining occupancy variation shown in Fig. 4. Going to the lowest stellar-mass threshold we in fact note that the trend reverses with the lower-concentration halos appearing at lower halo mass than the higher concentration ones. Our highest number density sample is just at the cusp of this change.

5. Galaxy Assembly Bias

Having investigated the behavior of the halo occupation functions and the resulting occupancy variations with ${\rm V_{max}}$, we now proceed to see its impact on galaxy clustering, namely GAB. For completeness, we also provide a brief investigation of halo assembly bias in Appendix A, showing the concentration-dependent halo clustering with ${\rm M_{vir}}$ and ${\rm V_{max}}$.

To investigate the impact of assembly bias on galaxy clustering we compare the clustering of galaxies in our sample to that of shuffled galaxy samples. The shuffling follows the methodology of Croton et al. (2007) (see also Z18 and C19), randomly reassigning the galaxy content of halos among halos of the same mass. More specifically, the central galaxies are randomly shuffled among halos within the same mass bin. The satellite galaxies are moved together with their original central galaxy, preserving the same distribution, and thus maintaining the same contribution to the correlation function from intra-halo pairs. In order to investigate galaxy assembly bias with respect to ${\rm V_{max}}$, we perform an analogous shuffling procedure in ${\rm V_{max}}$ bins. We present the results using 0.05 dex bins in ${\rm M_{vir}}$ and 0.044 dex bins in ${\rm V_{max}}$, but we checked different binnings and verified that this makes no difference to our results.

The shufflings remove any dependence of the galaxy population on secondary properties of the halos (or their assembly history) other than that inherent in halo mass or ${\rm V_{max}}$. The difference between the clustering of the original galaxy sample and the shuffled sample, then reflects the impact of assembly bias on galaxy clustering, namely the level of GAB. In the idealized extreme case, if ${\rm V_{max}}$ were able to encapsulate all the occupancy variations exhibited by galaxies, there would be no difference between the clustering of the original sample and the ${\rm V_{max}}$-shuffled one.

Figure 7 shows the resulting GAB signatures when shuffling either by ${\rm M_{vir}}$ or ${\rm V_{max}}$, for the three stellar-mass selected samples studied. The uncertainties on this clustering ratio, estimated from jackknife resampling, are again negligible over most of the range (and start becoming noticeable only above separations of $20\,h^{-1}\,{\rm Mpc}$; see, e.g., Fig. 10 of Z18 or Fig. 10 of C19). The results for ${\rm M_{vir}}$ are identical to the ones examined in C19, showing a $\sim$12%-16% excess clustering due to assembly bias. This arises from the combined effect of the occupancy variation and halo assembly bias: central galaxies preferentially occupy the more concentrated halos which are more strongly clustered. The surprising, and rather disappointing, results lies with the ${\rm V_{max}}$ case. We find that in general the GAB associated with ${\rm V_{max}}$ is comparable to that of ${\rm M_{vir}}$. This is somewhat puzzling given the significant reduction of the centrals occupancy variation when using ${\rm V_{max}}$.

To further assess the situation, we repeat this analysis but now considering only the central galaxies, the results of which are presented in the top panel of Figure 8. Note that we only show here larger scales in the 2-halo regime, since without satellites there are no intra-halo pairs. We see that when considering just the central galaxies there is a significant reduction in the GAB signal, of about $40\%$ at $10\,h^{-1}\,{\rm Mpc}$ scale. This reduction originates from the amount of assembly bias “captured” by ${\rm V_{max}}$, reflected in the much reduced occupancy variation with concentration (shown in Fig. 3). The occupancy variations however are not fully captured by the switch to ${\rm V_{max}}$ when considering other parameters, as demonstrated for halo age in Fig. 5 (see also Croton et al. 2007; Matthee et al. 2017) and thus most of the GAB effect remains.

The lack of overall reduction in GAB for ${\rm V_{max}}$ compared to the standard GAB signature for ${\rm M_{vir}}$ exhibited in Fig. 7 is likely due to the significantly increased occupancy variation for satellites, which also contribute to the signal via central-satellite pairs (Zu et al. 2008). It appears that the decreased centrals occupancy variation for ${\rm V_{max}}$ together with the increased satellites occupancy variation, and coupled with the roughly unchanged halo assembly bias (Appendix A), result in a comparable level of GAB to the original measurement.

To remedy the increased contribution from the satellites occupancy variation we attempt a hybrid shuffling scheme, where the central galaxies are shuffled by ${\rm V_{max}}$ while the satellites are shuffled by ${\rm M_{vir}}$. The result of this hybrid shuffling is shown as the green line in the bottom panel of Fig. 8, exhibiting a roughly $15\%$ reduction in GAB on large scales. The slight deviation of the hybrid case from a ratio of unity on small scales arises from changes to the contribution from central-satellite pairs due to the more complex shuffling procedure. One can envision improving this further by finding a different halo property that would diminish the satellites occupancy variations, perhaps by utilizing the number of substructures. Alternatively, one may be able to develop a composite new parameter that will simultaneously improve both the centrals occupancy variation and the satellites one. However, attempts to do this at least in the context of halo assembly bias (e.g., Mao et al. 2018; Villarreal et al. 2017; Xu & Zheng 2018; Han et al. 2019) have shown that this is largely unattainable, and such modeling approaches might be too complex to be practical in any case.

Another potential candidate is ${\rm V_{peak}}$, the peak value of ${\rm V_{max}}$ throughout the accretion history of the halo. This parameter is often used in the context of abundance matching methods (e.g., Conroy et al. 2006; Reddick et al. 2013; Chaves-Montero et al. 2016; Guo et al. 2016) to better connect galaxies to halo substructure. Its main advantage in such methodologies, however, is in connecting satellites to subhalos, which differs from the HOD approach relating both central and satellite galaxies to the main host halo. Still, recent claims have suggested that ${\rm V_{peak}}$ provides the tightest relation to the stellar mass of galaxies (He 2019) and is free from secondary dependences (Xu et al. 2019), touting its use as a better proxy for halo mass. Remaining scatter in the relation is attributed to stochastic baryonic effects (Matthee et al. 2017; Kulier et al. 2019).

Figure 9 shows the resulting GAB, when using ${\rm V_{peak}}$ as the primary halo property instead of ${\rm M_{vir}}$ or ${\rm V_{max}}$, for the three different galaxy samples. Namely, we calculate the clustering of the galaxies relative to a galaxy sample shuffled according to ${\rm V_{peak}}$, where ${\rm V_{peak}}$ is obtained from the merger tree of each halo. We see that using ${\rm V_{peak}}$ partially suppresses the GAB signature, mostly in the 1-halo to 2-halo transition region. This effect is more significant with increased number density, as more small galaxies are included. Appendix C shows further analysis of the relation of ${\rm V_{peak}}$ to other properties, suggesting that the differences arise primarily from low-mass, high-concentration and older halos where ${\rm V_{peak}}$ varies from ${\rm V_{max}}$. This may be related to the existence of splashback galaxies, i.e., galaxies whose halos had their mass accretion histories truncated due to a close encounter with a larger halo. Such halos are thought to have a significant role in halo assembly bias of low mass halos (Mansfield & Kravtsov 2019), and are likely to impact the intermediate “transition” scales, as these centrals likely reside just outside the virial radius of the larger halo. We defer further investigation of this phenomenon with higher-resolution simulations to future work.

We note, however, that even if ${\rm V_{peak}}$ can be utilized as a parameter which partially encapsulates assembly bias, determining it in simulations requires calculating the full merger trees of the halos with sufficient mass resolution. In such a scenario, one will likely also have access to other sophisticated modeling techniques which may be more useful. We focused in this work on a property like ${\rm V_{max}}$ that would be readily available in most N-body simulations commonly used for creating large mock catalogs and for constraining cosmology.

6. Conclusion

We use a state-of-the-art semi-analytic galaxy formation model applied to the Millennium simulation to study the prospects of a conceptual modification of the HOD approach, replacing the virial mass of the halo by its maximal circular velocity of the halo, ${\rm V_{max}}$. The motivation is that this revised halo occupation function may encapsulate the effects of assembly bias into the formalism, enabling more accurate modeling of the galaxy-halo connection and galaxy clustering, and allow us to produce realistic mock catalogs that incorporate this effect. We thus explore the different aspects of assembly bias, namely halo clustering dependence on secondary parameters (halo assembly bias; see Appendix A), the variation in the galaxy content of halos with these parameters (occupancy variation; § 3) and the their impact on galaxy clustering relative to a shuffled galaxy sample (galaxy assembly bias; § 5). We mostly use here galaxy samples with fixed number density ranked by stellar mass at the present epoch, and for the secondary halo parameters we investigate the variation with halo concentration and halo formation time. To get a broader understanding of the origins of the occupancy variation, we also examine the relation between stellar mass and the different halo properties (§ 4). Finally, we also investigate the potential of utilizing ${\rm V_{peak}}$, the peak value of ${\rm V_{max}}$ across the halo’s assembly history as the proxy for halo mass.

The main conclusions from our work are summarized as follows:

• •

Halo assembly bias, i.e. the dependence of halo clustering on concentration, is largely unchanged when examined in bins of fixed ${\rm V_{max}}$ versus fixed halo mass. The same holds for ${\rm V_{peak}}$.

• •

Employing ${\rm V_{max}}$ significantly reduces the occupancy variation with halo concentration for central galaxies, however, it increases the satellites occupancy variation.

• •

The centrals occupancy variation is partially reduced when using halo formation time as the secondary halo property.

• •

The change to ${\rm V_{max}}$ does not reduce the level of GAB, despite the reduction in the central occupancy variation. The GAB signature remains essentially the same when using ${\rm V_{max}}$ or ${\rm M_{vir}}$, irrespective of sample number density.

• •

${\rm V_{max}}$ does not prove to be a useful quantity also when examining samples selected by SFR or when looking at $z=1$, with varying results.

• •

Using ${\rm V_{peak}}$ slightly reduces the GAB signal, impacting mostly the 1-halo to 2-halo transition regime, an effect likely related to splashback galaxies. However, the large-scale effect remains largely unmitigated.

Assembly bias remains a challenge for contemporary models of galaxy clustering and the galaxy-halo relation. Perhaps a more intricate parameter can better encapsulate the effects of galaxy assembly bias, however given the complex nature of the phenomenon this might be hard to achieve and not-trivial to implement in practice. While the level of assembly bias in the Universe remains an open and debated issue, the results shown here can help inform theoretical modeling of it and attempts to determine it in observations.

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. We acknowledge useful discussions with participants of the Shanghai Assembly Bias workshop in the final stages of this work. IZ, SC, EJ & NP acknowledge the hospitality of the ICC at Durham University. IZ and SEK acknowledge support by NSF grant AST-1612085. SC is supported by a Juan de la Cierva Formacion Fellowship (FJCI-2017-33816). EJ acknowledges support from “Centro de Astronomía y Tecnologías Afines” BASAL 170002. NP acknowledges support from Fondecyt Regular 1191813. This project received financial support from the European Union’s Horizon 2020 Research and Innovation programme under the Marie Sklodowska-Curie grant agreement number 734374.

Appendix A A. Halo Assembly Bias

We examine the concentration dependence of halo clustering in the simulation, namely halo assembly bias. Following Gao et al. (2005), we bin the halos in discrete bins of halo mass and calculate the auto-correlation function of these halo samples and of the 20% most concentrated and 20% least concentrated halos. These are shown in the left-hand side of Figure 10, and exhibit the well-studied concentration-dependent halo clustering (e.g., Wechsler et al. 2006; Gao & White 2007; Mao et al. 2018). We see that for relatively-low halo masses, more concentrated halos are more clustered than less concentrated halos. This trend is the strongest for our lowest mass bin, decreases with increasing mass, then reverses sense and continues to increase in amplitude for the most massive halos.

The right-hand side of Figure 10 repeats this analysis, but now in fixed bins of ${\rm V_{max}}$. The bins are chosen to roughly match the halo mass bins according to the relation between ${\rm V_{max}}$ and ${\rm M_{vir}}$ (Fig. 2) and have identical number densities. We find very similar results for the concentration-dependent halo clustering as a function of ${\rm V_{max}}$, with comparable assembly bias amplitudes in all bins and a reversal of the trend at larger ${\rm V_{max}}$ values. Our results are in good agreement with Sato-Polito et al. (2018) who find a similar behavior with a reversal of the halo assembly bias effect at ${\rm M_{vir}}\sim 10^{13}h^{-1}{\rm M}_{\odot}$ and at ${\rm V_{max}}\sim 330{\rm kms}^{-1}$, using the MultiDark suite of simulations (Klypin et al. 2016). For completeness, we also investigate the concentration-dependence of clustering for fixed bins of ${\rm V_{peak}}$, the peak value of ${\rm V_{max}}$ across the halos’ assembly history (not shown here). We find nearly identical results for ${\rm V_{peak}}$ as for ${\rm V_{max}}$. We conclude that switching from ${\rm M_{vir}}$ to either ${\rm V_{max}}$ or ${\rm V_{peak}}$ has no significant impact on halo assembly bias.

Appendix B B. Results for SFR Selected Samples and for Redshift $z=1$

We report here an analogous analysis to the one in the main part of the paper but performed for galaxy samples selected by their SFR, which may be relevant for galaxy selections of upcoming surveys. We use the same three number densities. Figure 11 shows the occupancy variation for one representative case of the $n=0.01\,h^{3}\,{\rm Mpc}^{-3}$ number density sample. The left-hand side shows the standard occupation functions as a function of halo mass. This is the same case presented in Fig. 7 of C19. Of note is the characteristic shape of the occupation function, which is different than that of stellar mass selected samples, due to the paucity of star-forming galaxies residing as centrals in massive halos. The occupancy variation for SFR-selected samples is similar to that of stellar mass selected ones, with more concentrated halos preferentially hosting central galaxies at the “knee” of the occupation and having fewer satellites. In the “dip” of the centrals occupation at higher mass, a smaller fraction of the more concentrated halos tends to host star-forming central galaxies.

The right panel of Fig. 11 shows the occupancy variation in the ${\rm V_{max}}$ case. In this case, the occupancy variation for the central galaxies at the “turnover” reverses sense, with low-concentration halos starting to host central galaxies at lower halo mass. The variations in the central occupation at higher halo masses remain similar, but perhaps slightly reduced, while the satellite occupancy variations increase.

The net impact on galaxy clustering is shown in Figure 12, where we show again the ratio of correlation function to that of shuffled samples where the galaxy content was randomly reassigned to halos of the same ${\rm M_{vir}}$ or ${\rm V_{max}}$, which effectively erases the occupancy variations. The overall changes with number density are consistent with those found by C19. We find that for all SFR-selected samples, this ratio slightly decreases for ${\rm V_{max}}$ with respect to the one for ${\rm M_{vir}}$, likely due to the preferential occupation of centrals in low-concentration halos. Hence we find that, for all cases, the correlation function of the ${\rm V_{max}}$ shuffled samples is larger than that of the ${\rm M_{vir}}$ shuffled one. The resulting impact on GAB, however, changes with number density (i.e., SFR threshold). For the highest number density sample, utilizing ${\rm V_{max}}$ slightly decreases GAB (defined as the deviation from a ratio of unity); for the middle number density it nearly diminishes the GAB signal; while for the lowest number density it increases the GAB effect. Thus it is hard to draw any conclusions on the usefulness of switching to ${\rm V_{max}}$.

Finally, in the left-hand side of Figure 13 we explore the GAB signatures at a higher redshift of $z=1$. At $z=1$ we find that the clustering ratio is now larger for ${\rm V_{max}}$ than for ${\rm M_{vir}}$ for all cases, or rather that the clustering of the ${\rm V_{max}}$-shuffled sample is lower than that of the ${\rm M_{vir}}$-shuffled sample. And, once again, the specific impact on GAB depends on the number density. For completeness, we also show in the right-hand side of Fig. 13 the galaxy assembly bias results for the stellar mass selected samples at $z=1$. We see that in this case as well, for all number densities, the clustering with respect to the ${\rm V_{max}}$-shuffled sample is increased relative to the ${\rm M_{vir}}$ case, resulting in increased GAB signatures. Overall, the results presented here strengthen our conclusion that, despite its claimed “potential”, ${\rm V_{max}}$ is unable to encapsulate galaxy assembly bias effects.

Appendix C C. Additional results with ${\rm V_{peak}}$

In this appendix we provide supplementary information regarding using ${\rm V_{peak}}$, the maximum value of ${\rm V_{max}}$ for each halo over cosmic history, as our proxy for halo mass. Figure 14 presents the relation of ${\rm V_{peak}}$ to ${\rm M_{vir}}$ color coded either by halo concentration (left) or formation time (right). These relations can be compared to the analogous ones for ${\rm V_{max}}$ shown in Fig. 2. We find that the relations are similar, but with a larger scatter in the ${\rm V_{peak}}$ case, especially for the low ${\rm M_{vir}}$ / low ${\rm V_{peak}}$ range. This arises from the scatter between ${\rm V_{peak}}$ and ${\rm V_{max}}$, primarily for low-mass halos where ${\rm V_{peak}}$ varies from ${\rm V_{max}}$ and extends to larger values. It is noteworthy that the low-mass halos with the largest ${\rm V_{peak}}$ values tend to have higher concentrations and, strikingly, earlier formation times. This supports our hypothesis that these are splashback halos that had a larger mass in the past.

Figure 15 shows the relation between stellar mass and ${\rm V_{peak}}$ for the SAM central galaxies. This relation is similar to the analogous one for ${\rm V_{max}}$ shown in Fig. 6. However, it exhibits less secondary dependences on concentration – note the lack of noticeable extremes of the most concentrated halos – in particular for lower values of stellar mass. This contributes to the slight improvement in GAB seen in Fig. 9 and its dependence on number density.

Finally, Figure 16 shows the occupancy variations with concentration for an HOD computed as a function of ${\rm V_{peak}}$. Again, we find very similar results to the occupancy variations with ${\rm V_{max}}$ (Fig. 3 and 4). The centrals occupancy variation improves (i.e., decreases) to a varying degree for the different samples while the satellites occupancy variation increases with respect to the “standard” occupancy variations with ${\rm M_{vir}}$. There are slight differences in the level of occupancy variation for ${\rm V_{peak}}$ with respect to those for ${\rm V_{max}}$, with the highest number density sample (lowest stellar mass threshold) exhibiting now nearly diminished centrals occupancy variation, likely correlated to the relatively bigger improvement seen for the GAB measurement in that case (bottom panel of Fig. 9).