The Three Causes of Low-Mass Assembly Bias

TL;DR

This paper investigates the physical causes of low-mass halo assembly bias, identifying three main processes: large-scale tidal effects, gravitational heating, and splashback subhaloes, through detailed analysis and elimination of other factors.

Contribution

It provides a comprehensive analysis of the physical mechanisms behind low-mass halo assembly bias, highlighting three key contributing processes and ruling out others.

Findings

Splashback subhaloes account for two thirds of the bias.

Remaining bias is due to a small fraction of haloes in dense regions.

Large-scale tidal fields and gravitational heating are primary causes.

Abstract

We present a detailed analysis of the physical processes that cause halo assembly bias -- the dependence of halo clustering on proxies of halo formation time. We focus on the origin of assembly bias in the mass range corresponding to the hosts of typical galaxies and use halo concentration as our chief proxy of halo formation time. We also repeat our key analyses across a broad range of halo masses and for alternative formation time definitions. We show that splashback subhaloes are responsible for two thirds of the assembly bias signal, but do not account for the entire effect. After splashback subhaloes have been removed, we find that the remaining assembly bias signal is due to a relatively small fraction ($\lesssim 10\%$) of haloes in dense regions. We test a number of additional physical processes thought to contribute to assembly bias and demonstrate that the two key processes are the slowing of mass growth by large-scale tidal fields and by the high velocities of ambient matter in sheets and filaments. We also rule out several other proposed physical causes of halo assembly bias. Based on our results, we argue that there are three processes that contribute to assembly bias of low-mass halos: large-scale tidal fields, gravitational heating due to the collapse of large-scale structures, and splashback subhaloes located outside the virial radius.

Tables

Table 1. Table 1: The fraction of haloes in the Bolshoi simulation which are removed by the different cuts described in the text. The first three rows show the subhalo fraction, f subhalo , subscript 𝑓 subhalo f_{\rm subhalo}, for the different subhalo cuts described in section 2.3 . The last six rows correspond to the assembly-bias-removing cuts described in section 3.2 for different proxies and show f removed , subscript 𝑓 removed f_{\rm removed}, the fraction of haloes which must be removed after splashback subhaloes have been cut from the sample, and ρ S , subscript 𝜌 𝑆 \rho_{S}, the Spearman correlation coefficient between c vir subscript 𝑐 vir c_{\rm vir} and a given proxy. Note that f removed subscript 𝑓 removed f_{\rm removed} and ρ S ​ ( c vir , X ) subscript 𝜌 𝑆 subscript 𝑐 vir 𝑋 \rho_{S}(c_{\rm vir},X) are completely uncorrelated, as discussed in section 2.8.1 .

Subhalo definition	$f_{\mathrm{subhalo}}$		Section
$R_{\mathrm{vir}}$ subhaloes	0.27		§2.3
flyby subhaloes	0.33		§2.3.1
splashback subhaloes	0.37		§2.3.2
Removal criterion	$f_{\mathrm{removed}}$	${\rho }_S$	Section
	$0.48\pm 0.04$	-0.47	§2.6
	$0.16\pm 0.03$	-0.16	§2.5.1
$N_5>18\pm 2$	$0.14\pm 0.02$	0.09	§3.2
	$0.09\pm 0.01$	-0.23	§2.5.2
	$0.09\pm 0.01$	-0.19	§2.5.2
	$0.056\pm 0.007$	-0.36	§2.6

Table 2. Table 2: The same as Table 1 , but for the BolshoiP simulation

Subhalo definition	$f_{\mathrm{subhalo}}$		Section
$R_{\mathrm{vir}}$ subhaloes	0.28		§2.3
flyby subhaloes	0.33		§2.3.1
splashback subhaloes	0.38		§2.3.2
Removal criterion	$f_{\mathrm{removed}}$	${\rho }_S$	Section
	$0.54\pm 0.07$	-0.47	§2.6
	$0.18\pm 0.03$	-0.18	§2.5.1
$N_5>21\pm 2$	$0.16\pm 0.03$	0.09	§3.2
	$0.11\pm 0.02$	-0.24	§2.5.2
	$0.11\pm 0.01$	-0.20	§2.5.2
	$0.062\pm 0.009$	-0.36	§2.6

Equations

$$M_{\Delta }=\frac{4\pi }{3}\Delta {\bar{\rho }}_{\mathrm{m}}{R}_{\Delta }^3,$$

$$\Delta =\frac{18{\pi }^2+82\tilde{x}-39{\tilde{x}}^2}{{\Omega }_{\mathrm{M}}},$$

$$\rho (r)=\frac{{\rho }_s}{x(1+x{)}^2},$$

$$D_{\mathrm{vir},i}$$

$$=$$

$$\begin{array}{rl}{T}_{ij}=\sum _k& \frac{m_k}{(x_k^2+y_k^2+z_k^2{)}^{5/2}}\times \\ & \left(\begin{array}{ccc}{y}_k^2+z_k^2-2x_k^2& -3x_k{y}_k& -3x_k{z}_k\\ -3x_k{y}_k& x_k^2+z_k^2-2y_k^2& -3y_k{z}_k\\ -3x_k{z}_k& -3y_k{z}_k& x_k^2+y_k^2-2z_k^2\end{array}\right).\end{array}$$

$${\Phi }_{\mathrm{steepest}}(r_1)=\frac{1}{2}{\alpha }_1{r}_1^2,$$

$$R_{\mathrm{tidal}}={\left(\frac{G{M}_{\mathrm{vir}}}{{\alpha }_1}\right)}^{1/3};M_{\mathrm{tidal}}=M(<R_{\mathrm{tidal}})$$

$$M_{\mathrm{tidal},\mathrm{b}}={\int }_0^{R_{\mathrm{tidal}}}dR{\int }_0^{V_{\mathrm{esc}}(R)}dV\frac{dM}{dRdV}$$

$$V_{\mathrm{esc}}=V_{\mathrm{vir}}{\left\{2\frac{(1+c_{\mathrm{vir}})\ln (1+c_{\mathrm{vir}}x)}{x\left[(1+c_{\mathrm{vir}})\ln (1+c_{\mathrm{vir}})-c_{\mathrm{vir}}\right]}\right\}}^{1/2}.$$

$$M_{\beta ,\mathrm{b}}={\int }_0^{\beta R_{\mathrm{vir}}}dR{\int }_0^{V_{\mathrm{esc}}(R)}dV\frac{dM}{dRdV}$$

$$\mathcal{M}(r)=\frac{⟨m_i{m}_j{⟩}_{i,j\in P(r)}-⟨m{⟩}^2}{⟨m^2⟩-⟨m{⟩}^2}.$$

$$M_{\mathrm{vir}}(a_{1/2})$$

$$R_{\mathrm{tidal}}^{(2)}={\left(\frac{3}{2}\right)}^{1/3}{R}_{\mathrm{Hill}}^{(2)},$$

$$F_{\mathrm{tidal}}^{(\mathrm{ex})}(r)\propto {\left(\frac{r}{R_{\mathrm{tidal}}^{(2)}}\right)}^{1+\alpha },$$

$$\frac{R_{\mathrm{tidal}}^{(\mathrm{ex})}}{R_{\mathrm{tidal}}^{(2)}}={\left(\frac{F_{\mathrm{tidal}}^{(\mathrm{ex})}(L_i^{(2)})}{F_{\mathrm{tidal}}^{(2)}(L_i^{(2)})}\right)}^{-\frac{1}{3+\alpha }},$$

$$|R_{\mathrm{tidal}}^{(\mathrm{ex})}/R_{\mathrm{tidal}}^{(2)}-1|<|F_{\mathrm{tidal}}^{(\mathrm{ex})}/F_{\mathrm{tidal}}^{(2)}-1{|}^{-1/2}.$$

Full text

††pagerange: The Three Causes of Low-Mass Assembly Bias–D††pubyear: 2018

1 Introduction

The most visually striking feature of the large-scale structure of the universe is the clustered, web-like distribution of galaxies, with vast voids separated by walls and filaments (e.g., Bond et al., 1996). Understanding the clustering of galaxies within the context of the $\Lambda+$Cold Dark Matter ($\Lambda$CDM) model relies on the generic model in which galaxies are formed by the dissipation of diffuse baryon plasma within growing dark matter haloes (e.g., White & Rees, 1978). Galaxy clustering is then interpreted in terms of the clustering of dark matter haloes (e.g., see Desjacques et al., 2018; Wechsler & Tinker, 2018, for recent reviews), which is generally different from that of matter, i.e. the distribution of haloes is “biased,” relative to the mass distribution (Kaiser, 1984).

Halo bias depends primarily on halo mass (e.g., Mo & White, 1996; Sheth & Tormen, 1999) and this dependence is now both well-understood theoretically and well-calibrated numerically (Desjacques et al., 2018). It is also now known that halo bias has secondary dependences on other halo properties, such as formation time, concentration, spin, and ellipticity (Gao et al., 2005; Wechsler et al., 2006; Harker et al., 2006; Gao & White, 2007; Jing et al., 2007; Li et al., 2008; Faltenbacher & White, 2010; Villarreal et al., 2017; Sato-Polito et al., 2018; Han et al., 2018). The first such secondary dependence was found for halo formation time and its closely related proxy – halo concentration (Gao et al., 2005; Wechsler et al., 2006; Harker et al., 2006; Jing et al., 2007) and has become known as “assembly bias.” Specifically, the bias of “old” haloes (early formation time) is generally different than that of “young” (late formation time) haloes, with the difference depending on halo mass and the definition of formation time(e.g., Li et al., 2008).

Assembly bias is important for the theoretical interpretation of galaxy clustering and its potential to provide useful cosmological constraints (e.g., Abazajian et al., 2005). There have been significant observational efforts to detect the related but distinct phenomenon of galaxy assembly bias, the dependence of galaxy clustering on secondary halo properties, which themselves experience halo assembly bias (Zentner et al., 2014; Wechsler & Tinker, 2018). This is a difficult task, because the effect needs to be unambiguously disentangled from the dependence of galaxy properties on halo mass and satellite classification scheme, meaning that observational signatures of galaxy assembly bias have as yet proved to be elusive (Campbell et al., 2015; Lin et al., 2016). If halo assembly bias does have a signature in galaxy clustering, it would be important to understand its physical origin in order to construct robust and accurate models. Conversely, if halo assembly bias does not have observational signatures, it would be important to understand why tracers of halo age and tracers of galaxy age behave differently. We note, however, that this paper concerns itself exclusively with halo assembly bias and not with galaxy assembly bias.

The focus of this paper is to understand the physical origin of halo assembly bias, particularly in the regime of galaxy-scale halo masses. This is distinct from the origin of assembly bias at large masses, which is related to the properties of the peaks of the initial Gaussian density perturbations from which these massive haloes collapse (Dalal et al., 2008). Peaks with the same mass but different curvature will cluster differently because peaks with larger curvatures are located in lower-density environments, while peaks with smaller curvatures are in higher-density regions. This gives rise to assembly bias because peak curvature is directly related to a halo’s mass accretion history, which is also affected by tidal torques from the surrounding anisotropic mass distribution (Desjacques, 2008). Although this curvature-related bias can be reduced by compensating effects found in some proxies of halo age Zentner (2007); Sandvik et al. (2007); Mao et al. (2018), it is present for other age definitions and when more physical definitions of halo boundaries and masses are used (Chue et al., 2018).

At smaller halo masses, however, the physics of assembly bias is more complex because the mass evolution of haloes is determined by a combination of the properties of their initial density peaks, and also by non-linear processes (e.g., Wang et al., 2007; Hahn et al., 2009). The simple and striking manifestation of this is that the sign of assembly bias switches for small-mass haloes when $c_{\rm vir}$ is used as a measure of halo age (Wechsler et al., 2006; Dalal et al., 2008).

A number of studies have explored the physical processes that can give rise to halo assembly bias in the small-mass regime. One readily apparent process is the non-linear effects that a massive host halo can exert on its smaller-mass neighbours. In particular, “splashback” (often also called “backsplash”) subhaloes pass within the inner regions of a larger halo but are located outside its virial radius at the epoch of analysis. Such haloes appear isolated, but will have had their mass accretion histories truncated due to their previous close encounters with their hosts and have thus been studied as a potential source of low-mass assembly bias (Wang et al., 2009; Li et al., 2013; Wetzel et al., 2014; Sunayama et al., 2016).

Although splashback subhaloes are mostly found within three virial radii of their host halo, they can give rise to an assembly bias signal at much larger distances. This is because at large scales the spatial distribution of splashback subhaloes will track the distribution of their massive hosts and will therefore be more strongly clustered than that of distinct haloes. A similar effect would occur if subhaloes located within the virial radius of their host were included in the sample used to measure halo clustering and assembly bias. This is illustrated in Fig. 1, which compares the clustering of early- and late-forming haloes with splashback subhaloes included and removed, respectively. Removing splashback subhaloes significantly reduces the difference in clustering between the two halo samples, even on scales much larger than the virial radius of the most massive haloes within the volume. Nevertheless, multiple studies have demonstrated that splashback subhaloes alone cannot be responsible for the entire assembly bias signal (Wang et al., 2009; Sunayama et al., 2016), a fact that can be seen visually in Fig. 1. A similar conclusion was reached by Hearin et al. (2015), albeit in the related but distinct context of galactic conformity.

