Elucidating the impact of massive neutrinos on halo assembly bias

Massive neutrinos have non-negligible impact on the formation of large-scale structures. We investigate the impact of massive neutrinos on the halo assembly bias effect, measured by the relative halo bias $\hat{b}$ as a function of the curvature of the initial density peak $\hat{s}$, neutrino excess $\epsilon_\nu$, or halo concentration $\hat{c}$, using a large suite of $\Sigma M_\nu{=}0.0$ eV and $0.4$ eV simulations with the same initial conditions. By tracing dark matter haloes back to their initial density peaks, we construct a catalogue of halo twins that collapsed from the same peaks but evolved separately with and without massive neutrinos, thereby isolating any effect of neutrinos on halo formation. We detect a $2\%$ weakening of the halo assembly bias as measured by $\hat{b}(\epsilon_\nu)$ in the presence of massive neutrinos. As there exists a significant correlation between $\hat{s}$ and $\epsilon_\nu$ ($r_{cc}{=}0.319$), the impact of neutrinos persists at a reduced level~($0.1\%$) in the halo assembly bias measured by $\hat{b}(\hat{s})$. However, we do not detect any neutrino-induced impact on $\hat{b}(\hat{c})$, consistent with earlier studies and the lack of correlation between $\hat{c}$ and $\epsilon_\nu$ ($r_{cc}{=}0.087$). We also discover an analogous assembly bias effect for the neutrino haloes, whose concentrations are anti-correlated with the large-scale clustering of neutrinos.


