Measuring neutrino mass and asymmetry with matter pairwise velocities

Neutrinos are believed to be the most abundant fermions in the Universe, but their masses are unknown, except for being non-zero but much smaller than other fermions. Cosmological relic neutrinos could also have non-zero chemical potentials (or asymmetries). Using neutrino-involved N-body simulations, we investigate the neutrino effects on the matter pairwise velocity, which itself is an interesting probe of cosmology. We find that for light-halo ($[10^{11},10^{13}]\ M_\odot$) mean pairwise velocity, in the transition range ($[4,15]\ \mathrm{Mpc}$), the effects of neutrino masses overwhelm the effects of neutrino asymmetries, while in the two-halo-group range ($[25,50]\ \mathrm{Mpc}$), for both light and heavy haloes ($[10^{13},10^{15}]\ M_\odot$), the effects of neutrino asymmetries dominate, making it possible to disentangle the two effects. We provide fitting formulae to quantify the effects of neutrino mass and asymmetry on halo-halo pairwise velocities.


INTRODUCTION
Neutrino physics carries implications beyond the Standard Model of Particle Physics.Notably, neutrino oscillation experiments unequivocally establish that at least two of the neutrino mass eigenvalues are non-zero.However, it is not known whether neutrino masses follow the normal or inverted hierarchy.From neutrino flavour oscillation experiments, assuming the normal (inverted) hierarchy, we can infer that the sum of neutrino mass eigenvalues, denoted as   ≡    , is greater than or equal to 0.06 eV (0.10 eV) (Esteban et al. 2020).The Karlsruhe Tritium Neutrino (KATRIN) experiment, which measures the tritium beta decay spectrum, has established an upper bound on the effective electron anti-neutrino mass,  ν , at 0.8 eV (at a 90% CL) (Aker et al. 2021), or   ≲ 2.4 eV.Project 8 seeks to lower the limit on  ν to 0.04 eV (  ≲ 0.12 eV) through Cyclotron Radiation Emission Spectroscopy (CRES, Project 8 Collaboration et al. 2022).
Beyond laboratory-based investigations, cosmology is another promising arena for studying neutrino properties, as neutrinos constitute the most abundant fermions in the universe.Relic neutrinos are radiation-like in the early universe and delay the epoch of radiation-matter equality.Consequently, they imprint discernible signatures on the cosmic microwave background (CMB).In particular, an upper bound of   ≲ 0.26 eV (95% CL) has been presented based solely on the Planck CMB data, with a refinement ★ E-mail: 1155129240@link.cuhk.edu.hk of   ≲ 0.13 eV (95% CL) achieved with the baryon acoustic oscillation (BAO) data included (Planck Collaboration et al. 2020).
In a later cosmic epoch, relic neutrinos undergo a transition to become matter-like.Nevertheless, according to Shoji & Komatsu (2010), neutrinos do not cluster at scales smaller than their freestreaming length-a length inversely proportional to their masses.The high dispersion velocity inherent to neutrinos tends to wash out the formed structures, thereby imparting a notable imprint upon the large-scale structure (LSS).Combining local LSS surveys, such as the Sloan Digital Sky Survey (SDSS), Pantheon (type Ia supernovae) and Dark Energy Survey 3x2pt (DES, weak lensing) data, together with Planck CMB data mentioned above, previous works have derived an upper bound for   ranging from 0.09 to 0.115 eV (95% CL, depending on which specific data is used and how the analysis is done) (Alam et al. 2021;Di Valentino et al. 2021;Palanque-Delabrouille et al. 2020).
In the near future, several LSS surveys, such as DESI1 , Euclid2 , LSST3 , and SKA4 , among others, are poised to becoming invaluable tools for determining neutrino properties.While tools based on linear perturbation theory such as CLASS (Lesgourgues 2011) and CAMB (Lewis et al. 2000;Howlett et al. 2012) allow for qualitative understanding of neutrinos' effects on LSS, for a quanti-tative description of the fully non-linear evolution of LSS we have to employ N-body simulations.
One approach to incorporating neutrinos into N-body simulations involves treating them as an additional category of collisionless simulation particles, as elucidated in previous studies (Brandbyge et al. 2008;Viel et al. 2010).Moreover, neutrino particles possess initial velocities described by the Fermi-Dirac distribution (Brandbyge et al. 2008;Viel et al. 2010;Bird et al. 2012;Villaescusa-Navarro et al. 2014;Castorina et al. 2015;Inman et al. 2015;Villaescusa-Navarro et al. 2015;Adamek et al. 2017;Emberson et al. 2017;Banerjee et al. 2018).Owing to their minuscule masses and weak interaction, neutrinos exhibit substantial free-streaming behaviour, particularly at high redshifts.To accurately resolve the evolution of neutrinos in their 6D phase space, one needs to employ small time steps, which are computationally demanding.Moreover, neutrino particles are prone to multiple crossings of the simulation boundary, leading to their stochastic dispersion within the simulation box.Consequently, particle-based neutrino simulations inevitably face the shot-noise issue, as reflected in the Poisson noise on the neutrino power spectrum, scaling as ∼ 1/  , where   is the number of simulated neutrino particles.A direct approach to mitigate shot noise involves increasing   , as exemplified by the adoption of almost O (10 12 ) neutrino particles in the work by Emberson et al. (2017); Yu et al. (2017).Another approach involves regular sampling of neutrino particles' velocities (Banerjee et al. 2018), a method implemented in the Quĳote simulations (Villaescusa-Navarro et al. 2020).
Although particle-based simulations can capture the non-linear evolution of neutrinos in the late universe, current methods for shot noise reduction impose substantial computational and memory requirements.Thus, an alternative approach employs the particlemesh (PM) grid, where neutrinos are elaborated on the PM grid whenever long-range forces are calculated (Brandbyge & Hannestad 2009).Within this framework, neutrinos either evolve under a fluid Boltzmann hierarchy expansion (Viel et al. 2010;Hannestad et al. 2012;Archidiacono & Hannestad 2016;Banerjee & Dalal 2016;Dakin et al. 2019;Tram et al. 2019;Inman & Yu 2020;Chen et al. 2021a) or respond linearly to a non-linear dark matter field (Ali-Haïmoud & Bird 2013;Liu et al. 2018;McCarthy et al. 2018;Zeng et al. 2019;Partmann et al. 2020;Chen et al. 2021b;Wong & Chu 2022).In particular, Yoshikawa et al. (2020Yoshikawa et al. ( , 2021) ) directly integrate the 6D Vlasov equation of neutrinos on the grid.As there are no actual neutrino particles, grid-based neutrino-involved simulations do not have shot noise problems and also require less memory since there is no need to store the information of one additional type of particles.However, in grid-based neutrino-involved simulations, the tree force among neutrinos is neglected.Fortunately, this sacrifice is justifiable for two reasons: first, neutrinos are exceptionally light, limiting their over-density growth to around  (1), and second, the scales of interest significantly exceed the mesh resolution.
Moreover, a hybrid approach, initially explored by Brandbyge & Hannestad (2010), offers an opportunity to take advantage of both grid-based and particle-based methods.In this approach, neutrinos are classified as either 'fast' or 'slow', each evolving under grid-based and particle-based methods, respectively (Banerjee & Dalal 2016;Bird et al. 2018).Another hybrid technique, also adopted in the MillenniumTNG (Hernández-Aguayo et al. 2023) and FLAMINGO (Schaye et al. 2023) simulations, solely captures deviations from a background neutrino phase space model, such as the thermal Fermi-Dirac distribution, during the evolution (Elbers et al. 2021).
Building upon contemporary methodologies and the capacity to execute high-resolution neutrino-involved simulations, numerous investigations have explored neutrino effects in real space.These investigations include the matter power spectrum (Dakin et al. 2019;Zeng et al. 2019;Elbers et al. 2021), the marked power spectrum (Massara et al. 2021), the void size function (Bayer et al. 2021;Verza et al. 2022), the halo bias (Villaescusa-Navarro et al. 2014;Castorina et al. 2015), the halo mass function (especially at the high-mass end, see e.g.Brandbyge et al. 2010;Castorina et al. 2015;Adamek et al. 2017;Liu et al. 2018;Bayer et al. 2021), the halo merger tree (Liu et al. 2018;Wong & Chu 2022), and the cosmic neutral hydrogen (HI) distribution (Villaescusa-Navarro et al. 2015).However, the matter velocity field, which conveys half of the phase space information was less studied.Several previous investigations have studied the impact of neutrinos on the velocity power spectrum (Inman et al. 2015;Howlett et al. 2017;Zhou et al. 2022).In this paper, we focus on the matter velocity field, particularly the mean peculiar pairwise velocity  12 () and the mean peculiar pairwise velocity dispersion  12 (), which have been shown to be powerful probes in constraining cosmological models.
In this paper, we study not only the effects of   on  12 () and  12 () but also the impact of the neutrino chemical potential   , where  labels the neutrino mass eigenstates.Since the distribution function of relic neutrinos is frozen after decoupling, we define the degeneracy parameters as   ≡   /  , where   is the relic neutrino temperature.Traditionally, regardless of whether it pertains to earlyuniverse observations, such as the CMB, or late-universe non-linear constraints on neutrinos, the analysis has assumed   = 0 (Planck Collaboration et al. 2020;Dark Energy Survey and Kilo-Degree Survey Collaboration et al. 2023).However, recent Planck CMB data still allow   to assume values on the order of unity (Yeung et al. 2021).More importantly,   may alleviate the Hubble tension (Yeung et al. 2021;Barenboim et al. 2017a).Additionally,   play a non-negligible role in the late universe, impacting the matter power spectrum (Zeng et al. 2019) and the halo merger tree (Wong & Chu 2022), although exhibiting a high degree of degeneracy with   .In this study, we focus on investigating neutrino effects using a gridbased method, following the approach outlined by Ali-Haïmoud & Bird (2013) and Zeng et al. (2019).
The subsequent sections of this article are organized as follows.Section 2 gives a brief review of relic neutrinos, accompanied by an elucidation of our simulation framework.Section 3 discusses the theoretical underpinnings of pairwise velocity and pairwise veloc-ity dispersion.In Section 4, we present our simulation results and quantify the neutrino effects, and Section 5 summarizes our findings and discusses their implications.Appendices A and B present the overview of our neutrino-involved simulation method and simulation convergence tests, respectively.