Another process that could contribute to assembly bias is the truncation of a halo’s mass growth by the tidal force generated its most gravitationally-dominant neighbor (Hahn et al., 2009; Behroozi et al., 2014; Hearin et al., 2016b; Salcedo et al., 2018). Even though single-halo tidal forces become small beyond $\approx 3\--5\times R_{\rm vir}$ of the host, the haloes truncated by these forces can give rise to large-scale assembly bias in a way similar to splashback subhaloes. A similar truncation of halo mass growth can be caused by the overall tidal force from all of the surrounding haloes and structures in the matter distribution (Hahn et al., 2009; Wang et al., 2011; Paranjape et al., 2018; Musso et al., 2018), as the largest filaments and sheets generate strong tidal fields throughout their volumes. Since these structures can be several tens of Mpc in size, they can comfortably give rise to assembly bias on large scales. This effect has been characterized in terms of both the tidal force and the anisotropy of the tidal field, although, in practice, a high degree of anisotropy tends to correlate with the magnitude of the tidal force, so it is not clear that the two effects can be separated cleanly.

Finally, the gravitational heating of matter within large-scale structure structure has been proposed as a process that can contribute to assembly bias (Wang et al., 2007; Dalal et al., 2008). For example, matter within the deep potentials of filaments can acquire large velocities during accretion, and thus cannot be accreted by small-mass haloes located within the filament. The smaller accretion rates of such haloes would thus give rise to assembly bias. Note that although gravitational heating and strong tidal forces generally happen in similar regions, they are physically distinct phenomena: tidal forces arise via large second derivatives in the gravitational potential, while gravitational heating is caused by the potential depth.

Although significant effort has been devoted to studying these effects (Wang et al., 2007; Dalal et al., 2008; Wang et al., 2009; Hahn et al., 2009; Wang et al., 2011; Li et al., 2013; Wetzel et al., 2014; Sunayama et al., 2016; Hearin et al., 2016b; Paranjape et al., 2018; Salcedo et al., 2018; Musso et al., 2018), their relative importance and a coherent physical picture for the origin of low-mass halo assembly bias has not yet been established. The primary goal of this paper is to rectify this. To this end, we define a set of quantitative proxies for each of the different processes outlined above and use them to investigate the relative contribution of these processes to the low-mass assembly bias signal. Specifically, we examine how efficiently sample cuts defined by each proxy can remove the signal.

The paper is organized as follows. In section 2 we describe basic definitions and measurements and describe our cosmological simulations and halo sample, with sections 2.5-2.6 focusing on the proxies of the processes described above, and section 2.8 describing the core methodology of this paper. In section 3 we present measurements and estimates of the relative contribution of different processes to low-mass assembly bias. We discuss topics related to the interpretation of this work in section 4 and summarize our results in 5. The key results of this study are presented in Fig. 3.

2 Methods

2.1 Simulations and codes

In this paper we use halo catalogues and particle snapshots from the Bolshoi and BolshoiP cosmological dark matter-only $N$-body simulations, made available through the CosmoSim111https://www.cosmosim.org and Skies and Universes222http://skiesanduniverses.iaa.es databases (Klypin et al., 2017). Both simulations followed the evolution of 20483 particles in boxes of size 250 $h^{-1}$ Mpc using the ART code, with force resolutions of $\approx 1\,h^{-1}$ kpc. The Bolshoi simulation assumed $\Lambda$CDM cosmology with parameters consistent with the WMAP 9 year constraints: $\Omega_{\rm M}=0.270,$ $\Omega_{\rm B}=0.047,$ $\sigma_{8}=0.820,$ $n_{s}=0.95,$ and $H_{0}=70.0\ {\rm km\ s^{-1}\,Mpc^{-1}}$. The BolshoiP simulations assumed parameters consistent with the constraints from the Planck observatory: $\Omega_{\rm M}=0.307,$ $\Omega_{\rm B}=0.048,$ $\sigma_{8}=0.823,$ $n_{s}=0.96,$ and $H_{0}=67.8\ {\rm km\ s^{-1}\,Mpc^{-1}}$. The corresponding particle masses are $1.3\times 10^{8}\ h^{-1}\,M_{\odot}$ and $1.5\times 10^{8}\ h^{-1}\,M_{\odot}$, respectively. Additional technical information, such as time stepping criteria and initial conditions, can be found in Klypin et al. (2011) and Klypin et al. (2016). Although all plots shown in this paper use Bolshoi data, all analyses were repeated for BolshoiP with similar results. We also use the simulation suite described in Diemer & Kravtsov (2015) for some convergence and performance testing.

Haloes in the Bolshoi and BolshoiP simulations were identified using version 0.99RC2+ of the Rockstar halo finder (Behroozi et al., 2013c), and version 1.0+ of the related consistent-trees method (Behroozi et al., 2013b) was used to construct halo merger trees. The catalogues and merger trees we use were downloaded from the CosmoSim database. We use the Shellfish algorithm to identify splashback shells – the 3D surfaces formed by the outermost apocentres of accreted matter (Mansfield et al., 2017)333https://github.com/phil-mansfield/shellfish. We use the Colossus python package (Diemer, 2017a)444http://www.benediktdiemer.com/code/colossus/ to calculate various relevant cosmological quantities and statistics and the halotools package (Hearin et al., 2016a)555https://halotools.readthedocs.io to calculate correlation functions efficiently .

2.2 Basic halo properties

Throughout this paper, we use halo masses, $M_{\Delta},$ defined within a radius, $R_{\Delta}$, enclosing a specified density contrast $\Delta=\rho(<R_{\Delta})/\bar{\rho}_{\rm m}$:

[TABLE]

where $\bar{\rho}_{\rm m}$ is the mean density of matter. We use the “virial” density contrast defined as

[TABLE]

where $\widetilde{x}=\Omega_{\rm M}(z)-1$. This definition formally corresponds to the virialization density of a spherical tophat perturbation collapsing in a flat $\Lambda$CDM universe at the redshift of analysis (Bryan & Norman, 1998). This choice is dictated by the fact that Rockstar calculates properties of haloes using particles contained within friends-of-friends (FOF) groups identified using the linking length of $b=0.28$ (in units of mean interparticle separation) in the Bolshoi and BolshoiP halo catalogues. Such linking length ensures full percolation only out to density contrasts $\gtrsim\Delta_{\rm vir}$ (Behroozi et al., 2013c), which means that not all particles are included in the groups at distances corresponding to lower $\Delta$. This, in turn, biases halo properies such as concentration. Additionally, density profiles of haloes begin to deviate from the NFW (Navarro–Frenk–White) form beyond $R_{\rm vir}$ (e.g., Becker & Kravtsov, 2011; Diemer & Kravtsov, 2014), which means that NFW profile fits at $R\gtrsim R_{\rm vir}$ are generally biased (see Appendix A for an extended discussion on this issue).

We note, however, that $R_{\rm vir}$ does not correspond to an actual physical boundary of the halo (for detailed discussion see, e.g., More et al., 2015 as well as section 2.3 below). As such, we would have been free to select many other equally valid choices of $\Delta$ for our baseline definition, so long as they are cross-matched with profile properties measured from a catalogue which uses $\Delta\geq\Delta_{\rm vir}$. We discuss the significance of this choice when interpreting our results in section 4.2 and find that it is not a significant issue.

We assume that the radial density profiles of dark matter haloes are well-approximated by an NFW profile within $R_{\rm vir}$:

[TABLE]

where $x\equiv r/R_{\rm s}$ (Navarro et al., 1997). This profile is fully specified by a combination of $R_{\rm vir}$ and the “concentration,” $c_{\rm vir}=R_{\rm vir}/R_{\rm s}$. We use the $c_{\rm vir}$ values measured by the Rockstar code, which calculates these values by breaking particles within $R_{\rm vir}$ into radial bins containing at least 15 particles, calculating densities in each bin, and performing a $\chi^{2}$-minimization fit to Eq. 3

We also make use of halo circular velocity profiles, $V(r)=\sqrt{GM(<r)/r}$, and its maximum values, $V_{\rm max}$. The largest value of $V_{\rm max}$ during the evolution of a halo is denoted as $V_{\rm peak}$. $V_{\rm max}$ and $V_{\rm peak}$ are useful quantities because they allow us to define halo samples by potential depth without any dependence on our choice of halo boundary. Additionally, some halo samples will contain subhaloes, and the high ambient density around these objects can mean that $R_{\rm vir}$ is not a meaningful property. $V_{\rm peak}$ is useful in particular because it is a proxy for halo potential depth prior to mass loss or tidal stripping (e.g., Nagai & Kravtsov, 2005), physical processes which both play significant roles in our analysis. Furthermore, models of galaxy clustering indicate that “peak” mass definitions are better correlated with the observable properties of galaxies than values at the current epoch (Reddick et al., 2013).

There are a number of definitions of halo age used in the literature: single-epoch accretion rates (e.g., Lacey & Cole, 1993; Li et al., 2008), current halo properties - such as concentration - related to a halo’s mass accretion history (e.g., Wechsler et al., 2006; Villarreal et al., 2017; Sato-Polito et al., 2018), the epoch at which a halo first achieved half of its current mass (e.g., Gao et al., 2005), or a characteristic timescale of an analytic fit to halo mass accretion history (e.g., Wechsler et al., 2002; Zentner, 2007). In this paper, we primarily adopt $c_{\rm vir}$ as a tracer of halo age, with older haloes having larger concentrations. We briefly explore the effect of using different definitions in section 4.2.

We focus on $c_{\rm vir}$ for several reasons. First, $c_{\rm vir}$ has been demonstrated to strongly correlate with a number of explicit indicators of halo age (Bullock et al., 2001; Wechsler et al., 2002; Zhao et al., 2003; Lu et al., 2006; Ludlow et al., 2013, 2014). Second, the connection between accretion history and $c_{\rm vir}$ has a solid theoretical underpinning (Zhao et al., 2003; Lu et al., 2006; Dalal et al., 2010), as demonstrated by the accuracy of the concentration models based on halo mass accretion history (e.g. Zhao et al., 2003, 2009; Dalal et al., 2010; Ludlow et al., 2014; Diemer & Joyce, 2018). Third, the convergence criteria for halo density profiles (e.g. Power et al., 2003; Navarro et al., 2004; Springel et al., 2008), and for concentrations (see, e.g., Section 3.2 in Diemer & Kravtsov, 2015) are well studied and it is thus relatively straightforward to identify regimes in which numerical concentrations can be trusted.

2.3 Definition of halo boundaries and subhaloes

Throughout this paper, we define subhaloes as the haloes located within the boundary of a larger “host halo,” and refer to all non-subhaloes as “distinct haloes.” Of course, this classification depends on the definition of halo boundary and will have a clear qualitative meaning only if we use halo boundary definition that corresponds to an actual physical boundary.

Traditionally, spheres of radius $R_{\rm vir}$ (or some other overdensity radius) are used as halo boundaries, but this choice has a number of issues (see, e.g., Diemer et al., 2013b; More et al., 2015). The first issue is that that there is no commonly-used overdensity radius that corresponds to any physical change or feature in the radial profiles of various halo properties (see, e.g., fig. 3 of Diemer et al., 2013a). The second issue is that many studies have established that a substantial fraction of bound subhaloes and matter have first orbits whose apocentres take objects out to as far as $\approx 2-3\times R_{\rm vir}$ of the host halo (Gill et al., 2005; Ludlow et al., 2012; Mansfield et al., 2017; Diemer, 2017b).

Fortunately, haloes do have unambiguous edges manifested as sharp drops in density and caused by the pileup of particles at the apocentres of their first orbits. These edges in form 3D surfaces called “splashback shells,” and enclose almost all matter and subhaloes ever accreted by a halo. Haloes outside $R_{\rm vir}$ of their host, but within its splashback shell are called “splashback subhaloes.”666The terminology used to refer to these objects is varied: different authors refer to them as “backsplash suhaloes” or “splashback subhaloes,” and often refer to them as “haloes” instead of “subhaloes.” All these terms refer to the same concept. Some authors may use the term “flyby [sub]haloes” interchangeably with “splashback [sub]haloes,” although the former term generally implies that merger tree analysis has been used.

Splashback subhaloes can be identified and removed in one of two ways. The first is a classification based on the past halo trajectories, where merger trees are used to determine whether a halo has ever been within a larger host (e.g., Ludlow et al., 2009; Wang et al., 2009; Diemer, 2017b; Diemer et al., 2017). The second is to directly identify splashback shells of haloes and flag all haloes wthin them as subhaloes. We adopt the second approach as our fiducial classification method, but employ both throughout the paper to ensure that our results are robust and do not rely on the specifics of either approach.

For lexical clarity, we refer to subhaloes identified through merger tree analysis as “flyby subhaloes” and subhaloes identified through the construction of splashback shells as “splashback subhaloes.”

2.3.1 Flyby subhaloes