INTRODUCTION
Massive neutrinos play a key role in the formation of large-scale structures (LSS) of the Universe (Lesgourgues & Pastor 2006;Wong 2011).Their lack of clustering on scales below the free-streaming length induces a scale-dependent growth of the LSS, which enables stringent constraints of the sum of neutrino masses Σ  with cosmological probes (Cuesta et al. 2016;Vagnozzi et al. 2017;Brinckmann et al. 2019;Dvorkin et al. 2019;Giusarma et al. 2016;Planck Collaboration et al. 2020;Alam et al. 2017;Ivanov et al. 2020;Palanque-Delabrouille et al. 2020;Di Valentino et al. 2021).Another important effect caused by the massive neutrinos is the differential growth of dark matter haloes between environments of different neutrino-todark matter density ratios (Yu et al. 2017).Through this effect, the presence of massive neutrinos may cause an impact on the so-called "halo assembly bias" (hereafter referred to as HAB; Gao & White 2007), potentially opening a new avenue for measuring Σ  with LSS observations.However, such impact can only be weak at best, as it has so far evaded detections in the simulations (Lazeyras et al. 2021).In this paper, we focus on the "twin" haloes collapsed from the same initial density peaks but evolved separately with and without massive neutrinos, aiming to elucidate the impact of massive neutrinos on the assembly bias of massive haloes.
HAB has been robustly predicted by the hierarchical structure formation of Λ+cold dark matter (CDM) simulations (Sheth & Tormen 2004;Gao et al. 2005;Wechsler et al. 2006;Harker et al. 2006; Jing ★ E-mail: yingzu@sjtu.edu.cnet al. 2007;Li et al. 2008).It has increasingly been referred to as the "secondary bias" (Salcedo et al. 2018;Mao et al. 2018;Chue et al. 2018;Montero-Dorta et al. 2020;Contreras et al. 2021;Lin et al. 2022;Lazeyras et al. 2023;Wang et al. 2023) -the large-scale bias of haloes depends not only on halo mass, but also on other halo properties such as concentration, formation time, and spin.At any given halo mass  ℎ , the strength of HAB is often characterised by the relative halo bias b as a function of some "indicator" property , where (| ℎ ) is the bias of haloes with  and ⟨| ℎ ⟩ is the overall bias at  ℎ .Contreras et al. (2021) found that HAB is independent of cosmological parameters in the ΛCDM, at least for the high-mass haloes that we focus on in this paper.Such an independence on cosmology could be useful for measuring Σ  , if HAB turns out to be sensitive to massive neutrinos.
Observationally, HAB at the high mass end has been at the brink of detection recently (Yang et al. 2006;Lin et al. 2016;Miyatake et al. 2016;Zu et al. 2017;Lin et al. 2022;Zu et al. 2021;Sunayama et al. 2022;Zu et al. 2022).Using cluster weak lensing, Zu et al. (2021) found that the stellar mass of the central galaxies is a good observational proxy for concentration.Combining with the b measurements from cluster-galaxy cross-correlations, Zu et al. (2022) detected a possible evidence of the cluster assembly bias, thereby demonstrating a viable path to accurate HAB measurements using future cluster surveys.In particular, ΛCDM simulations predict that the cluster-mass haloes with higher concentrations have a lower bias than their low-concentration counterparts.Dalal et al. (2008, here-after referred to as D08) explained this phenomenon using the peak formalism (Kaiser 1984;Bardeen et al. 1986;Sheth & Tormen 1999), which predicts that the curvature of initial density peaks is correlated with the large-scale overdensity at fixed peak height.Using simulations, D08 demonstrated that the rare peaks with steep curvatures are more likely to collapse into clusters with high concentrations in low-density environments, and vice versa for those with shallow curvatures.In this paper, we follow D08 to trace haloes back to their initial density peaks and employ their peak curvature as a key indicator of HAB.
Although the presence of massive neutrinos delays the collapse of initial density peaks into haloes, the collapse time remains accurately predicted by the excursion set theory (Press & Schechter 1974;Bond et al. 1991) if the variance of the CDM (assuming baryons follow the CDM) rather than the total mass field is considered (Ichiki & Takada 2012).Therefore, it is perhaps unsurprising that Lazeyras et al. (2021) did not find any impact of neutrinos on HAB.However, Yu et al. (2017) discovered that the mass of haloes in neutrinorich environments tends to be systematically higher than those in neutrino-poor regions, probably due to the higher capturing rate of the slow-moving neutrinos in neutrino-rich environments.If halo assembly histories are coherently modulated by neutrinos across different environments, it is reasonable to expect such a modulation to leave some discernible imprint on the HAB in simulations.To isolate the impact of neutrinos, we construct a catalogue of "twin" haloes that originated from the same initial density peaks but evolved separately in massless vs. massive neutrino cosmologies.Therefore, any difference of HAB between the two twin samples would provide a smoking gun evidence of the impact of neutrinos on HAB.
This paper is organized as follows.We describe the construction of our peak-based twin halo catalogue in Section 2 and present the HAB comparison using different indicators in Section 3. We summarize our results and look to the future in Section 4. Throughout this paper, we adopt a spherical overdensity-based halo definition so that the inner overdensity within halo radius  200 is 200 times the mean density of the Universe, and the halo mass  ℎ is defined as the total mass enclosed by  200 .