NEUTRINO-INVOLVED SIMULATIONS
In this section, we will give a brief introduction to relic neutrinos and our neutrino-involved simulations.

Neutrinos
We assume that the relic neutrinos would be frozen in the Fermi-Dirac distribution after their decoupling from the primordial thermal bath.Consequently, the neutrino energy density is where we adopt the natural units  = ℏ = 1 and assume three mass eigenstates of neutrinos labelled by .Here,   = √︃  2 +  2  is the neutrino energy, and the two bracketed terms represent neutrinos and their corresponding anti-neutrinos.Furthermore,   is proportional to   .We therefore introduce the constant factor   ≡   /  , the degeneracy parameter.
Subsequently, we will perform simulations using an array of combinations of   and  2 to systematically investigate their impact on the LSS, particularly the matter pairwise velocity  12 () and its corresponding dispersion  12 ().

Simulations
Given that only   is of relevance, we consider three degenerate neutrino mass eigenvalues.We have modified Gadget-2 (Springel 2005;Zeng et al. 2019) to incorporate neutrino effects.Our current version is almost as fast as the pure ΛCDM simulation without neutrinos.The details of the integration of neutrino effects into Nbody simulations are provided in Appendix A. The initial conditions for the simulations are generated using 2LPTic (Crocce et al. 2006), and to be consistent, we have modified the Friedmann equations used in original 2LPTic to account for neutrino effects.The initial matter power spectra at redshift  = 99 are derived using CAMB (Lewis et al. 2000), with the neutrino asymmetry added (Yeung et al. 2021).
We have studied a total of 10 cases, which includes an even distribution of   and  2 values, coupled with the corresponding cosmological parameters as detailed in Table 1.To ensure consistency, for each combination of   and  2 , the associated cosmological parameters have been obtained by refitting the Planck CMB and BAO data using CosmoMC, with neutrino asymmetry added (Yeung et al. 2021), with the Planck 2018 plikHM_TTTEEE likelihood.We perform two distinct simulation sets, denoted as  1 and  2 , both with a cold dark matter (CDM) particle count of   = 1024 3 and PM grid number  grid = 1024 3 , but with different box sizes of  box = 1000 ℎ −1 Mpc ( 1 ) and 250 ℎ −1 Mpc ( 2 ), respectively.Here, ℎ is the dimensionless Hubble parameter at  = 0, i.e.  0 = 100ℎ km s −1 Mpc −1 .The CDM particle mass and the gravitational softening length are around 8 × 10 10 ℎ −1  ⊙ (1 × 10 9 ℎ −1  ⊙ ) and 24.4 ℎ −1 kpc (6.1 ℎ −1 kpc) for the simulation set  1 ( 2 ), respectively.The identification of dark matter haloes is performed using ROCKSTAR (Behroozi et al. 2013), which uses adaptive hierarchical refinement of the friends-of-friends groups in six phase-space dimensions and one temporal dimension.We only consider distinct host haloes with mass defined following the virial threshold (Bryan & Norman 1998), i.e.Δ vir = 102 with respect to the critical density  crit .Moreover, we follow the definition of halo centre and velocity used in ROCKSTAR (Behroozi et al. 2013).The halo centre is the mean value of particles in the inner subgroup of this halo, and the halo velocity is averaging particle velocity within the innermost 10% of the halo virial radius (Behroozi et al. 2013).Similar to 2LPTic, we modify the expansion rate in ROCKSTAR to account for the neutrino contribution.
In Figure 1, we present a 2D projected density field  from the simulation set  1 .The left, middle, and right columns, respectively, show the cases of   =  2 = 0 (A0 in Table 1),   = 0.24 eV,  2 = 0 (A7), and   = 0.24 eV,  2 = 0.8 (A9).There are broad similarities in LSS in the different columns.However, upon closer inspection, visible differences become apparent.In particular, the finite neutrino mass serves to smooth out structures, effectively delaying structure formation.On the contrary, the introduction of a finite asymmetry results in sharper structures and an acceleration of the structure formation.This degeneracy between the effects of   and  2 on structures has been investigated in the context of the matter power spectrum (Zeng et al. 2019) and the halo merger tree (Wong & Chu 2022).In the present study, we focus on the impact of neutrinos on the matter velocity fields, representing the other half of phase space.