To identify flyby subhaloes, we use the following procedure for each halo in the $z=0$ halo catalogue. First, using consistent-trees (Behroozi et al., 2013b), we identify the main-line branch for the halo, labelling the $z=0$ halo the “root halo” and all other haloes on the branch its “progenitor haloes.” If any haloes on the branch are within the virial radius of another halo at any redshift, the root is flagged as a flyby subhalo.

This process is complicated by the fact that during major mergers the virial radii of both merging haloes fluctuate significantly and it is common for both host haloes to be at least temporarily identified as subhaloes of one another. This can lead to the final host halo being misidentified as a flyby subhalo of an object that no longer exists once the merger is complete. To rectify this, if the search of a root halo’s progenitors reveals that some progenitor, $P,$ is within $R_{\rm vir}$ of a host halo, $H,$ we only classify the root halo of $P$ as a flyby subhalo when the following three conditions are met:

$H$ must have a root halo at $z=0$. 2. 2.

The root halo of $H$ must not be within $R_{\rm vir}$ of the root halo of $P.$ 3. 3.

The root halo of $H$ must have a strictly larger mass than the root halo of $P.$

Our tests indicate that just enforcing conditions 1 and 2 is sufficient to correct the overwhelming majority of false classifications. This procedure can be extended to root redshifts other than $z=0$.

Although the identification of flyby subhaloes is well-defined and only requires the use of a merger tree, it is not without drawbacks. First, the method uses $R_{\rm vir}$, which as we discussed above does not correspond to a physical halo boundary. Second, this approach does not distinguish between ordinary subhaloes with apocentres outside $R_{\rm vir}$, and subhaloes that may have undergone dynamical three body interactions that resulted in their unbinding and ejection and are a qualitatively distinct population from splashback subhaloes. Although a substantial fraction of subhaloes may have undergone such interactions (Sales et al., 2007; Ludlow et al., 2009), we find that haloes which have been ejected from the splashback shell are rare and do not have an impact on our analysis (see section 3.5). Third, this method does not count haloes within the splashback shell on their first infall as subhaloes, even though this population is similar to first-infall haloes within $R_{\rm vir},$ which this method does classify as subhaloes.

2.3.2 Splashback shell subhaloes

The simplest way to estimate the size of a halo’s splashback shell it to approximate it as a sphere and estimate its radius from the location of sharp steepening it causes in the halo’s density and subhalo number density profiles (e.g. Fillmore & Goldreich, 1984; Bertschinger, 1985; Diemer & Kravtsov, 2014; Adhikari et al., 2014; More et al., 2015; Diemer, 2017b; Diemer et al., 2017). This radius is then called the “splashback radius,” $R_{\rm sp}.$ However, the application of this method for individual haloes is not straightforward (see Mansfield et al., 2017). In addition, actual splashback shells are not spherical and spherical approximation may result in misclassification of a certain fraction of subhaloes. For this reason, we use the Shellfish algorithm (Mansfield et al., 2017) to identify fully 3D splashback shells.

The Shellfish algorithm identifies splashback shells by measuring sharp density drops in many 1D density profiles along tens of thousands of lines of sight around a halo and fits a flexible smooth 3D surface to their location (Mansfield et al., 2017). Once Shellfish has identified splashback shells, we use the efficient intersection-checking method described in Appendix B to flag all haloes within the splashback shell of any larger halo as splashback subhaloes.

There are three complications to using Shellfish which must be addressed before it can be used to construct subhalo catalogs: its $N_{\rm 200m}$ convergence limit, the occurrence of rare but catastrophic fitting failures, and its behavior for low-accreting hosts. We perform extensive tests on all three issues and find that once accounted for in the ways described below, they do not have a significant effect on our results.

First, Shellfish has a rather stringent convergence limit and requires that haloes have more than $5\times 10^{4}$ particles within $R_{\rm 200m}$, the overdensity radius corresponding to $\Delta=200\,\rho_{\rm m}$, to achieve $R_{\rm sp}$ measurements with accuracy better than $5\%$. This corresponds to the $M_{\rm 200m}\gtrsim 7-8\times 10^{12}\ h^{-1}M_{\odot}$ or $V_{\rm peak}\gtrsim 280\ {\rm km\ s^{-1}}$ in the Bolshoi and BolshoiP boxes. Below this mass, we use the fitting formula for the median $R_{\rm sp}$ provided in Mansfield et al. (2017), and flag haloes within spheres of radius $R_{\rm sp}$ instead. Tests using the higher resolution L0063_CBol box from Diemer & Kravtsov (2014) indicate that this results in a negligible number of subhalo misclassifications compared to using real Shellfish-identified splashback shells because the majority of the splashback subhaloes in our mass range have hosts larger than 280 km s*-1*.

The second complication is that for a small number of host haloes ($\approx 1\%$, Mansfield et al., 2017), irregularities in the local density field cause Shellfish to fail to identify the correct surface shape, adopting a barbell-shaped surface instead, which can cause subhaloes well within $R_{\rm vir}$ to be misclassified as distinct haloes. To mitigate this, we mark haloes as splashback subhaloes if they fall within either their host’s splashback shell or within a sphere centred on that host of radius $R_{\rm vir}.$ We analysed the distribution of the minimum radii of Shellfish shells in haloes which were visually-identified to be unaffected by this surface fitting failure and found that the minimum radii are generally larger than $R_{\rm vir}$. Thus, the procedure we adopt is unlikely to result in misclassification of distinct host haloes as subhaloes.

The third complication is that the Shellfish algorithm underestimates the size of splashback shells for haloes that are accreting slower than the baseline pseudo-evolution accretion rate (Mansfield et al., 2017). However, this only lowers the splashback radius by $\approx 10\%$ and few haloes massive enough to host subhaloes in our target mass range accrete this slowly, so it is not expected to be a significant issue. Empirically, we find that virtually all flyby splashback subhaloes whose hosts are in this accretion regime are also within the splashback shells of their hosts (see section 3.5), so we do not explicitly account for this effect.

2.4 Halo sample

Although we will examine the mass-dependence of assembly bias in section 3.4, the majority of our analysis focuses specifically on low-mass haloes. Our primary concern when defining a halo sample is to prevent the inclusion of haloes whose convergence radii are large enough that they introduce numerical effects into Rockstar’s measurements of $c_{\rm vir}.$ As mentioned above, the numerical reliability of density profiles has been well studied, but for cosmological simulations with small softening scales the exact convergence properties are covariant with particle count, softening scale, halo mass, and time stepping scheme (Power et al., 2003; Ludlow et al., 2018), so determining convergence limits for an individual simulation should always be done through the comparison of carefully constructed multi-box suites.

Because there is only a single Bolshoi box, we place an upper bound on the convergence limit using the CBol simulation suite described in Diemer & Kravtsov (2015). Of particular note is the box CBol_L0125, which has the same particle mass to Bolshoi, but which has much larger timesteps within halo centers, implying that the convergence radius of Bolshoi should be smaller than that of CBol_L0125.777The difference in softening scale between these boxes makes an exact comparison difficult without a detailed analysis beyond the scope of this paper. See Diemer & Joyce (2018) for some additional discussion on the subtleties of comparing Bolshoi to this simulation suite. We find that when using the same Rockstar version and configuration variables as our Bolshoi catalogue, the $V_{\rm peak}-c_{\rm vir}$ relation for the CBol_L0125 box agrees with the higher resolution CBol_L0063 box above $V_{\rm peak}=120$ km s${}^{-1},$ corresponding to a somewhat conservative cutoff particle count of $N_{\rm peak}\approx 1.3\times 10^{3}$.

Our low-mass halo sample includes haloes with $120\ {\rm km\ s^{-1}}<V_{\rm peak}<220\ {\rm km\ s^{-1}}$ (approximately $1.7\times 10^{11}\ h^{-1}\ M_{\odot}<M_{\rm peak}<1.2\times 10^{12}\ h^{-1}\ M_{\odot}$). Due to the slope of the halo mass function, the majority of haloes will be close to the lower mass limit, making the choice in upper mass limit less important. We chose the upper mass limit so our sample spans roughly a factor of eight in $M_{\rm peak}$ and find that our results are not particularly sensitive to this choice.

2.5 Measuring tidal force strength

Tidal forces have been proposed as a potential cause of assembly bias (Hahn et al., 2009; Wang et al., 2009; Hearin et al., 2016b; Salcedo et al., 2018; Paranjape et al., 2018) because they can slow down, stop, or reverse mass accretion. These fields are strongest in dense environments, such as within large-scale filaments or near the outskirts of massive haloes, allowing distant haloes in similar environments to have correlated accretion histories. Below, we describe methods for measuring the strength of both the single-halo tidal field, and the large-scale tidal field.

2.5.1 Tidal force from a single halo

A typical simplifying assumption when calculating the tidal force felt by a halo is to assume that it is primarily caused by a single massive halo. If one also assumes that the point of interest is orbiting around that halo on a circular orbit, one can compute the tidal Hill radius, $R_{\rm Hill},$ corresponding to the distance to the nearest two Lagrangian points when the effective potential is approximated to second order. However, the assumptions that are made in calculating $R_{\rm Hill}$ are not correct for distinct haloes in a $\Lambda$CDM cosmology. This is because these haloes are almost never on circular orbits around each other and, as we discuss in section 2.5.2, the tidal force generally has a significant contribution from multiple haloes and from the large-scale matter distribution. Thus, formally, the Hill radius is not a physically meaningful quantity for distinct haloes. Nevertheless, the classical Hill radius can be used to estimate the tidal force of a halo’s most gravitationally-dominant neighbour.

As a simple and definitionally robust proxy for $R_{\rm Hill}$ we use the virial radius-scaled distance, $D_{{\rm vir},i}$, for every distinct halo $i$:

[TABLE]

where $j$ runs over all distinct haloes within some search radius, $R_{0},$ which are more massive than the halo, and $R_{ij}$ is the distance between haloes $i$ and $j$. Haloes with smaller $D_{\rm vir}$ experience larger tidal forces and haloes with larger $D_{\rm vir}$ have smaller tidal forces. As Eq. 5 shows, $D_{\rm vir}$ is proportional to $R_{\rm Hill}$. This means that a rank-ordering by $D_{\rm vir}$ is equivalent to a rank-ordering by $R_{\rm hill}/R_{\rm vir}$, while formally $D_{\rm vir}$ is always a well-defined quantity and also allows for easy comparison with other assembly bias studies (e.g., Villarreal et al., 2017; Salcedo et al., 2018)

Our tests indicate that $D_{\rm vir}$ is well-converged for haloes in the mass range $120\ {\rm km\ s}^{-1}<V_{\rm peak}<220\ {\rm km\ s}^{-1}$ for $R_{0}\approx 100\,R_{\rm vir}.$

2.5.2 Large-scale tidal radius and mass

Although the single-source approximation is reasonably accurate for subhaloes, our tests indicate that most distinct haloes have multiple neighbors which contribute significantly to the tidal forces they feel. Moreover, we found that large-scale structures in mass distribution, such as filaments can contribute to the tidal force experienced by haloes substantially. For example, by combining the assumption of cylindrical symmetry with the radial density profiles of filaments reported in Cautun et al. (2014), we construct a toy model for filament potentials. Applying this model, we find that even in moderate-sized filaments with $R_{\rm filament}\gtrsim 3\,h^{-1}\,{\rm Mpc}$, the tidal force generated by the filament is comparable to or stronger than the typical tidal force generated by a halo’s single most gravitationally dominant neighbor.

For this reason, we compute the tidal radius of a halo calculated from the overall matter distribution around a halo, $R_{\rm tidal}$, as a proxy for the combined tidal force from all neighbour haloes and structures. To compute $R_{\rm tidal}$, we first construct the tidal tensor, $T_{ij}$, the Hessian of the external potential:

[TABLE]

Here, $k$ runs over all particles between two search radii, $R_{\rm min}$ and $R_{\rm max},$ $m_{k}$ is the mass of particle $k,$ and $x_{k},$ $y_{k},$ and $z_{k}$ are the components of the displacement vector from the halo centre to particle $k.$ The tidal radius lies along the steepest repulsive axis of the tidal field, and since the tidal tensor, like all Hessians, equivalently describes the second derivatives at the origin of a paraboloid with eigenvectors pointing along the paraboloid’s axes, the tidal field along the chief repulsive axis is given by

[TABLE]

where $\alpha_{1}$ is the most negative eigenvalue of $T_{ij},$ and $r_{1}$ is the radial distance along the corresponding eigenvector. We then assume that all non-tidal pseudo-forces (most notably the centrifugal force) are small and that at large distances the halo’s mass is well-approximated by $M_{\rm vir}$, making the tidal radius and the corresponding tidal mass

[TABLE]

To increase computational efficiency, we make two further approximations. First, we do not add the tidal contribution from any particles further than $100\,R_{\rm vir},$ and second, we subsample particles by a factor of 64 and multiply $m_{k}$ by 64 in Equation 6. Our tests indicate that the combined effects of both these approximations on $R_{\rm tidal}$ are at the sub-percent level. We set $R_{\rm min}=10\,R_{\rm vir}.$ This choice is discussed in detail in Appendix C.