PEAK-BASED TWIN HALO CATALOGUE
We compare the HAB effects with and without massive neutrinos using the Quijote simulations (Villaescusa- Navarro et al. 2020), the same suite used by Lazeyras et al. (2021).Among the vast Quijote ensemble, we select one "Fiducial" set of 40 ΛCDM simulations with Σ  =0.0 eV, and another " +++  " set of 40 massive neutrino simulations with Σ  =0.4 eV, hereafter referred to as the 0.0 and 0.4 simulations, respectively.The 0.4 simulations bear the highest Σ  within the Quijote suite, so that any possible neutrino-induced HAB would exhibit the strongest signal.Both the CDM and neutrinos are particle-based, with 512 3 CDM particles in a comoving volume of 1 ℎ −3 Gpc 3 of each simulation and an additional 512 3 neutrino particles for each 0.4 run.The two sets of simulations have the same cosmological parameters, with Ω  =0.3175, Ω  =0.049, ℎ=0.6711,   =0.9624,  8 = 0.834, and =−1.More important, identical CDM initial conditions (with the same phases) at =127 are shared by each of the 40 pair of 0.0 and 0.4 simulations, so that the haloes found at =0 can always be mapped back to the initial density peaks in the =127 snapshots.We identify dark matter haloes in the =0 outputs of the 40 pairs of simulations using ROCKSTAR (Behroozi et al. 2013).We calculate halo concentration  using the maximum circular velocity method (Prada et al. 2012), and use the relative concentration ĉ, defined as as an indicator of HAB, where ⟨| ℎ ⟩ and   ( ℎ ) are the average concentration-mass relation and its 1− scatter in either the 0.0 or 0.4 simulations.
To construct a high-purity twin halo catalogue from the paired simulations, we proceed step by step as follows.For each halo with  ℎ ≥10 14 ℎ −1  ⊙ at =0 in the 0.0 simulations, we collect all its CDM particles within  200 , retrieve their initial positions at =127, and compute the barycentre of those particles in the initial density field.We then identify the highest peak within a search radius of  pk = 200 1/3  200 about the barycentre as the true density peak that collapsed into that halo.The peak identification is performed over the smoothed density field derived with Nbodykit (Hand et al. 2018).Since the exact density peak also exists in the initial condition of the paired 0.4 simulation, we collect all the CDM particles within  pk about the peak in the initial condition, and identify its matched counterpart as the halo that inherited the highest fraction of those particles at =0 in the 0.4 simulation.To make sure the matches are unique, we repeat the above procedures but start from the haloes in the 0.4 simulations.We only keep the 97% twins that are reproduced during the reversed search in our final twin catalogue, leaving 3458747 unique halo twins with  ℎ ≥10 14 ℎ −1  ⊙ from the 40 paired simulations.
For the density peak associated with each twins pair, we characterise its peak height by the halo mass  ℎ at =0 in the 0.0 simulation to avoid any confusion in the notations (as  is reserved for neutrinos rather than peak height).The peak curvature  is computed as where Δ cdm is the CDM overdensity difference between 0.75 pk and  pk , and Δ log  is the logarithmic differences between the enclosed CDM masses within the two radii.Note that our curvature definition differs from that of D08 in two aspects: 1) we compute  using the gradient of the profile of initial density peaks at =127, instead of the main branch of the halo merger trees used by D08; 2) Unlike D08, our curvature definition has an extra minus sign so that steeper peaks have larger positive values of .As a sanity check, Figure 1 shows the distribution of 1000 randomly selected haloes on the  ℎ vs.  plane, colour-coded by  according to the colourbar on the right, where  is the linear theory-extrapolated Lagrangian overdensity at = 0. Similar to the Fig. 1 in D08, while the overall distribution of  centred around the critical overdensity   =1.68, steeper peaks generally collapsed earlier and have higher values of  than the shallower ones.To remove the average halo mass dependence of  in our analysis, we introduce the relative peak curvature ŝ as an HAB indicator, defined as where ⟨| ℎ ⟩ is the average curvature as a function of halo mass.
Inspired by neutrino differential condensation discovered in Yu et al. (2017), we measure the neutrino-to-CDM ratio   for each halo in the 0.4 simulations.In particular, we define   as where Ω cdm (Ω  ) is the average CDM (neutrino) density of the simulation, and  cdm (  ) is the local CDM (neutrino) density averaged between =5−12 ℎ −1 Mpc around each halo.The minimum aperture of 5 ℎ −1 Mpc is for excluding the halo region; We have tested different maximum apertures between 10−20 ℎ −1 Mpc and verified that the impact on our results is neglible.The radial range is set so that it is below the neutrino free-streaming length  fs (≃20 ℎ −1 Mpc for Σ  =0.4 eV) and the distance scales that we adopt for calculating halo bias.The average neutrino-to-CDM ratio ⟨   | ℎ ⟩ deceases slightly with increasing halo mass (black dashed curve), with a mass-dependent scatter    ( ℎ ) (gray shaded band).Similar to Equation 4, we define the neutrino excess parameter   as another indicator of the HAB, so that Our peak-based twin halo catalogue has a unique advantage for isolating the subtle effect of massive neutrinos on HAB.Since each twins pair is tied to the same density peak in the initial condition of both 0.0 and 0.4 simulations, we regard them as the same physical object that merely exhibits different behaviors with the massive neutrinos on and off.Therefore, if we measure two separate HAB signals using the relative peak curvature ŝ as the indicator in both 0.0 and 0.4 simulations, the difference between b( ŝ|0.0)and b( ŝ|0.4)can only be induced by massive neutrinos.More important, we can now measure b(  |0.0), which is the HAB of haloes in the 0.0 simulations, but using their   in the 0.4 simulations as the HAB indicator.By comparing this HAB signal with that measured directly from the 0.4 simulations b(  |0.4),we will be able to ascertain the direct impact of massive neutrinos on HAB.
The left column of Figure 2 shows the initial CDM overdensity  ini cdm at =127 surrounding two representative density peaks with equal peak height, one with a shallow ( ŝ=−0.51; top) curvature and the other steep ( ŝ=0.40; bottom).In each inset panel, we compare the 1D overdensity profiles of the two peaks, with the profile of that shown in the main panel indicated by the solid curve.Each peak subsequently collapsed into the twin haloes in the 0.0 (second column) and 0.4 (third column) simulations, with their halo radii  200 and characteristic scales   indicated by the outer orange and inner lime circles, respectively.The distributions of CDM particles (black dots) are very similar between the two simulations, and the neutrino den- . 3D cross-correlation functions between haloes and CDM  hm (left), haloes and neutrinos  h (middle), and the ratio between the two correlation functions  h /  hm (right), for haloes divided into four quartiles by their relative peak curvature ŝ (listed in the top right corner of the left panel).Solid and dashed curves show the measurements from the 0.0 and 0.4 simulations, respectively.The bottom panels show the ratio between profiles of the individual quartiles to that of the overall population.The inset panels in the bottom left panel provide a zoom-in view of the difference between the 0.0 and 0.4 measurements above 15 ℎ −1 Mpc.
sity 1 +   roughly follows the CDM distribution on large scales in the 0.4 simulation.The rightmost column shows the neutrino-to-CDM density ratio   / cdm (normalised by Ω  /Ω cdm ) surrounding the halos formed in the 0.4 simulation.Each panel of Figure 2 corresponds to a region of 30 ℎ −1 Mpc×30 ℎ −1 Mpc size centred on the peak/halo, projected along a 12 ℎ −1 Mpc length in the -direction.
All the distributions are colour-coded according to the respective colourbars on top.As expected, the steep peak collapsed into haloes with higher concentration than the shallow one.The concentrations of haloes in the 0.4 simulation is in general higher than those in the 0.0 one (Brandbyge et al. 2010;Villaescusa-Navarro et al. 2013;Lazeyras et al. 2021).This overall concentration difference is removed from the HAB comparison by using ĉ.In the third column, the two haloes in the 0.4 simulation share very similar neutrino content within the halo radii (orange circles), but differ significantly in their   distributions in the outside.However, the two   / cdm distributions shown in the fourth column are very similar on scales below 5 ℎ −1 Mpc, while exhibiting strong discrepancies on scales above 5 ℎ −1 Mpc.Therefore, within the 0.4 simulation, the halo formed from the high-ŝ peak has a larger   =0.57than that from the low-ŝ one (−0.75).This signals the existence of a positive correlation between ŝ and   , hence a potential neutrino-induced HAB.