VELOCITY FIELDS
In this section, we present the theoretical framework underlying the matter pairwise velocity  12 (), and its corresponding dispersion  12 ().

Pairwise velocity
The mean pairwise peculiar velocity  12 () is defined as where  ≡  1 −  2 represents the distance vector between a pair of particles in the comoving frame, and  1,2 is the peculiar velocity of the respective particle.The average ⟨...⟩ is computed over all pairs separated by  ≡ ||, and  12 () depends only on  due to isotropy.The mean pairwise velocity  12 () characterizes the peculiar velocity difference between pairs of particles along their line of separation, revealing their tendency to approach each other.It is evident from Eq. ( 2) that when a pair of particles move towards each other,  12 < 0 (typically due to gravitational attraction).At large separation, particles have minimal correlation, so  12 () tends towards zero.Conversely, for small separations, the pair of particles is either bound together or in a state of stable clustering.In this case, their mean relative velocity, denoted as  12 in physical coordinates, approaches zero, implying (Mo et al. 2010, §6) where  () represents the Hubble parameter and  the scale factor.
In accordance with the pair conservation equation (Mo et al. 2010, §6), we have where  (, ) denotes the average number of particles within a comoving distance  of particle 2. The second term in the equation represents the average particle flux per unit time passing out of the spherical shell with radius  centred at particle 2. Here, n() is the mean particle number density at time , and Ξ(, ) is the two-point correlation function.By definition, We can simplify Eq. ( 4) by utilizing the conservation of the total No. number of particles, n() 3 = constant, to obtain where the averaged correlation function Ξ(, ) is defined as Thus, Eq. ( 6) reveals that  12 () encompasses not only the information of the real space two-point correlation function Ξ() but also its temporal evolution.