Some authors have suggested that the primary feature of interest in the tidal field is its anisotropy, which can be defined in a number of ways (Wang et al., 2011; Paranjape et al., 2018). We chose to use $R_{\rm tidal}$ as a proxy for the total tidal force for two reasons. First, there are a number of different proxies for anisotropy and it is not clear a priori which definition is optimal. Second, we carried out analysis of assembly bias described in sections 2.7 and 2.2 using $\alpha_{R}$ and $q_{R}$ from Paranjape et al. (2018) and $t$ from Wang et al. (2011) as proxies for the tidal anisotropy and found that all of these proxies were not as efficient at removing assembly bias as $R_{\rm tidal}$.

2.6 Measuring gravitational heating

To gauge the contribution of gravitational heating to assembly bias, we use the mass of bound matter within the tidal radius, $R_{\rm tidal}$, defined in the previous section:

[TABLE]

Here, $V$ is the absolute velocity of a particle relative to the halo centre, while $V_{\rm esc}(R)$ is the escape velocity at a radius $R$ from the halo centre computed assuming that the halo is well-approximated by an NFW profile:

[TABLE]

Here, $x=r/R_{\rm vir},$ $V_{\rm vir}=\sqrt{GM_{\rm vir}/R_{\rm vir}},$ and $c_{\rm vir}$ is halo concentration. To speed up particle containment checks when computing mass profiles, we apply the algorithm described in Appendix B.

We also construct the variable

[TABLE]

for some constant $\beta.$ $M_{\beta,{\rm b}}$ allows us to isolate the effect of gravitational heating from the effect of external tidal fields because it does not include a dependence on $R_{\rm tidal}.$ Although a range of $\beta$ were used in our analysis, our results are primarily reported in terms of $\beta=3,$ for reasons we describe in section 3.2.

While these approximations are standard practice for computing particle boundedness, it is likely that they break down significantly in the outskirts of haloes. We discuss this in greater depth in Appendix D and argue that this should not have a significant effect on our results in section 4.1.

2.7 Assembly bias statistics

To study assembly bias, one must have a statistic that measures how clustering strength depends on a halo age proxy, $c_{\rm vir}$ in our case. The most direct approach is to split haloes into high-$c_{\rm vir}$ and low-$c_{\rm vir}$ samples, measure the clustering strength of each sample independently using correlation functions, and compare them. There are multiple ways of doing this, ranging from measuring the two-point correlation function of haloes, $\xi_{\rm hh}$, in each $c_{\rm vir}$-selected subsample to measuring the bias function, $b(r)=\xi_{\rm hm}/\xi_{\rm mm}$ (e.g, Gao et al., 2005; Gao & White, 2007; Faltenbacher & White, 2010). While this family of approaches is a valid and commonly-used, there are a number of associated issues. First, the definition of subsamples is arbitrary, and the strength of the measured signal depends on this definition somewhat. Second, if small $c_{\rm vir}$ ranges are chosen to maximize signal strength, statistical errors increase due to the comparatively small number of haloes used.

We use an alternative statistic – the marked correlation function (the MCF, Beisbart & Kerscher, 2000; Gottlöber et al., 2002) – which avoids this issue and which has been used in a number of assembly bias studies (e.g., Wechsler et al., 2006; Villarreal et al., 2017). For a sample of objects with assigned mark, $m$, the MCF is computed as:

[TABLE]

Here, $P(r)$ is the set of all pairs which are separated by a distance within the same radial bin as $r$. Following Villarreal et al. (2017), we define concentration marks for haloes in narrow circular velocity bins as their percentile within the $c_{\rm vir}$ distribution of that bin. Specifically, we use ten logarithmic bins in $V_{\rm peak}$ from 120 km s*-1* to 220 km s*-1*. The narrow bin width is required because the $c_{\rm vir}$ distribution is mass-dependent. This, combined with the mass-dependence of clustering, would result in illusory assembly bias signals in any halo sample defined over a sufficiently large mass range.

2.8 Measuring the connection between assembly bias and other variables

To evaluate the relative contribution of different physical processes to assembly bias, we need a way to gauge how strongly proxies for these processes, such as $D_{\rm vir},$ $R_{\rm tidal},$ $M_{\rm\beta,b},$ $M_{\rm tidal},$ or $M_{\rm tidal,b},$ are related to assembly bias. One simple way to do this is to measure the correlation coefficient between $c_{\rm vir}$ and each variable. However, as discussed in section 2.8.1, any approach that relies on measuring the connection between a proxy and formation time has serious issues.

Instead, in this paper, we follow an approach similar to that of Villarreal et al. (2017). We determine the strength of the connection between assembly bias and a proxy $X$ by finding the percentage of haloes ranked by $X$ that need to be removed from the sample to eliminate the assembly bias signal. For example, if 30% of haloes must be removed according to $X$ before the assembly bias signal is eliminated, but only 5% of haloes must be removed to achieve this for another proxy, $Y,$ we conclude that the physical process traced by $Y$ has a more significant contribution to assembly bias than the process traced by $X$.

Specifically, we first sort distinct haloes according to a proxy $X$, then remove a fraction of haloes $f=N_{\rm removed}/N_{\rm tot}$ for a series of $f$ values ranging from $0.01$ to $N_{\rm distinct}/N_{\rm tot}$ in steps of $0.01$. We then define $f_{\rm removed}$ as the minimum $f$ for which the MCF is within $1-\sigma$ of zero. The sample variance of the MCF is estimated by dividing the simulation box into eight octants, computing the MCF in each octant and finding the standard error on these MCFs at a constant $f_{\rm removed}$. Note that $f_{\rm removed}$ is normalized by the total number of haloes and not by the number of distinct haloes to make it easier to combine with different subhalo classification schemes.

We use a similar method to estimate the sample variance of $f_{\rm removed}$ itself, computing $f_{\rm removed}$ for each octant independently and finding the standard error on these values. Note that these errors on $f_{\rm removed}$ account for contributions from sample variance computed using the same octants, which means that while the uncertainties accurately estimate the scatter on measurements in independent boxes, there is likely covariance between the $f_{\rm removed}$ errors measured for different proxies within the same simulation. This means that the uncertainty on the relative ordering of $f_{\rm removed}$ values for multiple proxies within a single simulation is likely to be smaller than these errors would estimate. We discuss this further in section 3.2.

When calculating $f_{\rm removed},$ we compute the MCF in the radial range $[4,8]$ comoving $h^{-1}\,$Mpc. We have repeated all analysis in this paper with several other choices of radial ranges and did not find any significant qualitative difference in results. The primary result of moving to larger radii is that the amplitude of the reference MCF becomes smaller relative to the error, meaning that smaller cutoffs are able to make the signal consistent with zero. Thus, to be conservative, we use a relatively small-radius cutoff. We illustrate this in Fig. 2, which shows the MCF after distinct haloes below the $f_{\rm removed}$ cutoff for $M_{\rm tidal,b}/M_{\rm vir}$ have been removed from the sample: the MCF is consistent with zero out to $18\ h^{-1}$ Mpc. We have repeated all analysis in this paper using several different radial ranges and results remain qualitatively similar.

We note that this method is effective only for assembly bias models in which haloes are initially unbiased or negatively biased but where a small subset of haloes in extreme environments are pushed to older ages by some non-linear process. If, instead, assembly bias is strongly present in all environments, there will be no value of $f_{\rm removed}$ which can remove it. It is known that assembly bias is present across all halo ages (e.g., see fig. 3 in Wechsler et al., 2006), so a finding that there are variables with small values of $f_{\rm removed}$ would already put interesting constraints on the physics of assembly bias. We discuss this in more depth in section 3.3.

2.8.1 Difference between $f_{\rm removed}$ and age correlation

A number of previous studies evaluated the contribution of a given physical process with an associated proxy, $X,$ by measuring the correlation between $X$ and a proxy of halo age, $A$ (e.g., Hahn et al., 2009; Wang et al., 2011; Hearin et al., 2016b; Salcedo et al., 2018). This can be done using the Spearman’s rank coefficient, $\rho_{S}(A,X),$ or by measuring the slope of the average trend $X(A).$ While this approach provides indications of which proxies correlate well with halo age, by itself it cannot be used to gauge the relative contribution of different physical processes to assembly bias. This is because a correlation between age and proxy can only lead to assembly bias if clustering strength also varies strongly as a function of $X$. Comparison of the proxy–halo age correlation strength thus does not provide enough information to unambiguously gauge the contribution of the corresponding process to assembly bias. For example, $D_{\rm vir}$ and $R_{\rm tidal}/R_{\rm vir}$ have roughly the same level of correlation with $c_{\rm vir},$ but haloes experience wildly different differential clustering with with respect to both variables. Consequently, assembly bias is not connected to these two variables with the same strength.

As an illustration, Table 1 lists values of $f_{\rm removed}$ and the Spearman’s rank correlation coefficient, $\rho_{S},$ between $c_{\rm vir}$ and several different proxies and shows that these two quantities are almost completely unrelated. We therefore strongly recommend against drawing conclusions about assembly bias from measurements of correlation with halo age (see Mao et al., 2018, for additional discussion and caveats associated with using correlation coefficients in the context of assembly bias).

3 Analysis

3.1 Splashback subhaloes and assembly bias

We first test whether splashback subhaloes misclassified as distinct haloes by standard subhalo definitions (i.e., splashback subhaloes outside $R_{\rm vir}$ of a larger host) are responsible for low-mass halo assembly bias. The number of haloes removed by our different subhalo definitions is shown in Table 1. Our results are shown in Figure 2, using both methods discussed in section 2.7 for measuring assembly bias. The figure shows that splashback subhaloes cannot account for the entirety of assembly bias, although they contribute about two thirds of the signal. This is consistent with conclusions of the previous studies (Wang et al., 2009; Sunayama et al., 2016). The novel feature of this analysis is that we find a similar effect for two independent definitions of the splashback haloes: using evolutionary trajectories (§2.3.1) and using non-spherical 3D splashback shells identified using the Shellfish code (§2.3.2).

Note also that although results for the two definitions are similar, the two samples of haloes are not identical. Flyby subhaloes identified using merger trees are guaranteed to have passed their orbital pericentre and thus likely have experienced a strong tidal interaction with the host. On the other hand, when we use Shellfish all subhaloes within the splashback shell are classified as splashback subhaloes, including those haloes that have entered the shell, but have not yet passed their pericentre. Given that both halo samples have exactly the same local environments, the fact that removal of infalling splashback subhaloes results only in a small decrease of the assembly bias signal means that this portion of the assembly bias signal is due to the stripping subhaloes experience during their pericentre passage. Conversely, any mass growth suppression subhaloes experience on their way to pericentre is comparatively unimportant important.

We further compare the splashback subhaloes and flyby subhaloes in section 3.5.

3.2 Contribution of tidal truncation and gravitational heating to assembly bias

We now investigate how the truncation of halo mass growth by the tidal forces, both from a halo’s most gravitationally-dominant neighbour and from the entire large-scale matter distribution, contributes to assembly bias. We also investigate the contribution of dynamical heating caused by the collapse of matter into sheets and filaments. To this end we use the five proxies of these processes defined in sections 2.5-2.6 – $D_{\rm vir},$ $R_{\rm tidal},$ $M_{\rm\beta,b},$ $M_{\rm tidal},$ and $M_{\rm tidal,b}$ – and evaluate what fraction of the distinct halo sample ranked by each of the proxies must be removed to eliminate the assembly bias signal.

$D_{\rm vir}$ is the $R_{\rm vir}$-normalized distance to the most tidally dominant halo. It is a proxy of the one-halo contribution to the tidal force proportional to the traditional Hill radius. $R_{\rm tidal}$ is the tidal radius calculated using only the distant matter distribution and $M_{\rm tidal}$ is the mass contained within the tidal radius. $M_{\beta,{\rm b}}$ is the bound mass within $\beta\,R_{\rm vir}$ for a specified constant $\beta$ and serves as a proxy of dynamical heating. Finally, $M_{\rm tidal,b}$ is the bound mass contained within the tidal radius and serves as a proxy for the combined effects of the total tidal force and gravitational heating.

Some care needs to be taken in setting $\beta$ for the proxy $M_{\rm\beta,b}.$ The most straightforward option would be minimize the value of $f_{\rm removed}$ across all values of $\beta,$ but this procedure selects $\beta\approx 1.5,$ which will typically be within the halo’s own splashback shell. $M_{\rm\beta,b}$ therefore correlates with $c_{\rm vir}$ simply because the latter determines the mass distribution within the halo. Indeed, we find the Spearman rank coefficient $\rho_{S}(M(<1.5\times R_{\rm vir}),c_{\rm vir})=-0.26,$ even before any unbinding procedure has been used. Instead, we choose to set $\beta=3.$ At this distance, correlations between the total enclosed mass and $c_{\rm vir}$ are negligible, and $\beta\,R_{\rm vir}$ will generally be larger than $R_{\rm sp}.$ This choice has little effect on $f_{\rm removed}$, which remains approximately the same for $\beta\gtrsim 2$.