RESULTS
To ascertain the correlation between ŝ and   , we show the 3D cross-correlation functions between haloes and CDM ( hm ; left), haloes and massive neutrinos ( h ; middle), and the ratio between the two ( h / hm ; right), for haloes in different quartiles of ŝ in Figure 3.
In each panel, the top main panel shows the profile comparison between four quartiles in ŝ, while the bottom sub-panel highlights their differences using the ratio between each quartile and the overall profile of that mass bin.The uncertainty bands are computed by Jackknife resampling over the 40 independent simulations.Note that we have combined the results from all halo masses above log  ℎ >14 in Figure 3, and we have carefully checked that any trend with halo mass has been removed by our adoption of the relative quantities.
In the left panel of Figure 3, solid and dashed curves are the  hm measurements from the 0.0 and 0.4 simulations, respectively.The  hm profiles exhibit the typical one-and two-halo components that can be described using the combination of NFW (Navarro et al. 1997) profiles and biased versions of the matter-matter auto-correlation function  mm (Hayashi & White 2008;Zu et al. 2014).A strong HAB effect is present in both simulations, as shown by the bottom subpanel, where haloes collapsed from steep peaks (red and brown) are more concentrated within the halo radius and have lower relative clustering amplitude on scales above 15 ℎ −1 Mpc than their shallow-peak counterparts (blue and green).However, the HAB in the 0.4 (dashed) simulations is slightly weaker than that in the 0.0 ones (solid), in that the relative clustering on scales above 15 ℎ −1 Mpc of the steepest (shallowest) quartile is lower (higher) in 0.4 than in 0.0 simulations.This weakening of ŝ-indicated HAB in the 0.4 simulations, albeit very small in amplitude, is robustly detected at ∼1− level across all radial bins above 15 ℎ −1 Mpc (see the zoomin inset panels).This suggests that the HAB is likely modified by the presence of massive neutrinos.The middle panel of Figure 3 is similar to the left, but for the  h profiles measured from the 0.4 simulations.Consistent with the results from LoVerde & Zaldarriaga (2014), the profiles of the "neutrino haloes" are shallower than that of the CDM haloes.Interestingly, the  h profiles exhibit its own HAB phenomenon analogous to the CDM version -the steep density peaks attract more concentrated "neutrino haloes" that exhibit lower relative clustering of neutrinos on large scales than the shallow ones.This to our knowledge is the first detection of the "neutrino halo" assembly bias effect in simulations.
The right panel of Figure 3 shows the ratio profiles between  h and  hm in the 0.4 simulations.Within  200 , the ratio profile rises sharply as ≃ 3/2 , but levels off as ≃ 2/3 beyond  200 .Though not shown here, the ratio profile should reach a plateau of unity at scales ≫ fs ≃20 ℎ −1 Mpc where neutrinos follow CDM.More important, the bottom panel reveals that ŝ indeed correlates with the neutrino-to-CDM ratio on quasi-linear scales.In particular, haloes collapsed from the steepest (shallowest) quartile exhibit a ∼10% enhancement (deficit) in their neutrino-to-CDM ratio between  200 and  fs compared to the average population.Therefore, we expect a positive correlation between ŝ and the previously defined neutrino excess parameter   .Given that massive neutrinos seem to cause an impact on the ŝ-indicated HAB (left panel of Figure 3), such a ŝ-  correlation implies that the impact on the   -indicated HAB would likely be greater.
While being theoretically informative, neither ŝ nor   is an ob-