Pairwise velocity dispersion
Similarly, the mean pairwise peculiar velocity dispersion is defined as By integrating over momenta of the second Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) equation and following a series of simplifications, we arrive at (Peebles 2020, §72) where  13 ≡  1 −  3 and 12 ≡  1 −  2 . (1, 2, 3) is the threepoint correlation function and  the particle mass.This equation describes how velocity dispersion and gravity affect the evolution of Ξ.The third to last and final two terms are the gravitation from the particle pair and all neighboring particles, respectively.The mean pairwise peculiar velocity dispersion  12 bears a direct relation to the term ⟨  12   12 ⟩ within Eq. ( 10).This quantity can be expressed as where the 2⟨ 2 1 ⟩/3 is due to the uncorrelated motion between the particles, while Π and Σ characterize the effects of correlated motions along the components of dispersion parallel and perpendicular to the separation vector , respectively.Thus, when projecting Eq. ( 11) on the separation line, we obtain  2 12 = 2 3 ⟨ 2 1 ⟩ + Π.At small separation scales, the system assumes a statistically static state, allowing for an estimation of  12 () via the Cosmic Virial Theorem (CVT) derived from Eq. ( 10) under the condition  ≫ Ξ ≫ 1 (Mo et al. 2010, §6;Peebles 2020, §75): Therefore, the CVT serves as a pivotal link between the mean relative velocity squared (representing kinetic energy) and the mean gravitational potential energy originating from all neighboring pairs.When  is large, where the factor 2 arises from the pairing, and 1/3 signifies that only the direction along the line of separation is considered.In parallel, ⟨ 2 1 ⟩ can be estimated using the Cosmic Energy Equation (CEE; Mo et al. 2010, §6;Peebles 2020, §74;Mo et al. 1997), which stems from the integration of the first BBGKY equation over the momentum space, where I 2 ≡ ∫ ∞ 0 Ξ(, ).Integrating the above equation yields Notably, the CEE establishes a connection between the meansquared peculiar velocity (representing kinetic energy) and its associated mean gravitational potential energy.In essence, the CEE embodies the conservation of Newtonian energy within the expanding universe.
As shown in Zeng et al. (2019); Wong & Chu (2022), the neutrino mass and asymmetry affect the large-scale matter distribution and structure growth rates, based on the theoretical framework outlined above, we expect that they will also affect  12 and  12 .In the following section, we will explicitly test this theoretical framework in the neutrino case by comparing the theoretical calculations with the simulation results.