The proxies described above are strongly (anti-)correlated with local matter density. Thus, when we rank-order haloes using these proxies and make cuts, we need to distinguish this procedure from simple density cuts, which do not differentiate between particular physical processes that operate in high-density regions. To this end, we use the number of distinct haloes with 120 ${\rm\,km\,s^{-1}}<V_{\rm peak}<220{\rm\,km\,s^{-1}}$ located within $X$ comoving $h^{-1}\,$Mpc of the centre of a halo, $N_{X}$, as a proxy of the density of the local environment. We tested radii ranging from $1-10\,h^{-1}\,{\rm Mpc}$ and found that the assembly bias signal can be eliminated by removing the smallest fraction of haloes for $X=5$. We thus use $N_{5}$ as our fiducial local environmental density proxy.

In Fig. 3 we show the fraction, $f_{\rm removed},$ of all haloes rank-ordered by different proxies that must be removed to eliminate the assembly bias signal (see section 2.8). The corresponding $f_{\rm removed}$ thresholds for each proxy are presented in Table 1, and the red curves in Fig. 2 show clustering strength as a function of distance after such a cut has been made to $M_{\rm tidal,b}$. Note that statistical errors on the MCF are smaller relative to its amplitude than errors on the $\xi_{\rm high}(r)/\xi_{\rm low}(r)$ curve, which is one of the chief reasons that we use the former in calculations of $f_{\rm removed}$.

As discussed in section 2.8, the $1-\sigma$ uncertainties of the MCF shown in Fig. 3 are the sample variance uncertainties. Therefore, errors on the relative ordering of $f_{\rm removed}$ for multiple proxies within a single box are likely to be smaller than these error estimates. This can be also be seen by comparing the Bolshoi and BolshoiP points in Fig. 3. Although the difference between these two boxes is consistent with the estimated sample variance, the relative-ordering of the proxies by $f_{\rm removed}$ is quite similar between the boxes, with BolshoiP consistently having $f_{\rm removed}$ values one or two percentage points higher than Bolshoi for all proxies other than the high-scatter $M_{\rm\beta,b}.$ From this we can comfortably infer that the non-systematic error on $f_{\rm removed}$ is on the order of $1\%.$

The first feature apparent in Fig. 3 is that it is possible to remove assembly bias by making a cut on the local density, meaning that assembly bias is only present in high-density regions. This is consistent with models which predict that low-mass assembly bias is caused by non-linear processes, but is not necessarily a generic prediction of such models, as one could imagine assembly bias existing in all regions to different degrees of severity. Fig. 3 also shows that the portion of assembly bias which is not caused by misclassified splashback subhaloes is due to a small number of haloes in extreme environments: the cut $M_{\rm tidal,b}/M_{\rm vir}<1.36$ removes only 6% of all haloes but reduces assembly bias to statistically undetectable levels. For comparison, the cut to the density proxy $N_{5}$ removes assembly bias when 14% of haloes are removed.

Further testing shows that there are two reasons why assembly bias can be eliminated by removing only a small fraction of haloes. First, the mean value of $c_{\rm vir}$ ceases to be a strong function of these proxies once the haloes below the $f_{\rm removed}$ cutoff have been removed from the sample. Second, halo clustering strength varies strongly as a function of proxy value within the cutoff range, but is almost constant throughout the remaining sample.

Finally, Fig. 3 shows that $f_{\rm removed}$ for both $D_{\rm vir}$, a proxy for the single-halo tidal force, and $M_{\rm\beta,b}$, a proxy for dynamical heating, are at least as large as $f_{\rm removed}$ for $N_{5}$. Even if $M_{\rm\beta,b}$ uses values of $\beta$ small enough that it is primarily picking up features in the halo’s own density profile, $f_{\rm removed}$ stays above 0.25. This means that the effect of single-halo tidal forces and dynamical heating on assembly bias cannot be distinguished from the trivial effect of environmental density on halo bias, which means that neither can account for the assembly bias on their own.

In contrast, $f_{\rm removed}$ for $R_{\rm tidal}$, $M_{\rm tidal}$, and $M_{\rm tidal,b}$ are smaller than for $N_{5},$ indicating that these proxies are more closely connected to assembly bias than local density. The fact that $R_{\rm tidal}$ and $M_{\rm tidal},$ which are calculated using only the large-scale contribution to the tidal field, have lower $f_{\rm removed}$ than $D_{\rm vir}$ shows that it is the tidal force from large-scale structures, not from individual haloes, that play the dominant role in the assembly bias. $M_{\rm tidal,b}$ has the lowest $f_{\rm removed}$ and is thus the most closely connected to assembly bias of all the proxies we consider.

To summarize, the results of this and previous subsections show that $\approx 70\%$ of the low-mass assembly bias signal in $c_{\rm vir}$ is due to splashback subhaloes. The remaining $\approx 30\%$ of the signal is due to 10% of distinct haloes (6% of all haloes) that are affected by a combination of the truncation of their mass growth by large-scale tidal fields and dynamical heating caused by the collapse of sheets and filaments. There are thus three different physical processes that affect halo mass growth which all contribute significantly to the assembly bias signal.

3.3 The spatial and concentration distributions of the haloes responsible for assembly bias

In Fig. 4 we show the spatial distribution of splashback subhaloes located outside $R_{\rm vir}$ of their hosts in a 25 $h^{-1}$ Mpc thick slice of the Bolshoi simulation volume. In the same volume, we show the sets of distinct haloes that are removed under the criteria $N_{5}>18$ and $M_{\rm tidal,b}/M_{\rm vir}<1.36$, which each independently eliminate assembly bias. We also show the spatial distribution of a random $\times 5$ subsample of the set of distinct haloes that were not removed by the $M_{\rm tidal,b}/M_{\rm vir}<1.36$ cut.

Fig. 4 shows that both splashback subhaloes and haloes with low $M_{\rm tidal,b}/M_{\rm vir}$ values are strongly clustered in the fabric of the cosmic web: they lie within filaments, sheets, and nodes with characteristic scales of tens of $h^{-1}$ Mpc. Splashback subhaloes cluster strongly because they trace the spatial distribution of their massive host haloes, which are predominantly found in these dense environments. Haloes with low $M_{\rm tidal,b}/M_{\rm vir}$, on the other hand, are strongly clustered because the two physical processes that reduce $M_{\rm tidal,b}/M_{\rm vir}$ – strong tidal forces and gravitational heating – are strongest in similarly dense regions. The distributions of haloes with larger values of $M_{\rm tidal,b}/M_{\rm vir}$ or smaller values of $N_{5}$ are less clustered. We also provide a zoomed-in view of the distribution of these different groups in Figure 5.

Fig. 6 shows concentration distributions for different groups of haloes: all haloes outside $R_{\rm vir}$ of any host, splashback subhaloes outside $R_{\rm vir}$ that have not passed through pericentre of their orbit, splashback subhaloes outside $R_{\rm vir}$ which have passed their first pericentre, and distinct haloes outside the splashback shell of any host which have low $M_{\rm tidal,b}/M_{\rm vir}.$ This figure shows that the $c_{\rm vir}$ distribution of post-pericentre splashback subhaloes is biased to much larger values and are responsible for almost the entire high-$c_{\rm vir}$ tail of the overall concentration distribution. This indicates that the concentrations of such halos are affected substantially by the strong tidal interaction they experienced during their pericentre passage, which strips mass preferentially at the outskirts of haloes, thereby increasing their concentration (e.g., Kazantzidis et al., 2004). In contrast, splashback subhaloes that are on their first infall and distinct haloes with $M_{\rm tidal,b}/M_{\rm vir}<1.36$ have comparable concentration distributions and are only slightly shifted relative to the overall distribution of concentrations. The modest shift in $c_{\rm vir}$ is consistent with an older age of these haloes, rather than the large concentration boost in halos that have experienced tidal stripping. We note that the strength of the high $c_{\rm vir}$ tail becomes weaker if the halo sample is defined by $M_{\rm vir}$ or $V_{\rm max}$. This is because haloes that lost mass after their first pericentre passage drop out of the $M_{\rm vir}$-defined sample, but stay within the $V_{\rm peak}$ defined sample.

3.4 Time and mass dependence of assembly bias

As discussed in section 1, $c_{\rm vir}$ halo assembly bias has opposite signs at masses above and below the non-linear collapse mass scale, $M_{*}.$ Fig. 7 shows the dependence of assembly bias on $M_{\rm vir}/M_{*}$ in the Bolshoi simulation for the entire sample of distinct haloes (black line and shading), and samples in which subsets of haloes have been removed using different criteria discussed earlier in this section (coloured lines). We first divide haloes into logarithmic bins of $V_{\rm peak}$ with 0.08 dex width. We use bins of $V_{\rm peak}$ to be consistent with the rest of our analysis, although we show the assembly bias signal as a function of the median $M_{\rm vir}/M_{*}$ within each bin. For each bin above the convergence limit of $V_{\rm peak}=120\,\rm km\,s^{-1}$ we measure the MCF in the separation range of $4-8$ comoving $h^{-1}$ Mpc and split the simulation into eight equal-size sub-boxes to estimate the 1-$\sigma$ error on the MCF amplitude. To probe a wide range of $M_{\rm vir}/M_{\rm*}$ values, we use the $z=0$, $1$, $1.4$, and $3$ Bolshoi snapshots, with the $z=0$ snapshot giving us access to the lowest values of $M_{\rm vir}/M_{*}$ and $z=3$ giving us access to the highest.

The dependence of assembly bias on $M_{\rm vir}/M_{*}$ for distinct haloes outside $R_{\rm vir}$ of any larger host is consistent with the results of Wechsler et al. (2006). Removing splashback subhaloes reduces the assembly bias substantially at $M_{\rm vir}/M_{*}\lesssim 10$, but does not eliminate it completely, and does not affect assembly bias at larger masses. Removing further distinct haloes using $M_{\rm tidal,b}/M_{\rm vir}<1.36$ cut eliminates assembly bias entirely at $M_{\rm vir}/M_{*}\lesssim 5$, but likewise does not affect the assembly bias at larger $M_{\rm vir}/M_{*}$. This illustrates that the physical origin of assembly bias in the high-mass regime is not related to tidal forces or dynamical heating.

Interestingly, Fig.7 also shows that removing haloes using environmental density, $N_{5},$ does remove assembly bias at all $M_{\rm vir}/M_{*}$. Visual inspection reveals that this is because this cut removes the same spatial regions across time.

Given that halos and large-scale structure evolve with time, we also redid this analysis by removing a constant fraction of distinct haloes ranked by $M_{\rm tidal,b}/M_{\rm vir}$ and by $N_{5}$ at different redshifts rather than using a fixed cut as in Fig. 7. The results of such analysis are almost identical, albeit with slightly higher $f_{\rm removed}.$

3.5 Sensitivity to splashback subhaloes identification method

As discussed above, we use two different methods to identify splashback subhaloes: 1) haloes that move within $R_{\rm vir}$ of a larger halo at some point during their evolution and (“flyby subhaloes”) 2) haloes located within the splashback shell identified by the Shellfish algorithm (“splashback subhaloes”). It is clear that the samples of subhaloes identified using these methods cannot be identical: haloes that are on their first approach to a host and are already within the splashback shell but are still outside $R_{\rm vir}$ will be classified as splashback subhaloes by the second method, but not the first. Conversely, haloes that previously passed within $R_{\rm vir}$ of the host, but are now outside of the splashback shell identified by Shellfish will be identified by the first method, but not the second.

We find that $\approx 40\%$ of splashback subhaloes ($\approx 4\%$ of all haloes) are not identified as flyby haloes, but only $\approx 6-8\%$ of flyby subhaloes ($0.4-0.5\%$ of all haloes) are not identified as splashback subhaloes. The latter subhaloes are misidentified largely around host haloes below convergence limit of the Shellfish algorithm: if we restrict this analysis to host haloes that meet the convergence requirements of $N_{\rm 200m}>5\times 10^{4}$ and $\Gamma_{\rm DK14}>0.5$ (see Mansfield et al., 2017, for details), we find that only $1-2\%$ of flyby haloes ($\lesssim 0.1\%$ of all haloes) are not identified by the second method. This small fraction indicates that the splashback shells identified by Shellfish for well-resolved haloes capture the vast majority of the splashback subhaloes identified by the traditional subhalo trajectory method. This also indicates that the fraction of subhaloes ejected by three-body interactions via the slingshot process beyond the splashback shell (Kravtsov et al., 2004; Sales et al., 2007) is quite small and that most of the subhaloes outside $R_{\rm vir}$ are on their natural dynamical orbit around their host halo. We note that this conclusion should not be extended to $V_{\rm peak}<120$ without further testing: it is plausible that slingshot processes become more significant at lower masses.