Figure 1 .
Figure 1.Distribution of 1000 haloes randomly drawn from the twin halo catalogue on the CDM Lagrangian overdensity  (linearly extrapolated to =0) vs. halo mass  ℎ plane, colour-coded by their peak curvature  according to the colourbar on the right.

Figure 2 .
Figure 2. Distributions of the initial CDM overdensity  ini cdm (leftmost), current-day CDM density (2nd from left), neutrino density (2nd from right), and neutrino-to-CDM density ratio   / cdm (normalised by Ω  /Ω cdm ; rightmost) around two pairs of twin haloes with the same peak height, one collapsed from a shallow ( ŝ=−0.51; top) density peak and the other steep ( ŝ=0.40; bottom).Each panel shows a slice of dimension 30 ℎ −1 Mpc×30 ℎ −1 Mpc×12 ℎ −1 Mpc centred on the peak/halo, with the key halo/peak properties annotated in the bottom left corner and the simulation/redshift information in the top left.The maps are colour-coded according to the respective colourbars on top.The inset panel in each of the leftmost panels shows the radial profiles of the CDM overdensity around the two peaks, with the solid curve corresponding to the  ini cdm distribution shown in the main panel.The inner lime and outer orange circles on the rest of the panels denote the scale radii and  200 of the haloes, respectively.