RESULTS
In this section, we will discuss our results of  12 () and  12 (), both at  = 0, obtained with the data sets  1 ( box = 1000 ℎ −1 Mpc) and  2 ( box = 250 ℎ −1 Mpc), as described in Section 2. There are three effects of relic neutrinos: the modified expansion history, the neutrino free-streaming effects, and the refitted cosmological parameters.According to Zeng et al. (2019), the modified expansion history has a negligible, sub-percentage impact on the matter power spectrum.On the other hand, the effects from neutrino free-streaming and the refitted cosmological parameters are more important.
In Appendix B3, we investigate the impact of cosmic variance (i.e.different random seeds to generate initial conditions) on the particle-particle (halo-halo) pairwise velocity,  pp ( hh ), and its associated dispersion,  pp ( hh ).We find that even though there are fluctuations in the magnitudes of  pp/hh and  pp/hh , the effect of neutrinos remains robust and unaffected by cosmic variance.Our results presented below are based on simulations using the same random seed.

Particle-particle case
As demonstrated in Appendix B1, it is accurate enough to randomly select 1% of the total CDM particles to calculate the particle-particle pairwise velocity ( pp ) and dispersion ( pp ), with uncertainties staying within 0.5% when  ≳ 4 Mpc.Additionally, from Appendix B2, the neutrino effects based on simulation sets  1 and  2 are comparable.Therefore, we will only present  pp and  pp using 1% of randomly chosen CDM particles from  1 .
Figure 2 displays  pp ( pp ) in the top (bottom) row.The subpanels show the percentage fractional deviations from the reference scenario with massless neutrinos (A0 in Table 1).The first column keeps  2 = 0 while changing   , and the second column varies  2 while keeping   = 0.06 eV.
The dashed lines in the top row of the subpanel are derived from the linear approximation of  pp given by Eq. ( 8), where Ξ is obtained by the inverse Fourier transformation (IFT) of the linear matter power spectrum at  = 0 (generated from CAMB) with the same cosmological parameters listed in Table 1.Our analysis shows that the simulation result (the solid black line) is in close agreement with the linear approximation (the dashed black line) for large  in the upper left subpanel.Additionally, in the lower subpanels of the top row, the neutrino effects from simulations (solid coloured lines) agree with the corresponding linear approximations (dashed coloured lines) at large .The linear theory consistently underestimates the behaviour of  pp and fails to provide meaningful insight into the non-linear regime, which encompasses nearly all of the  values depicted in Figure 2.
The coloured stars in the lower left panel of Figure 2 represent the results obtained through the CEE, as expressed in Eq. ( 13).Here, Ξ in I 2 (i.e. in Eq. ( 15)) is the IFT of the non-linear halofit matter power spectrum (Bird et al. 2012) (chopping  lower than 2π/ box to account for the finite box size effect).Due to the lack of a halofit model that takes into account the finite asymmetry  2 , we only include the cases with  2 fixed at zero.The lower left panel shows that the values of  pp calculated from CEE are in close agreement with our simulation results.In particular, the effects of neutrinos, as shown in the lower subpanel, are in good agreement with the predictions of the CEE.
The upper panels of Figure 2 show that the minimum of  pp is at  ∼ 4.5 Mpc, which is the transition scale between the onehalo and two-halo regions.As  increases or decreases, the pairwise velocity tends to zero, which is consistent with our theoretical predictions.At large separations, CDM particles move randomly with little correlation, whereas at smaller separations, they tend to form stable clusters.
In the first column of Figure 2, it is clear that as   increases, both the magnitudes of  pp and  pp decrease.This is in line with the idea that neutrinos tend to smooth out LSS.Additionally, when comparing the left and right columns of the figure, it is evident that   and  2 have opposite effects on  pp and  pp , except for  pp at large .In this case, a larger  2 reduces  pp , same as increasing neutrino masses.In summary,  2 increases the magnitude of  pp in the region between the one-halo and two-halo regimes and smooths out  pp in the two-halo region.However, interestingly, the effects of  2 on  pp remain in the opposite direction to   throughout the range we show.
To conduct a quantitative exploration analysis of these phenomena, we introduce the ratio of  pp , denoted as   pp , where  pp (0, 0, ) is from the fiducial massless scenario (A0 in Table 1).Subsequently, we can define the averaged value of   pp (  ,  2 , ) within a certain  interval   as Δ  p  (  ,  2 ), where  denotes different  intervals.
The effects of neutrinos are seen mainly at  beyond the onehalo region.Also, as Appendix B2 shows, the effects of neutrinos do not converge within the one-halo range under the current simulation M ν = 0.06 eV Figure 2. Particle-particle mean pairwise peculiar velocity ( pp , upper panels) and dispersion ( pp , lower panels), respectively, for  2 = 0 and different   (left panels), as well as   = 0.06 eV and different  2 (right panels).The subpanels below each graph show the percentage fractional deviations between the coloured cases and the reference case with massless neutrinos (A0 in Table 1).The upper panels show the linear approximation in dashed lines, Eq. ( 8), and the coloured stars in the lower left plot show the results from the CEE, as in Eq. ( 13).
resolution.Therefore, we consider two different  intervals: Furthermore, we apply the two-variable linear regression on Δ  p  (  , ), where   p  and   p  are two fitting parameters, and  ′  ≡   /0.06 eV denotes the dimensionless neutrino mass.
Similarly, in the case of  pp , we can define Subsequently, we employ the same linear regression approach for The results of the fitting are also shown in Table 2.In the transition region [4, 15] Mpc, when  ′  ∼  (1) and  2 ∼  (1), the neutrino mass and asymmetry have similar effects on  pp (  p1 (%) = −1.538± 0.006 and   p1 (%) = 1.569 ± 0.032), but with opposite signs.In the two-halo region [25, 50] Mpc, the effect of  2 is slightly more prominent than that of   (  p2 (%) = −1.254± 0.001 and   p2 (%) = 1.775 ± 0.004).Additionally, the effects of  2 on  pp remain consistent in direction for both the transition and two-halo regions, unlike the case for  pp .The upper left and right panels of Figure 3 show the fitting results for Δ  p1 and Δ  p2 .