This is consistent with earlier studies that analysed the radial distribution of flyby subhaloes (Ludlow et al., 2009; Wang et al., 2009; Li et al., 2013) and found that flyby subhaloes are common at distances up to $2\,R_{200c}$, with numbers decreasing quickly at larger radii but with a small population present out to $\approx 4\,R_{\rm 200c}$. We find that the radial distribution of flyby subhaloes is due to large size of the splashback shell relative to $R_{\rm 200c}$, its non-spherical shape, and the substantial scatter between $R_{\rm 200c}$ and the maximum radius of the splashback shell. For our sample of distinct haloes, the mean value of $R_{\rm sp}/R_{\rm 200c}$ is 2.16 and the mean value of $R_{\rm sp,max}/R_{200c}$ is 2.80 with a $1-\sigma$ scatter of $\approx 0.6$, where $R_{\rm sp,max}$ is the maximum radius of any point on the splashback shell.

Lastly, Figure 2 compares the MCF after both methods have been used to remove splashback subhaloes. The difference is small relative to the overall amplitude of the signal. We also find that when the procedure described in sections 2.8 and 2.8.1 is used, both classification schemes find similar cutoff values. However, $f_{\rm removed}$ is necessarily $\approx 3\%$ larger when flyby flagging is used to remove splashback subhaloes because these cuts must also remove infalling splashback subhaloes. The exception to this is the $M_{\rm tidal}$ proxy, which is higher for almost all splashback subhaloes than it is for almost all distinct haloes. Thus, $M_{\rm tidal}$ cannot remove assembly bias without removing the entire sample. This leads us to conclude that our general results are robust to differences in the subhalo classification scheme.

3.6 Comparison of Bolshoi and BolshoiP simulations

All analysis presented above was done using the Bolshoi simulation with cosmological parameters consistent with the final WMAP mission constraints (see section 2.1). To test the dependence of our results on the assumed cosmology, we repeated all analysis using the BolshoiP simulation, which assumes cosmological parameters consistent with the Planck mission constraints and found that all of the results are qualitatively consistent. The difference in $\Omega_{\rm m}$ in the Bolshoi and BolshoiP simulations leads to small changes in the cutoff values for $M_{\rm\beta,b}$, $D_{\rm vir},$ $M_{\rm tidal}$, $R_{\rm tidal},$ and $M_{\rm tidal,b}$, but the values of $f_{\rm removed}$ are within 0.01 of the values found for the Bolshoi simulation for all cuts, which the exception of the high-error $f_{\rm removed}$ value for our least efficient proxy, $M_{\rm\beta,b}.$

4 Discussion

4.1 Issues associated with proxy definitions

In this study we define and use several proxies of physical processes that could conceivably contribute to assembly bias. Of these, $D_{\rm vir},\,R_{\rm tidal},\,M_{\rm tidal},\,M_{\rm\beta,b},$ and $M_{\rm tidal,b}$ require estimating the local tidal force and/or determining whether a given particle is bound or unbound. However, it is not trivial to accurately determine whether a particle is bound in the outskirts of haloes (see, e.g., Behroozi et al., 2013a, for an extended discussion of related issues), and strong assumptions and approximations must be employed in the estimates of tidal forces. Errors made in estimating a particular proxy should result in additional scatter in its correlation with $c_{\rm vir}$ and should increase the uncertainty in our estimate of its contribution to the assembly bias. As a corollary, improvements in proxy definitions should only decrease the measured $f_{\rm removed}$ values. In practice, only $M_{\rm tidal,b}$ is strongly affected both by uncertainties in the tidal force estimate and by issues of identifying bound particles, which means that improvements in proxy estimates would primarily reduce $f_{\rm removed}$ for $M_{\rm tidal,b}$, while having an equal or lesser effect on our other proxies. This means that such improvements would not change our conclusions.

A detailed analysis of the errors associated with the approximations necessary for tidal force calculation can be found in Appendix C. A discussion of the issues related to identification of bound particles can be found in Appendix D.

4.2 Sensitivity of results to definitional choices

In this section we discuss the impact of the choices and assumptions made in our fiducial analysis on our results. We have already discussed how our choice of clustering statistic used to estimate assembly bias affects out results in sections 2.7 and 3.2 (see Fig. 2), so here we focus on the effect of our choices of $V_{\rm peak}$ for defining halo samples, $R_{\rm vir}$ as our reference halo boundary, and halo concentration as our formation time proxy. Although we did present justifications for our choices in section 2, it is important to assess how sensitive our results and conclusions are to these choices.

As an alternative to $V_{\rm peak}$, we could define halo samples using $M_{\rm peak},$ $V_{\rm max}$, or $M_{\rm vir}$. We find that sample selection by $M_{\rm peak}$ leads to results similar to our fiducial case, but using $V_{\rm max}$ or $M_{\rm vir}$ leads to a somewhat different behaviour. The amplitude of the MCF with only $R_{\rm vir}$ subhaloes removed is closer to the amplitude of the MCF with splashback subhaloes removed for a $V_{\rm peak}$ or $M_{\rm peak}$ cut. However, further removing splashback subhaloes with a $V_{\rm max}$ or $M_{\rm vir}$ cut results in only a small decrease in amplitude. This is because subhaloes generally experience significant mass loss and therefore sample selection based on their peak mass or circular velocity results in larger subhalo fractions compared to selection on current mass (cf., also Nagai & Kravtsov, 2005). The large-$c_{\rm vir}$ tail seen in Fig. 6 is weaker when $V_{\rm max}$ or $M_{\rm vir}$ is used to define the halo samples for the same reason. Our other results, such as the values of $f_{\rm removed}$ or the spatial distribution of different halo subsets, are largely unaffected. This is because the haloes removed by these cuts have merely had their accretion histories slowed: they have not experienced significant mass loss.

Although most analysis in this paper uses splashback shells as halo boundaries, we use spheres of radius $R_{\rm vir}$ as halo boundaries when we compute fiducial $\mathcal{M}(r)$ curves and when we classify “flyby” subhaloes. We have repeated these analyses using other commonly used values of $\Delta\bar{\rho}_{\rm m}$, and found that the main difference, unsurprisingly, is in the change of the amplitude of the reference MCF. Definitions with high values of density contrast, such as $\Delta\bar{\rho}_{\rm m}=200\rho_{\rm c}$ or $\Delta\bar{\rho}_{\rm m}=500\rho_{\rm c}$, result in a modest increase of the reference MCF amplitude, while definitions with low density contrasts, such as $\Delta\rho=200\rho_{\rm m}$, result in a modest decrease of the MCF amplitude. To decrease the reference MCF amplitude to the level of the MCF after removal of splashback subhaloes requires $\Delta\bar{\rho}_{\rm m}\approx 100\rho_{\rm m}$ – close to the typical density contrast enclosed by the splashback shell (see Fig. 12 in Mansfield et al., 2017). Changes in $\Delta$ used for radius definition have little effect on the amplitude of the MCF when flyby subhaloes have been removed because most subhaloes have their first pericentres at radii well within all of the commonly-used definitions of halo radius.

We have chosen to use $c_{\rm vir}$ as a proxy of halo age, but assembly bias behaves differently for different proxies (e.g., Villarreal et al., 2017; Mao et al., 2018; Salcedo et al., 2018), so one would reasonably wonder if halo removal criteria also depend on this choice. A full investigation of different definitions is beyond the scope of this work, but as a preliminary discussion, we repeat our analysis for the most commonly used alternative age proxy – the expansion factor at which the virial mass of the main progenitor of a halo was half of the halo’s current mass, $a_{1/2}$:

[TABLE]

Note that although large (small) $c_{\rm vir}$ values correspond to old (young) haloes, the opposite is true for $a_{1/2}.$ To simplify comparison with $c_{\rm vir}$-based results, we use $-a_{1/2}$ as the formation time proxy, so the sign of the MCF retains the same qualitative meaning.

Results for the $-a_{1/2}$ proxy are shown in Fig. 8, where we use the same cuts that removed the assembly bias in the $c_{\rm vir}$-based analysis. The figure shows that stricter cuts are required to eliminate $-a_{1/2}$ assembly bias. When we follow the procedure described in section 2.8, the $M_{\rm tidal,b}$ cut that removes $-a_{1/2}$ assembly bias results in $f_{\rm removed}=0.14,$ compared to $f_{\rm removed}=0.06$ for $c_{\rm vir}$ assembly bias. Other proxies experience similar increases in $f_{\rm removed},$ with the exception of $N_{5},$ which removes assembly bias for both definitions in almost all bins. The fact that Bolshoi and BolshoiP measurements of $f_{\rm removed}$ agree to within 0.01 for all variables (see §3.6) indicates that the differences in $f_{\rm removed}$ between the $c_{\rm vir}$ and $-a_{1/2}$ definitions are significant.

Another difference is that in contrast to the $c_{\rm vir}$ MCF, the amplitude of the $-a_{1/2}$ MCF does not reverse sign at large halo masses (see Fig. 8 and several previous studies: Gao et al. 2005; Gao & White 2007; Wetzel et al. 2007; Mao et al. 2018; Sato-Polito et al. 2018). This puzzling behaviour at first appears to be inconsistent with the physical origin of the high-mass assembly bias argued for by Dalal et al. (2008). However, Chue et al. (2018) showed that $a_{1/2}$ and similar measures of formation time can be problematic if measured relative to a standard overdensity mass. Because these definitions do not account for mass in the splashback shell, haloes measured at a constant $M_{\rm vir}$ actually have a range of “true” splashback-enclosed masses, and the high-mass haloes will preferentially have early $a_{1/2},$ which increases the level of bias measured for early $a_{1/2}$ haloes. The intersection of this behaviour with Fig. 8 requires further study.

4.3 Comparison with previous work

The effect of splashback subhaloes on assembly bias was investigated in a number of recent studies (Wang et al., 2009; Li et al., 2013; Sunayama et al., 2016), which concluded that splashback subhaloes contribute significantly to low-mass halo assembly bias, but cannot account for the entire signal. These studies used the “flyby” approach to classify subhalos (similar to the method described in 2.3.1), which can include bona fide splashback subhaloes, as well as subhaloes ejected via slingshot effect after dynamical interactions with other subhalos. However, this classification cannot account for a large number of subhaloes within splashback shells that are on their first infall. This left open the question of the contribution of such infalling subhaloes on assembly bias. In this study we answer this question in section 3.5.

Additionally, these studies have not demonstrated the physical origin of the remaining assembly bias signal. Our results differ from those of Sunayama et al. (2016), which find that splashback subhaloes have little effect on the MCF at large distances ($R\,\gtrsim\,10\,h^{-1}\,$Mpc). Sunayama et al. (2016) used the same simulation and underlying halo catalogues as this paper, so this difference is likely due to the fact that their samples are defined by $M_{\rm vir}$ (see section 4.2) and their use of halo bias ratios to measure assembly bias. This statistic results in larger errors than the MCF, as we discussed in § 2.7 and 3.2. For example, comparison of the bin-to-bin scatter in Fig. 3 and Fig. 4 of Sunayama et al. (2016) to the 1-$\sigma$ error contours in the right panel of our Fig. 2 indicates that their measurements may not have been sensitive enough to probe large-scale assembly bias.

The conjecture that non-linear tidal and dynamical heating effects can be responsible for low-mass halo assembly bias was discussed in a number of studies (e.g. Wang et al., 2007; Dalal et al., 2008; Hahn et al., 2009; Wang et al., 2011; Hearin et al., 2016b; Paranjape et al., 2018; Salcedo et al., 2018; Musso et al., 2018). Often, arguments for a particular process are based on establishing existence of a correlation between halo formation time or halo bias and a proxy for a particular process, such as $R_{\rm hill}$, bound mass fraction, the magnitude of tidal eigenvectors, or various measures of tidal anisotropy. While such correlations provide useful information, by themselves they are not sufficient to establish that a given physical process is responsible for assembly bias. This is because the proxies of all these processes are all strongly correlated: an explicit comparison, such as that shown in Fig. 3, is more direct and compelling in identifying the responsible process.

Furthermore, looking at the global connection between halo formation time and a proxy is problematic for two reasons. First, we show that after splashback subhaloes are removed low-mass halo assembly bias is due to only a small fraction of distinct haloes. Thus, analysis relying on the global correlation strength is not optimal. Second, as discussed in section 2.8.1, the correlation between a proxy and halo formation time by itself contains no information about how closely that proxy is related to assembly bias: a strongly correlated proxy which experiences weak differential clustering, such as $M_{\rm\beta,b}$, will not contribute to assembly bias. We avoid both these issues with the procedure described in section 2.8.

We find a strong connection between tidal forces from the large-scale mass distribution and assembly bias in agreement with the conclusions of Hahn et al. (2009), Hearin et al. (2016b), and Salcedo et al. (2018). However, in contrast with these studies, we find that this this effect cannot be effectively approximated by assuming that haloes only feel the tidal force of their most gravitationally dominant neighbor. In fact, we find that when such an approximation is made, the connection is sufficiently weak that it is likely caused simply because the Hill radius is a crude estimate of local density (see section 3.2). This discrepancy is due to two factors. First, some of these studies do not perform the type of multi-variate analysis that would be necessary to differentiate between different contributing physical processes. Second, while these studies effectively map out the the connection between formation time and single-halo tidal proxies, this is unrelated to the connection between assembly bias and these proxies, an argument we make in section 2.8.1.

Also, contrary to the conclusion of Paranjape et al. (2018), we do not find a compelling evidence that large-scale tidal anisotropy contributes significantly to assembly bias beyond what is expected from its correlation with tidal field strength. Nevertheless, given the inaccuracies associated with all methods based on second-order approximations to the tidal field (see Appendix C), a more detailed study of tidal field anisotropy could prove fruitful, especially in the context of explicitly studying the tidal environments within structures like sheets and filaments.

Although our results are in qualitative agreement with the conjectures of Wang et al. (2007) and Dalal et al. (2008) that gravitational heating is a significant component of assembly bias, we also find that this connection only becomes strong when tracers simultaneously incorporate both gravitational heating and a halo’s zone of influence over the local tidal field, such as $M_{\rm tidal,b}.$

Our work uses an approach similar to that of Villarreal et al. (2017), so we have performed an in-depth comparison with their results. We find broad qualitative agreement between our $D_{\rm vir}$ results and the results of Villarreal et al. (2017), but find that quantitatively the $\Delta$ values they report imply $D_{\rm vir}$ values smaller than our findings by $\approx 25\%$.

This difference is due to two factors. First, we find that the sample variance in the boxes used by Villarreal et al. (2017) is larger than than they estimated. We estimate the sample variance using subvolumes of the BolshoiP simulation, which has nearly identical mass resolution and cosmology to the CPla_L0125 box used by Villarreal et al. (2017) and find that the actual variance is larger than the uncertainty they estimated by repeatedly shuffling marks among haloes. This means that MCFs in the CPla_L0125 box could be lower due to sample variance, which could thus result in less aggressive conditions for the removal of assembly bias. Second, while the Rockstar halo finder used in Villarreal et al. (2017) and in this paper is a state-of-the-art tool for measuring the properties of haloes with density contrasts of $\Delta\gtrsim\Delta_{\rm vir}$ (e.g., Knebe et al., 2013), it cannot effectively measure halo properties at lower density contrasts, such as the $\Delta=20$ contrast used by Villarreal et al. (2017). This is because there is no FOF linking length which can fully percolate all matter out to such large overdensity radii, while also allowing for efficient load-balancing. This leads to underestimates of halo masses and artefacts in the density profile. Contrary to the findings of Villarreal et al. (2017), we find that even with an unusually large linking length of $b=0.5,$ virtually all haloes have underestimated $M_{\rm 20m}$ masses. The ratio $M_{\rm true}(<R_{\rm 20m,Rockstar})/M_{\rm 20m,all,Rockstar}$ has $1-\sigma$ contours of $1.04-1.13$, with $2\sigma$ fluctuations reaching $\approx 2.5.$ The magnitude of underestimation is significantly worse at more commonly-used linking lengths. This and the fact that haloes no longer follow NFW profiles at large radii (e.g., Becker & Kravtsov, 2011; Diemer & Kravtsov, 2014), adds biases and noise to the $R_{\rm s}$ and $R_{\rm s{\rm Klypin}}$ values measured by Rockstar. This, in turn, artificially reduces the amplitude of the MCF. We discuss this issue in greater depth in Appendix A.

Our interpretation is consistent with the test presented in Figure 11 of Villarreal et al. (2017), which shows that when $R_{\rm 20m}$ is used to exclude subhaloes, but concentrations are measured from the halo catalogues constructed using larger $\Delta$, the MCF is not consistent with zero. We find that when we replicate their analysis using manually-constructed overdensity profiles, excluding “subhaloes” by spheres of radius $R_{\rm 20m}$ is no longer capable of mitigating assembly bias. Larger overdensity radii that are comparable with our reported $D_{\rm vir}$ cutoff values are required.

4.4 Directions for future work

In this paper, we focus on the dependence of halo bias on $c_{\rm vir}$, but galaxy properties are likely related to a number of halo properties. This means that the effects of secondary biases on galaxy clustering may not be confined to the $c_{\rm vir}$ bias dependence and may remain even if $c_{\rm vir}$ dependence of bias is removed. As discussed in 4.2, this is true even for an alternative choice of halo formation time proxy, $-a_{1/2}$. Furthermore, Villarreal et al. (2017) show that mitigating secondary biases in axis ratio and spin parameter is more difficult than removing bias in $c_{\rm vir}$, although the discussion in section 4.3 should be kept in mind when assessing these results. It would be useful to perform analysis comparable to the one presented here for a number of other key halo properties to build a more complete understanding of the physical origin of the corresponding dependencies of halo bias.

One of our key results is that despite the large contribution of tidal forces to $c_{\rm vir}$ assembly bias, this cannot be shown conclusively when using rough and inaccurate estimates of the tidal force, such as $D_{\rm vir}$ or the single-halo $R_{\rm hill}$. Although $R_{\rm tidal}$ defined in section 2.5.2 accounts for tidal forces from multiple haloes, it still is rather inaccurate, as we show in Appendix C and discuss in 4.1. The accuracy of the tidal force estimate can be improved by using a higher order approximation of the tidal field, by iteratively recalculating $R_{\rm tidal}$ and removing nearby sources accordingly, or by explicitly evaluating the tidal field outside the halo and identifying the turnover associated with the tidal radius directly. More accurate estimate of the tidal radius could result in a better identification of haloes responsible for assembly bias. This effort could also be aided by incorporating cosmic web classifiers (see review by Libeskind et al., 2018), which would allow higher accuracy analytic calculation of the tidal fields associated with nearby large-scale structure, rather than the low-order approximations that are required for generic point distributions.

Finally, all of the tidal field and dynamical heating proxy estimators discussed and used in this paper are computed from simulated quantities and cannot be immediately be applied to observations. A follow-up exploration of possible observable proxies that can remove particular flavors of secondary halo bias using mock catalogues will be a useful future avenue of research.

5 Summary and Conclusions

In this study, we present a detailed analysis of the physical causes of assembly bias - the dependence of halo clustering on proxies of halo age. We focus on the origin of halo assembly bias for haloes within the mass range that hosts typical galaxies and on ages traced by $c_{\rm vir}$, but we also present some of our key results across a broad range of masses and for other definitions of age (see Fig. 7 and 8).

We first explore the contribution of “splashback subhaloes” to assembly bias, where splashback subhaloes are defined either as halos that have passed within the virial radius of a larger halo at some point in the past or as haloes are located within the splashback shell of a larger halo, as determined by the method of Mansfield et al. (2017). Assembly bias is measured both before and after the removal of these subhaloes. We show that splashback subhaloes are responsible for about two thirds of the assembly bias signal, but do not account for the entire effect. Moreover, it is the subhaloes that have passed the pericentre of their orbit at least once that are responsible for the contribution of subhaloes to assembly bias. In addition, we find that the high-$c_{\rm vir}$ tail of the distinct halo distribution is due almost entirely to these same post-pericentre subhaloes.

At the mass ranges considered in this paper, we find that the fraction of haloes which have passed within the splashback shells of their hosts but are later located outside them is small, which indicates that the fraction of haloes ejected beyond the splashback shell due to three-body interactions is small.

We then investigate which additional physical processes contribute to assembly bias. We do this by constructing proxies of these processes for each halo, ranking distinct haloes according to each proxy, and measuring what fraction of the ranked halos need to be removed in order for the assembly bias signal to be statistically consistent with zero. We find that assembly bias is caused by a relatively small number of haloes in dense regions. These haloes have had their accretion histories truncated by a combination of large-scale tidal fields and the high velocities of ambient particles. We also demonstrate that neither process can cause assembly bias on its own and that these tidal fields are not well-modeled by assuming that the dominant tidal contribution comes from a single massive neighbor. We further argue that the commonly-used approach of measuring the correlation between a physical proxy and halo age cannot be used on its own to draw conclusions about the strength of the connection between that proxy and assembly bias.

A key finding of this study is that after splashback subhaloes are removed, the residual assembly bias is due to only 10% of distinct haloes (5% of all haloes). To summarize, 27% of haloes are removed due to a traditional $R_{\rm vir}$-based subhalo cut, a further 10% are removed due to a splashback subhalo cut, and finally 6% of all haloes are removed due to the cut based on $M_{\rm tidal,b}/M_{\rm vir}.$ These low $M_{\rm tidal,b}/M_{\rm vir}$ haloes are located within the largest filaments and are only slightly more concentrated than the general population. However, their strong spatial clustering results in an outsized effect on the global assembly bias signal.

We find that in the WMAP cosmology, the removal of haloes above a certain local density, as measured by the number of haloes within $5\,h^{-1}$ comoving Mpc, $N_{5}>18$, removes assembly bias for both $c_{\rm vir}$ and $a_{1/2}$ across all mass scales and redshifts and that a similar cut exists in a Planck cosmology. Such a cut removes a much larger fraction of haloes from the sample than the cut in $M_{\rm tidal,b}/M_{\rm vir}$, and thus does not correspond to a real physical process contributing to assembly bias. Nevertheless, this result indicates that it may be fruitful to explore whether density-based cuts on mock galaxy catalogs can be used to remove assembly bias from galaxy samples and motivates further studies in this direction.

Acknowledgements

We would like to thank Antonio Villarreal and Andrew Zentner for useful discussions during this work and Neal Dalal, Andrew Hearin, and Surhud More for comments that improved the paper. We would also like to thank organizers and participants of the "Quantifying and Understanding the Galaxy–Halo Connection" conference held at the Kavli Institute for Theoretical Physics in May 2017 for stimulating discussions and interactions that motivated this work. PM would like to thank Alice Mansfield for moral support while writing this paper. AVK would like to thank Institute of Astronomy at Cambridge University and its Sackler visitor program for the warm hospitality during the completion of this paper. The work presented in this paper was supported by the NSF grant AST-1714658 and by the Kavli Institute for Cosmological Physics at the University of Chicago through grant PHY-1125897 and an endowment from the Kavli Foundation and its founder, Fred Kavli. Analyses presented in this paper have been carried out using the Midway cluster at the University of Chicago Research Computing Center, which we acknowledge for support. The CosmoSim database used in this paper is a service by the Leibniz-Institute for Astrophysics Potsdam (AIP). The Bolshoi simulations have been performed within the Bolshoi project of the University of California High-Performance AstroComputing Center (UC-HiPACC) and were run at the NASA Ames Research Center. Analyses presented in this paper were greatly aided by the following free software packages: NumPy (Van Der Walt et al., 2011), SciPy (Jones et al., 2001-2016), Matplotlib (Hunter, 2007), and GitHub. We have also used the Astrophysics Data Service (ADS) and arXiv preprint repository extensively during this project and the writing of the paper.

Appendix A Effects of halo definition on concentration in the Rockstar halo finder

The Rockstar halo finder works by dividing the simulation into 3D friends-of-friends (FOF) groups, adaptively creating smaller 6D FOF groups in phase space, placing halo centres at the most refined 6D FOF groups, and finally calculating halo properties relative to those centres (Behroozi et al., 2013c). The size of the initial 3D FOF groups is set by the input linking length in units of the mean interparticle separation, $b$. The accuracy of the halo properties computed by Rockstar depends on the original 3D FOF groups percolating out to the baseline overdensity radius of the corresponding haloes (More et al., 2011; Behroozi et al., 2013c). This is particularly important when fitting halo density profiles: if halo boundaries extend into the unpercolated regions of the FOF group, the density in the outermost radial bins will be systematically underestimated, shifting the location of profile features. In this Appendix, we examine the effect of using different halo definitions on the measured concentrations.

Behroozi et al. (2013c) perform convergence tests that show that using linking length of $b=0.28$ leads to full percolation within $R_{\rm vir}$. They note that when one defines halo boundary larger than $R_{\rm vir}$, a larger linking length should be used and additional tests should be performed to ensure full percolation within such boundary. We test the effect of the halo boundary choice by running Rockstar repeatedly on the CBol_L0125 simulation from Diemer & Kravtsov (2015) for a variety of overdensity radii, $R_{\Delta}$, with overdensity ranging from  $\Delta=20\rho_{\rm m}$ to $\Delta=1600\,\rho_{\rm m}$ for $b=0.28,0.5$ and for a reference catalogue with $\Delta=\Delta_{\rm vir}$ and $b=0.28$. We then matched haloes across the catalogues to our reference $R_{\rm vir}$ catalogue. Our tests indicate that results are not sensitive to the way this matching is done, so we use a simple procedure where a halo is considered a ‘‘match’’ if its centre lies within 0.25 kpc of the centre of a counterpart in the $R_{\rm vir}$ catalogue. This criterion is sufficient to unambiguously match most haloes, but in the event that multiple haloes meet it, we match to the halo in that group with the closest $M_{\rm vir}$ to the reference halo. Subhalo and distinct halo status are not factored in to this matching. We restrict our sample to haloes classified as hosts by the $\Delta_{\rm vir}$ catalogue with $10^{11.5}\,h^{-1}\,M_{\odot}<M_{\rm vir}<10^{12}\,h^{-1}\,M_{\odot},$ as measured by the same catalogue. The choice of mass range has only a slight effect on results.