Halo-halo case
The simulation sets  1 and  2 have different halo mass distributions.We choose host haloes from each set that fall within two particular virial mass intervals: [10 11 , 10 13 ]  ⊙ (referred to as "light haloes") from  2 , and [10 13 , 10 15 ]  ⊙ (referred to as "heavy haloes") from  1 .This selection yields a population of around 250, 000 haloes for each case, with each halo composing of at least 200 CDM particles.The upper panels of Figure 4 show that  hh for light range [Mpc]  haloes (dashed lines, mass range [10 11 , 10 13 ]  ⊙ ) is consistently less negative than for heavy haloes (solid lines, mass range [10 13 , 10 15 ]  ⊙ ).This may be due to two reasons: firstly, the stronger gravitational interaction among heavier haloes, secondly, halo bias, i.e. heavier haloes tend to be located in denser regions, resulting in larger magnitudes of pairwise velocities.The same pattern is visible in the lower panels of Figure 4, where  hh for light haloes is smaller than for heavy haloes.Furthermore, it is evident from the upper left panel of Figure 4 that, for two different ranges of halo mass, the minimum points of  hh (), denoted as ( min ,  min ), show considerable differences.For light haloes,  min ∼ 8 Mpc, which is the transition radius between the regions of one halo group and two halo groups.This is similar to the behaviour observed for  pp , although the transition radius is larger for haloes than for particles.In the vicinity of  min for light haloes, an increase in neutrino mass leads to an increase in  min , while  min is much less sensitive to  2 .These trends are also visible in  hh (dashed lines in the lower panels).
Additionally, from the subpanels of Figure 4, it is evident that the pairwise velocity of heavy and light haloes is more sensitive to  2 for  ∈ [30, 50] Mpc (the two-halo-group region).However, this sensitivity is reduced for small  centered around light haloes  min (the transition region).In particular,  hh remains rather insensitive to  2 on various  scales.For light haloes, these recognizable effects of the asymmetry, particularly in the two-halo-group region, and the directionality of its effects on  hh are consistent with the patterns observed in the context of  pp , as discussed in Section 4.1.
For heavy haloes (solid lines in Figure 4),  min ∼ 2 Mpc, where the effects of halo merging are dominant.This is only slightly larger than the average virial radius (approximately 1 Mpc) of these heavy haloes.Therefore, the result of multiple heavy haloes colliding and merging is observed.This behaviour is also reflected in  hh , which has a similar peak on the same scale.
In addition, the  hh for heavy haloes also show a transition point at approximately 8 Mpc, similar to light haloes.This suggests an intermediate zone between the one-halo-group and two-halogroup regions.When  is larger than this transition range,  hh and  hh have similar shapes as those of light haloes, but with larger magnitudes.In contrast, heavy haloes are less affected by neutrinos than light haloes, due to their greater gravitational attraction and resistance to smoothing.
We calculate the mean deviation of   hh from 1, denoted as where  = 1 (2) represents the interval  1 = [4, 15] Mpc ( 2 = [25, 50] Mpc).For light haloes, we study both intervals, which allows us to investigate the effects of neutrinos in the transition region and the two-halo-group region.For heavy haloes, we focus on  2 , since the  1 region has a high level of statistical noise, making it difficult to obtain reliable results with the current simulation resolution.
Similarly, we apply the two-variable linear regression model for both heavy and light haloes, which are expressed as Δ  h  =   h   ′  +   h   2 .This approach is also applied to  hh .The results of our analysis of Δ  h2 from light haloes and Δ  h2 from heavy haloes are shown in the lower left and right panels of Figure 3, respectively.The fitting coefficients are provided in Table 3.
Hence, ideally, we can independently measure  2 by accurately determining the halo mean pairwise peculiar velocity in the twohalo-group region for heavy haloes.This, combined with the  hh for light haloes in both the intermediate and two-halo regions, would allow us to gain a better understanding of the neutrino properties.This would also further validate the estimates of  2 derived from observations of heavy haloes.

CONCLUSIONS
In this study, we investigate and quantify the effects of neutrino mass   and asymmetry  2 on the particle-particle (halo-halo) mean pairwise peculiar velocity  pp ( hh ) and dispersion  pp ( hh ).
All our results are summarized in Tables 2 and 3.These effects are due to the combination of refitted cosmological parameters, neutrino free-streaming effects, and the modified expansion history of the universe.Given these results, the degeneracy of   and  2 found in the matter power spectrum (Zeng et al. 2019) and the halo merger tree (Wong & Chu 2022) can be disentangled by the halo-halo mean pairwise peculiar velocity in different regimes (transition and twohalo-group) for both heavy and light haloes.Our key conclusions are summarized below.
Moreover, in our upcoming paper, we measure the halo-halo mean pairwise peculiar velocity using galaxy observational data (Cosmicflows-4, Tully et al. 2023).We target to measure/constrain   and  2 .
Although we visually depict the one-halo region ( ≲ 2 Mpc) in our plots, we do not perform a quantitative analysis in this particular range due to the lack of resolution in our simulations.Further studies may benefit from neutrino-involved simulations with higher resolution, which could provide a more comprehensive understanding of the dynamics in this range.The lower subpanels, below  pp and  pp , show percentage fractional deviations with respect to the particle selection case 50%, and the gray-shaded regions show deviations within ±0.5%.
The code used in this paper is developed based on Zeng et al. (2019).The major difference is that we include one more term f ′  (1) in Eq. (A10) which is included in Ali-Haïmoud & Bird (2013).The technical improvement is that we have optimized the calculations as described above.

APPENDIX B: CONVERGENCE TESTS B1 The number of CDM particles
Given the large number of CDM particles (1024 3 ) used in our simulations, it is not feasible to include all of them in the calculation of particle-particle velocity and dispersion.Therefore, we investigate how our calculations are affected by the number of CDM particles used.Here, we take A1 in Table 1, with   = 0.06 eV and  = 0 in the simulation set  1 .
Figure B1 shows four different scenarios for  pp and  pp .These scenarios involve using 50% of the total CDM particles and three separate cases that involve randomly selecting 1% of the CDM particles.We find that the differences in fractional percentages between the 1% and 50% cases are kept within 0.5% in the region of interest ( ≳ 4 Mpc).
Therefore, in this study, we use 1% of the CDM particles for the  pp and  pp calculations.) and dispersion ( pp , lower panels), respectively.Solid and dashed lines are from simulation sets  1 and  2 , respectively.The various coloured lines correspond to  2 = 0 and different   (A1, A4, and A7 in Table 1).
The lower subpanels show the percentage fractional differences from the fiducial massless case (A0).

B2 Box size
In this subsection, we select 1% CDM particles for  pp and  pp calculations, as discussed in Appendix B1, to investigate how the effects of neutrinos are influenced by the box size.The analysis, as shown in Figure B2, keeps  2 = 0 while changing   .Comparing the results from the simulation sets  1 (solid lines) and  2 (dashed lines) shows that although the magnitudes of  pp and  pp are visibly different, the effects of neutrinos remain consistent beyond 4 Mpc, with only minor variations.These variations are mainly due to cosmic variance, with a slight contribution from the 1% CDM particle subsample used.This result is in line with our expectation, since  pp and  pp are statistical properties of CDM, and should not be too sensitive to the size of the simulation box.Similar results are seen when fixing   while varying  2 .