Fig. 9 shows the ratio of $R_{\rm s}$ measured in catalogues constructed for different values of $\Delta$, denoted $R_{\rm s,\,\Delta},$ to $R_{\rm s}$ measured in the reference catalogue, denoted $R_{\rm s,\,vir}$. We show this ratio as a function of $\Delta$ for two values of $b$. We find that creating catalogues with larger linking lengths takes an inordinate amount of time, presumably because a large fraction of the simulation is placed into the same FOF group. The figure shows that .$R_{\rm s}$ measurements for $\Delta\lesssim 200$ are biased low relative to the values found for $\Delta=\Delta_{\rm vir}$ and there is a significant scatter between the two. The bias is about twice larger for $b=0.28$ compared to $b=0.5$.

The primary implication of this result is that concentrations measured for haloes identified and analyzed by Rockstar using low $\Delta$, such as $\Delta=20\,\rho_{\rm m},$ should not be trusted due to large systematic bias and scatter. This is due to lack of FOF percolation in the outskirts of haloes, which biases densities in the outskirts low and this, in turn, biases the best-fit $R_{\rm s}$ values low. This is also true for other ways of estimating concentration, such as deriving it from $V_{\rm max}/V_{\rm vir}.$ The context of this fact in relation to our work on assembly bias is discussed 4.3.

Behaviour at commonly used choices, such as $\Delta=200\,\rho_{\rm m},\,200\,\rho_{\rm c},\,500\,\rho_{\rm c}$ is also noteworthy. For $b=0.28,$ the systematic biases on $R_{\rm s}$ for these three definitions relative to our reference $\Delta=\Delta_{\rm vir}$ catalogue are +6%, -5%, and -5%, respectively. While the difference for $\Delta=200\rho_{\rm m}$ has contributions from lack of percolation, the difference between $\Delta_{\rm vir}$ and higher density definitions must be due to a different effect, such as deviations of halo profiles from the fitted NFW form. Any attempt to compare, for example, mass-concentration relations to the $\approx 5\%$ level measured with different primary definitions should account for this effect.

Lastly, as discussed in section 2.3, overdensity radii are fundamentally unphysical choices for halo boundaries, and $\Delta_{\rm vir}$ cannot be thought of as a more ‘‘correct’’ choice than other nearby overdensities. Consequently, Fig. 9 should not be interpreted as showing deviations from the true value of $R_{\rm s}$, but merely deviations from a particular reference value where the FOF groups are known to be percolated.

Appendix B Fast Halo Containment Checks

Numerous components of the analysis presented in this paper rely on containment checks, particularly when computing subhalo status, computing $R_{\rm tidal}$, or computing $M_{\rm tidal}$. Out sample contains $\approx 300,000$ haloes and the Bolshoi simulations contain $2048^{3}$ particles each, so a naive $N^{2}$ check of every pair of objects would be prohibitively expensive. This is particularly true when identifying splashback subhaloes through the surfaces found by Shellfish because Shellfish represents splashback shells using third-order Penna-Dines surfaces, which take roughly fifty math library function calls to evaluate. In this Appendix, we describe our approach for computing containment checks.

First, suppose we are given a set of points $P$, which we must check for containment within a set of haloes, $H.$ First, we construct a uniform 3D grid spanning the simulation volume and place all elements of $P$ within lists associated with each cell in the grid. Then, for each halo in $H$ we construct a bounding box fully enclosing its boundary and compute containment checks for only the particles which reside in grid cells that intersect with it. Because the lists associated with each grid cell are created once and potentially iterated over many times, we represent lists as dynamically allocated arrays instead of as linked lists to increase cache locality. We find that for a grid with $250^{3}$ cells containment checks are no longer a significant component of the runtime cost of any analysis in this paper.

In the case where a halo boundary is determined by an expensive function $f(\phi,\theta)$, such as the Penna-Dines functions used by the Shellfish code to approximate splashback shells, we use the following procedure to accelerate containment checks. First, for every halo in $H,$ we compute the minimum and maximum values of $f(\phi,\theta)$, $f_{\rm max}$ and $f_{\rm min}$. Since a point at a distance $r$ is automatically contained if $r<f_{\rm min}$ and automatically not contained if $r>f_{\rm max}$, we only evaluate $f(\phi,\theta)$ if the points is at a distance, $r$, that satisfies $f_{\rm min}<r<f_{\rm max}$.

Appendix C Tidal Force Errors

In this Appendix, we investigate some of the error properties of the tidal radius and discuss an important approximation made in our calculation of $R_{\rm tidal}$, the inclusion of a minimum cutoff radius when adding contributions to a halo’s local tidal tensor. For the purpose of clarity, we will refer to the tidal radius calculated after a second-order approximation of the external gravitational potential has been made as $R_{\rm tidal}^{(2)}$ and the tidal radius when the exact tidal field is used as $R_{\rm tidal}^{(\rm ex)}$. Other quantities will use an analogous referencing scheme. Elsewhere in this paper, $R_{\rm tidal}^{(2)}$ is referred to as $R_{\rm tidal}$. We also take the convention that the Lagrangian point between a halo and an external source is $L_{1}$ and that the point on the opposite side of the halo is $L_{2}.$ In cases where analysis is performed on haloes without a single external source, $L_{1}$ is the Lagrangian point with the lowest external potential.

Like the classical $R_{\rm Hill}$ derivation, our calculation of $R_{\rm tidal}^{(2)}$ (see Equation 8) assumes that the external tidal field felt by the halo is well-approximated by a second-order hyper-paraboloid. This is necessary because the tidal tensor which is used to determine the principle components of the tidal field only contains second derivatives of the gravitational potential. Note that in the special case where there is only a single external point source,

[TABLE]

so the discussion below can be extended to error analysis on the classical Hill radius. Note that the factor of $(3/2)^{1/3}$ is because the derivation of $R^{(2)}_{\rm Hill}$ assumes that the halo is on a circular orbit around the external point source and thus experiences a centrifugal force in addition to a tidal force, while the derivation of $R_{\rm tidal}^{(2)}$ assumes that all non-tidal pseudo-forces are zero. Given the scale of errors discussed below, and the fact that this factor decreases as haloes deviate from circular orbits -- a configuration which is very rare for distinct haloes -- we do not consider this difference to be significant.

We perform two complementary tests on the accuracy of $R_{\rm tidal}^{(2)}.$ First, we analytically compute $R_{\rm tidal}^{(2)}/R_{\rm tidal}^{(\rm ex)}$ for a single source at a distance $r$ from a halo, and second, we measure exact the error on $F_{\rm tidal}^{(2)}$ for our halo sample and combine this with reasonable assumptions about the shape of the tidal field to estimate upper limits on $R_{\rm tidal}^{(2)}/R_{\rm tidal}^{(\rm ex)}.$

We show the results of this first calculation in the left panel of Fig. 10. Although there is no closed-form expression for $R_{\rm tidal}^{(\rm ex)},$ it can be found numerically by maximizing the effective potential. We parameterize the error as a function of $r/R_{\rm tidal}^{(2)}$ which also absorbs the dependence on the mass ratio. We recover the well-known fact that as the mass ratio between the halo and the external source decreases and the tidal radius increases, the two Lagrangian points become asymmetric and that errors become increasingly significant. This can also be interpreted as an estimate of the error associated with a particular value inner cutoff radius for $r=R_{\rm min}$ when following the procedure described in section 2.5.2. This can be considered a worst-case estimate of the error at a given $R_{\rm tidal}^{(2)}$ because the true matter distribution will generally contain many points at distances larger than $R_{\rm min}$ which contribute significantly to the tidal field.

We perform our second test by first computing $L_{1}^{(2)},$ $L_{2}^{(2)},$ $R_{\rm tidal}^{(2)},$ and $F_{\rm tidal}^{(2)}$ for every halo in our sample for a particular choice of $R_{\rm min}.$ Then, we use the raw particle data to compute the radial and tangential components of $F_{\rm tidal}^{(\rm ex)}(L_{1}^{(2)})$ and $F_{\rm tidal}^{(\rm ex)}(L_{2}^{(2)})$ for these haloes. Particles within $R_{\rm min}$ are not included in this calculation. To obtain an estimate of $R_{\rm tidal}^{(\rm ex)}$ from this, we make two simplifying assumptions about the the shape of the tidal field. First, we assume that the exact Lagrangian points lie along the same axis as a halo’s second-order Lagrangian points. We find that the tangential components of $F_{\rm tidal}^{(\rm ex)}(L_{1}^{(2)})$ and $F_{\rm tidal}^{(\rm ex)}(L_{2}^{(2)})$ are small compared to the radial components, implying that this is a reasonable assumption. Second, we assume that along the lines connecting $L_{1}^{(2)}$ to $L_{1}^{(\rm ex)}$ and $L_{2}^{(2)}$ to $L_{2}^{(\rm ex)},$ the tidal force varies slowly enough that it can be well approximated by

[TABLE]

Here, $\alpha$ is an arbitrary constant which varies from halo to halo and may be different for different Lagrangian points within the same halo. It represents the deviation from the scaling seen when the tidal potential is approximated to second order. In this case,

[TABLE]

where $i$ indexes over Lagrangian points. This assumption is informed by tests on single-source effective potentials, which find that for all but the smallest external point sources, $-1<\alpha<0$. For haloes where $R_{\rm min}>R_{\rm tidal}^{(2)},$ we would expect that at a constant $R_{\rm tidal}^{(2)}$ the tidal field would be varying more quickly at $L_{i}^{(2)}$ when the field is generated by a single point source than when it is generated by a more diffuse matter distribution, so it’s likely that this range of $\alpha$ values holds for our simulated haloes as well. For this reason we can place the following upper bound on the error in tidal radius:

[TABLE]

We show the fractional error in $R_{\rm tidal}$ using this limit in the right panel of Fig. 10 as a function of the adopted $R_{\rm min}.$ Errors balloon uncontrollably for $R_{\rm min}\lesssim 4\,R_{\rm vir}$ but are more well-behaved at larger radii, with errors dropping to the $\approx 10\%$ level at $\approx 10\,R_{\rm vir}.$

One interesting feature of this Figure is that for $R_{\rm min}\lesssim 2\,R_{\rm vir},$ the error on the location $L_{1}$ becomes positive. This is likely because this is the characteristic size of the splashback radius, meaning that the halo’s own particles will be incorporated into the calculation of $F_{\rm tidal}^{(\rm ex)}.$ Since the tidal force is repulsive, this inclusion of halo particles will reduce the apparent strength of the field and increase $R_{\rm tidal}.$

These tests indicate that for $R_{\rm min}=10\,R_{\rm vir},$ the errors in $R_{\rm tidal}^{(2)}$ which are specifically due to the second order approximation of the tidal field are small. However, this analysis is performed at a constant $R_{\rm min},$ so it doesn’t account for errors due to the removal of significant sources close to the halo. This is not an issue for our analysis because $R_{\rm tidal}$ is explicitly a proxy for the large-scale tidal field, and our proxy $D_{\rm vir}$ is better suited for close sources. This would, however, become a significant issue for studies which need $R_{\rm tidal}$ for purposes other than rank-ordering haloes. Further discussion on the impact that improvements in the accuracy of $R_{\rm tidal}$ would have on our results can be found in section 4.1.

More generally, while the issue of measuring tidal radii around haloes with only a single significant source is well-explored (see §2 in van den Bosch et al., 2018, for a review), and the tidal radius due to the large scale field can be measured effectively with the tidal tensor, there currently does not exist an effective method for combining these two regimes. We outline a number of potential approaches which could be used to address this issue in 4.4, but consider the testing and calibration of such methods to be beyond the scope of this paper.

Appendix D Identifying Bound Particles in Halo Outskirts

While the concept of gravitational binding is straightforward to define for particles near the center of a non-accelerating halo, the same is not true for particles in the outskirts of haloes, especially those experiencing a strong tidal force. These difficulties arise from two key areas: first, it is difficult to disentangle the potential caused by a halo from the potential of its surroundings. Although haloes have a non-trivial amount of mass stored outside $R_{\rm vir},$ the so-called ‘‘two-halo’’ term starts to dominate the density distribution at $r\gtrsim 1-2R_{\rm vir}$ (e.g. Diemer & Kravtsov, 2014), meaning that any calculation of the potential which is done directly from the density profile or from the particle distribution must be done with care. Second, for particles near the tidal radius, the effective potential due to the external tidal field becomes significant. While this issue could in principle be solved by defining escape velocities relative to the minimum potential at either Lagrangian point, it also means that boundedness calculations will suffer from the same accuracy issues as the tidal radius calculations (see Appendix C).

The effect of tidal forces on particle escape velocities presents another issue for the analysis in this paper, specifically. If this effect is taken into account, it means that gravitational heating and tidal forces can no longer be disentangled. Even a ‘‘control’’ variable like $M_{\rm\beta,b}$ would depend on the tidal field, and could potentially make gravitational heating appear to be more a more significant contributor to assembly bias than it actually is. Primarily because of this reason, and to a lesser extent because of the issues described in the previous paragraph, we take on the simple and standard boundedness condition given in Eq. 10, but note that the $f_{\rm removed}$ value for $M_{\rm tidal,b}$ could become even lower if more sophisticated approaches were used.