B3 Cosmic variance
In this paper, we use a random seed value of 23456 for our simulations.Additionally, we conduct a set of neutrino-involved simulations with a different random seed, 12345, to assess the stability of neutrino effects.As described in Appendices B1 and B2, only 1% of the CDM particles in the simulation set  1 are used to calculate  pp and  pp .
Figure B3 shows the results, with solid and dashed lines representing the random seeds 23456 and 12345, respectively.Comparing lines of the same colour, we find that the absolute values of  pp/hh ,  pp/hh may vary significantly, but the percentage fractional deviations (Δ/, Δ/) with respect to their respective fiducial massless cases (A0 in Table 1) remain stable between different random seeds within the range we are interested in ( ≳ 4 Mpc).Additionally, this stability is also observed when fixing   while varying  2 .
Consequently, in this paper, we will use seed 23456.
This paper has been typeset from a T E X/L A T E X file prepared by the author.1) presented in the subpanels below.The middle panels show the halo-halo mean pairwise peculiar velocity ( hh ) and the corresponding dispersion ( hh ) calculated from light haloes ([10 11 , 10 13 ]  ⊙ ), while the lower panels present  hh and  hh for heavy haloes ([10 13 , 10 15 ]  ⊙ ).Different random seeds are distinguished by solid and dashed lines.The colours in each panel represent different   with  2 = 0.

Figure 1 .
Figure 1.The 2D projected density distributions of CDM within the simulation set  1 at  = 0 are presented on various spatial scales.Each row displays a different scale, with details about the selected slices provided below.The left, middle, and right columns show the cases of   =  2 = 0 (A0 in Table1),   = 0.24 eV,  2 = 0 (A7), and   = 0.24 eV,  2 = 0.8 (A9), respectively.The second and third rows provide magnified views of the regions enclosed by white squares in the first and second rows, respectively.The last row has white circles that mark a host halo with its virial radius ( vir ), and the corresponding  vir and  vir are written on the upper-left corners.The density colour bars are placed to the right of each row.

Figure 3 .
Figure 3. Δ  p1 (upper left panel) and Δ  p2 (upper right panel) as a function of  2 for different   .The lower left and right panels show Δ  h2 (from light haloes) and Δ  h2 (from heavy haloes), respectively.The data points acquired from simulations, each with its own error bars, are fitted to the two-variable linear regression model (dashed lines), and the different colours show different   .

Figure 4 .
Figure 4. Halo-halo pairwise velocity  hh (upper panels) and dispersion  hh (lower panels), respectively, for heavy haloes ( 1 , solid lines) and light haloes ( 2 , dashed lines)  2 = 0 and different   (left panels), as well as   = 0.06 eV and different  2 (right panels).The first (light haloes) and second (heavy haloes) subpanels below each graph show the percentage fractional deviations between the massive and the corresponding massless neutrino cases.

Figure B1 .
Figure B1.Particle-particle mean pairwise peculiar velocity ( pp , upper panels) and dispersion ( pp , lower panels), respectively.The black line corresponds to a random selection of 50% of the total CDM particles, while various coloured lines are different sets of randomly chosen 1% particles.The lower subpanels, below  pp and  pp , show percentage fractional deviations with respect to the particle selection case 50%, and the gray-shaded regions show deviations within ±0.5%.

Figure B2 .
Figure B2.Particle-particle mean pairwise peculiar velocity ( pp , upper panels) and dispersion ( pp , lower panels), respectively.Solid and dashed lines are from simulation sets  1 and  2 , respectively.The various coloured lines correspond to  2 = 0 and different   (A1, A4, and A7 in Table1).The lower subpanels show the percentage fractional differences from the fiducial massless case (A0).

Figure B3 .
Figure B3.Particle-particle mean pairwise peculiar velocity ( pp , upper left panel) and dispersion ( pp , upper right panel) computed from the set  1 , respectively, with percentage fractional deviations (Δ/ and Δ/) relative to the reference massless case (A0 in Table1) presented in the subpanels below.The middle panels show the halo-halo mean pairwise peculiar velocity ( hh ) and the corresponding dispersion ( hh ) calculated from light haloes ([10 11 , 10 13 ]  ⊙ ), while the lower panels present  hh and  hh for heavy haloes ([10 13 , 10 15 ]  ⊙ ).Different random seeds are distinguished by solid and dashed lines.The colours in each panel represent different   with  2 = 0.

Table 1 .
Refitted cosmological parameters employed in simulations.  is the sum of the neutrino mass eigenvalues, while  2 is the neutrino asymmetry. 0 and Ω  are the Hubble parameter and the fractional energy density attributed to the   ℎ component, respectively, at redshift 0. Furthermore,   and   represent the spectral index and the scalar amplitude of primordial perturbations, respectively.

Table 2 .
quantity  Linear regression fitting results of Δ  p  and Δ  p  .

Table 3 .
Linear regression fitting results of Δ  h  and Δ  h  for different  vir .