Cosmological neutrino N-body simulations of dark matter halo

The study of massive neutrinos and their interactions is a critical aspect of contemporary cosmology. Recent advances in parallel computation and high-performance computing provide new opportunities for accurately constraining Large-Scale Structures (LSS). In this paper, we introduce the TianNu cosmological N-body simulation during the co-evolution of massive neutrino and cold dark matter components via the CUBEP$^3$M code running on the supercomputer Tianhe-2 and TianNu's connected works. We start by analyzing $2.537\times10^7$ dark halos from the scientific data of TianNu simulation, and compare their angular momentum with the matched halos from neutrino-free TianZero, revealing a dependence of angular momentum modulus on neutrino injection at scales below 50 Mpc and around 10 Mpc.


Neutrinos
Despite being one of the fundamental particles in cosmology, our knowledge of neutrinos' features and interactions is limited. The absolute mass of neutrinos remains unknown, and its determination is a challenging problem that has important implications for both cosmology and particle physics.
The methods for measuring neutrino mass include three technical approaches: double beta decay, single beta decay, and cosmological observations. The first method strongly depends on a specific assumption that neutrinos are Majorana-type particles, meaning they are their own antiparticles. The CUORICONO [1,2] and recent CUPID-0 [3] experiments employed this approach. The second approach is more direct and does not depend on any preset model, only on the relativistic energy-momentum relation. KATRIN experiments, as one of the approach, [4,5] have made great strides, documented in [6][7][8][9]. The third approach involves relic neutrinos, which are the second most abundant particles in the universe after photons. They emerged soon after the Big Bang. Cosmologists constrain the mass hierarchy of neutrinos by observing their effects on the Cosmic Neutrino Background (CNB) [10,11], Cosmic Microwave Background (CMB) [12,13], and Large-Scale Structure (LSS) [14,15]. Many recent works investigate these effects, including haloscope experiments for detecting Cold Dark Matter (CDM) candidates (such as axions and dark photons) [16,17], measurements of absolute neutrino mass in [18] and DUNE [19] experiments, and constraints on the CDM parameter(s) [20,21]. All these attempts reported at least two neutrinos having non-zero masses, with * The first two authors contribute equally to this work.
Version July 29, 2023 submitted to Journal Not Specified https://www.mdpi.com/journal/notspecified arXiv:2307.14621v1 [astro-ph.CO] 27 Jul 2023 a lower mass limit of M ν ≡ ∑ M ν,i ≥ 0.05eV from oscillation experiment [22] and an upper limit of M ν ≤ 0.12eV from Planck's CMB observation [23]. Therefore, one can estimate the total mass of neutrinos by observing its cosmological effects and check whether their gravitational effects on the LSS evolution at given M ν match the current simulations in the non-linear regime, where the previous analytical consideration of structure growth become invalid. Generally, large neutrino masses will reduce cosmological (<1Mpc) clustering and suppress the power spectrum features of CDM [24].

N-body simulations
Particles with mass interact through gravity. However, no analytical expression exists for the phase space of these particles when their number exceeds 2. Cosmological numerical evolution requires a large number of particles (beyond one billion), and provides descriptions of their states as an alternative method, which is referred to as the N-body simulation [25,26]. It has a well-established history in cosmology since the earliest realization in the 1970s [27][28][29][30].
Cosmological N-body simulations can model non-Gaussian effects and provide statistical samples for analyzing the universe at different stages. The simulation evolves depending on the composition of the universe, which can be determined by the famous Friedmann equation. Most matter in the universe is massive CDM, so point particle approximation is used to investigate their mutual gravity during evolution. At present, gravitational interactions are calculated employing particle-mesh (PM), particle-particle (PP), PP plus PM (P 3 M) and tree algorithms in cosmological simulations with extremely huge scales [31].
Recently, the largest particle-amount N-body simulation, cosmo-π [32], was conducted by the supercomputer π 2.0 at the weak-scaling efficiency of 95%, containing ∼ 4.39 trillion CDM particles in a (3.2/h Gpc) 3 cube with the CUBE code [33]. The Quijote-PNG N-body simulations [34] were performed to investigate the information on primordial non-Gaussianity that can be obtained beyond today's perturbative regime. The MillenniumTNG project [35] combined the Millennium simulation and IllustrisTNG hydrodynamical galaxy formation model to study matter clustering and halo statistics. And the SIMBA simulation [36] found that baryonic matter has tight relations with specific angular momentum j for stars, and the dispersion of stellar j-M relation is largely driven by HI content.

TianZero and TianNu
To investigate the effect of massive neutrinos on LSS, massive neutrino particles were incorporated into the TianZero CDM particle simulation, resulting in the TianNu co-evolution simulation.
The TianNu simulation, conducted on the Tianhe-2 supercomputer, is among the largest cosmological N-body simulations to date, surpassing its previous cosmological neutrino simulations by an order of magnitude in scale. Several of its key indicators remain higher than those of later experiments [32]. The enormous amount of particles, large box size, and high force resolution of TianNu make it an exceptional simulation, facilitating the study of the co-evolution of neutrino-CDM particles and the detailed structure of LSS.
Moreover, the clustering effects of neutrinos and the density fluctuation of CDM are easily affected by the Poisson noise, which may cover the actual signals in the neutrino quantum fluctuations [37]. Plenty of particles are needed to reduce the Poisson noise in the simulation, so the Tianhe-2 supercomputer can reach sufficient accuracy.
The paper is structured as follows. In Section 2, we provide an overview of the TianNu simulation. Section 3 presents the results of the process flow of the TianNu scientific data. In Section 4 we discuss the simulated impact of neutrinos on CDM halos. In Section 5, we summarize the conclusions of this paper.

Numerical Methods
We will review the TianNu simulation first because the data we analyzed in this paper are the CDM-halo part of scientific data from the TianNu simulation, and the simulation was fully performed only once at the maximum computational scale of the Tianhe-2 supercomputer.

Parameter preset
Using the customized parallel N-body simulation code CUBEP 3 M installed on the Tianhe-2, we perform two simulations inside a cube with a side length of 1200Mpc/h. One of them is TianZero only containing 6912 3 CDM particle-groups (every group has M CDMz = 7 × 10 8 M ⊙ as the mass resolution), and the other is TianNu including 6912 3 massive CDM particle-groups (M CDMν = 6.9 × 10 8 M ⊙ ) and 13824 3 neutrino (m ν = 0.05eV) particle-groups (every group has M ν = 3 × 10 5 M ⊙ as the mass resolution to keep Ω ν /Ω M ≈ 0.3711%). The only operative difference is that massive neutrinos are injected in TianNu at redshift z=5 with Ω M fixed, while TianZero is entirely massiveneutrino free. Moreover, according to [38], the initial positions of massive neutrinos were generated by a Gaussian noise map, and their velocity was preset by a linear Zeldovich part and random thermal part. Two species of massless neutrinos are both added in TianNu and TianZero via the CLASS transfer function [39] for the cosmological background, and they had no contribution to the CDM halo mass. To avoid Poisson noise and sampling variance [40], the selected box size is relatively less contaminated by noise at the scale of 0.05 < k/(h/Mpc) < 0.5.
The two simulations have the same cosmological model and seed initial condition, imposing a flat Universe by assuming that Ω M + Ω Λ = 1, where the matter density parameter Ω M = 0.29 (therefore Ω ν ≈ 0.001076), dark energy density parameter Ω Λ = 0.71. Moreover, reduced Hubble constant, initial tilt and mass variances at 8Mpc/h are set to be h = 0.73, n s = 0.95, σ 8 = 0.84 [41]. With the Spherical Overdensity algorithm to find halos, we identify CDM halos with at least 100 particles. Further information regarding the simulations and halo finder is available in [42,43].

Code overview
We recommend interested readers to look at previous works [38,44] for a detailed analysis of the code structure and technical algorithms relevant to both the pure CDM and neutrino cases. Below, we provide a brief introduction to the CUBEP 3 M code The CUBEP 3 M code gets its name from the cube space, PP force and two levels of the PM algorithm. The particles are currently stored in a computer-generated cube, while in its predecessor (the PMFAST code [45]), they were stored in slabs. PP indicates the interaction between particles is considered in a short-range manner, such as between several adjacent fine grids. PM means we place particles on meshes of varying roughness and use a coarse grid instead of a group of particles as the minimum unit to calculate the macro and long-range interactions with the gravitational potential [25] more efficiency.
CUBEP 3 M is an open-source Fortran code that includes Message Passing Interface (MPI) for inter-node communication and Open Multi-Processing (OpenMP) for intranode multi-threading [25]. Once installed, the codes can be compiled and run from the root directory and its subdirectories.
The code generates a snapshot of the simulated universe by recording the positions and velocities of all particles at a given epoch of cosmic evolution (redshift z) and saving the data in an xv.dat file. Other expected output from the code includes statistical data on star clusters (halofind.dat), 2D projections of the 3D velocity field (proj_xy.dat, proj_yz.dat, proj_xz.dat), and the power spectrum (cic_power.dat), among others. These output files are required when compiling the corresponding code module.

Halo data and parameters
Knowing the halo parameter format is essential to understanding CUBEP 3 M and to conducting our simulation. Therefore, we give our parameter table in Figure 1. The raw halo data are the snapshot of redshift z=0.01. In Figure 1, we can see that the hpos represents the final halo position (x,y,z), mass_vir denotes the halo Virial mass, r_vir is the virial radius of the halo, x_mean and v_mean are the position (x,y,z) and velocity (x,y,z) of halo mass center, l_CM is halo angular momentum modulus (x,y,z), v2_wrt is the component-squared vector (x 2 ,y 2 ,z 2 ) of the v_mean, var_x is the variance (x,y,z) of mass center position. In the second row, the new x_mean_nu and v_mean_nu are the position (x,y,z) and velocity (x,y,z) of the halo mass center, while N_nu is the number density of neutrinos.
Through position and mass data, we can match halos in TianNu and TianZero under specific conditions, such as percent mass variation less than 8% and mass center position less than 0.08Mpc. Using angular momenta components data, we can study their moduli changes and direction shifts. With velocities data, we can explore the velocity distribution of the halo mass center.

Data Analysis
Previous operations provide many halo samples, and have the chance to discover if neutrinos have any significant impact on the halo.
The study finds 2.560 × 10 7 halos with spherical over-density algorithm [42]. We match the same halos in Figure 2 between TianNu and TianZero within a variation tolerance of 10% (8%) TianZero mass M z and 0.1Mpc (0.08Mpc) center distance respectively. To tighten the match standard, we finally adopt the 8% level and acquire 2.537 × 10 7 halos, corresponding to the 98.68% of 10% level and 97.33% of whole halos. Adopting a comparable approach to Bullock et al. [46], we obtain the angular momentum: After the above process, we obtain halo data, x_m_l.dat (x1, x2, x3, M, J1, J2, J3), and calculate the modulus of J, With the three components of angular momentum, we can compute the variations in its magnitude and direction.
From Figure 3, we can obtain that the median value of relative modulus variation is -0.3655% and the most (89.71%) halos have less angular momenta moduli, indicating that neutrino-injection will reduce most moduli. We plot the relation between the 1σ region of ∆J/J z and M z in Figure 4, where every point (such as median, upper and lower limits) corresponding to one mass is extracted from one thousand dark halos with close masses and different points do not include same halos. But from Figure 5, the 95.44% of directions shift less than 0.65 degrees, indicating no significant impact of neutrinos. And we also plot the relation between the 1σ region of θ and M z in Figure 6 with the same 1000-bin statistic in Figure 4.  Now, we want to explore the universality of neutrinos' effect on halos in different environments. To expedite the search process, we utilize a two-level division procedure. Initially, we partition the entire box into adjacent coarse cubes, each with a side-length of 50 Mpc. Subsequently, we examine every halo within a fine cell (concentric circle) of radius 5 Mpc. We choose the top 20 densest and 20 sparsest cubes, and select 200 halos from the most compact cell (hereafter referred to as compact halos) and 200 halos from the most sparse cell (hereafter referred to as scarce halos) in each cube. For clarity purposes, we only use the terms dense or sparse to describe (coarse) cubes and the terms compact or scarce to depict (fine) cells.
In Figure 7, we plot the averaged ∆j/j z (the variation of mass-reduce angular momentum modulus, and j = J/ M, M is the CDM-particle group number of each halo,  namely the mass divided by its mass resolution) for all halos in each selected cube, and present the mean and stand deviation of all averaged ∆j/j z . The mean ∆j/j z of dense cubes is −0.0067 (11), while for sparse cubes it is −0.0053 (13). Although their error bars do not touch the other mean values, they overlap with the other bar in most areas. So their difference exists but is very small. To provide more information about Figure 7, we list the cubes' statistics obtained from CDM halos in Tables 1 and 2.
In Figure 8, we plot the averaged ∆j/j z for only 200 compact and 200 scarce halos in each selected cube, and handle them as what we do in Figure 7. The mean ∆j/j z of compact and scarce cells in dense cubes are −0.0077(17) and −0.0056 (11), while for the two groups in sparse cubes they are −0.0051(16) and −0.0042 (13). The error bars in the dense cube are more exclusive and distinct. Moreover, in the compact halo cells in dense and sparse cubes, the two standard deviations of ∆j/j z are always wider than the other two cases in scarce cells, implying different intensity of inner interactions according to their local densities.
For comparison, we plot the distribution and median-1σ region of ∆j/j z in Figures  9 and 10. 70.06% samples exhibit a decrease in ∆j after massive-neutrino injection. The median value of ∆j/j z is -0.0057, which is close to the mean value of 20 sparsest cubes (-0.0056), the 200 scarce halos of the 20 densest cubes (-0.0051).

Discussion
From Figure 2, the count of distance is discrete. This comes from relatively limited two decimal places extracted from data, and more digits may alleviate the situation.
The space distribution of LSS is typically irregular, and cannot be adequately represented by simple geometric shapes such as a cube or sphere. Thus, our local-  density statistic of halos is rudimentary and inaccurate (although it is still valuable for simplifying computations). One possible approach to further studying the LSS distribution is via clustering algorithms, such as K-means, density clustering (e.g., DBSCAN), hierarchical clustering, etc. However, these new methods face the same challenges, i.e. defining a standard density and addressing the presence of critical halos just outside the analyzed region. Further, incorporating Kernel Density Estimation may provide a unique solution.
Few exceptionally massive halos will contribute importantly to the density statistics, since coarse cubes with more massive halos will be more easily selected. Therefore, many low-mass halos inside the fine cell of each massive halo will be chosen.
The average ∆j/j z for 20 dense and 20 sparse cubes can be compared by examining Figure 7. The error bars for the mean values of these measurements indicate similar interaction intensities on the angular momentum at a cube scale. Figure 8 shows two panels with different average values and standard deviations. The right panel shows error bars that do not exclude the other mean value, unlike its left panel which exhibits a more noticeable difference in values. But the four cases show some consistency that in higher environment density region, halos incline to have more notable ∆j/j z . Higher density usually comes from some rare but very massive dark halos, which attract more low-mass halos and neutrino groups, naturally resulting in more intensive gravitation. From our results, this may lead to a potentially more significant loss of total angular momentum modulus. The two mean ∆j/j z of 200 halos in sparse cubes (the right panel of Figure 8) are both bit lower than the one of all halos in sparse cubes (Figure 7). Considering the different sample amounts (∼800 and 200) and the mutual coverage of error bars in the 200-halo case, the discrepancy is still acceptable. In Figure 8, although the average halo mass of scarce halos in dense cubes is one order of magnitude less than that of compact halos in sparse cubes, the former's mean ∆j/j z is double that the latter's, suggesting that the halos' local density environment, rather than their masses, predominantly influences the magnitude of ∆j/j z . Generally, it is more proper to treat the density-dependent j-suppression as a statistical rule instead of a certain one.  We select to re-scale J according to its mass resolution (a single particle group) instead of the unit mass M ⊙ , because the group is the least simulated object and the unit groups are mass-different in the two simulations. It is more intuitive to compare the performance of these distinct units. More interestingly, what differences are there among the j-distribution of CDM halos, neutrinos and their combination, and are these three distributions in one co-evolution simulation consistent with the one from the singlecomponent simulation? Confined by limited accessible data, these possible discrepancies are beyond our present ability. Therefore, we advocate more simulations to further explore these questions and examine our conclusions.
The severe non-linear effects in large-scale structure formation make it difficult to provide an accurate theoretical prediction regarding the variance magnitudes of halos' ∆J/J z , θ, and ∆j/j z after massive-neutrino injection. For this reason, N-body simulations are being utilized in this area. Our results require more simulations to validate this finding. However, there have been other studies on the angular momentum modulus of pure CDM halos, aside from the case of CDM-neutrino co-evolution. Liao et al. [47] examined the specific angular momentum j(r, θ) of halos obtained through a high-resolution Bolshoi simulation, and the resulting universal j(r, θ) profile confirmed the rigid shell model. Additionally, Li et al. [48] investigated the spin transfer from collisionless CDM to gas during halo formation using the non-radiative cosmological simulation SURFS, providing more detailed results.
Besides our content, several other developments in TianNu warrant a brief review. Inman et al. [38] implemented neutrino particles into the cosmology N-body code CUBEP 3 M [25], which is designed for dark matter simulation as one of its applications, thereby enhancing its functionality. This implementation involved updating the methods used to calculate density and velocity fields.
Emberson et al. [43] optimized the primary TianNu code [38], by introducing MPI decomposition for long-range PM force, adding OpenMP parallelism for short-range PM and PP force, and introducing a new data compression method from 24 bytes per particle to 9 bytes. These upgrades ensured the successful execution of the TianNu simulation on the Tianhe-2 supercomputer.
Yu et al. [42] used TianNu simulation to reveal the differential neutrino condensation effect: halo properties like its mass function vary depending on the local neutrino relative abundance (neutrino/CDM), which can be used to infer massive neutrino mass. Massive halos can capture more neutrinos in neutrinos-rich regions than those in neutrino-poor regions.
Inman et al. [40] proposed that the relative motion between neutrinos and CDM leads to a large-scale dipole on the matter density field which can offer orthogonal constraints to break the degeneracies in the two-point statistics caused by galaxy bias and optical depth in cosmological neutrino mass measurement. The linear response dipole (from relative velocity displacement treated in Fourier space) and N-body dipole (regarded as the result of real space) exhibited a perfect agreement, indicating the possibility of their conception.
Qin et al. [41] discovered a local effect from massive neutrinos on the spatial distribution of CDM halos concerning the Delaunay Triangulation (DT) voids with TianNu halo data. Small voids typically locate in regions with more massive neutrinos, making the halos around them more susceptible. The flux of dense massive neutrinos enlarges the void and pushes the neighboring halos away from the void center. Consequently, small voids would become larger than the case without these neutrinos, and large voids will be smaller. This feature impacts the halo two-point correlation at ∼1Mpc/h and notably disrupts the function of numerical DT voids. It will play a vital role in current or future galaxy surveys on measurable neutrino impacts.
Apart from TianNu simulation and CUBEP 3 M code, Yu et al. [49] utilized a new information-optimized CUBE code [33], which greatly optimized the memory utilization by reducing the memory space per particle from 24 bytes to 6 bytes, and proposed a parity-odd gravitational effect from massive neutrinos, which uniquely torques the direction of the combined angular momentum field of the galaxies and halos. More importantly, this torque detection is free of contamination from linear perturbation and easily separated from other non-gravitational influences. A 21cm HI survey (like HSHS [50]) aiming at galaxies with z ≈ 1 and M ν =0.05eV, can reach a 5-σ significance of neutrino mass detection.
Furthermore, Kohji et al. [51] employed a combination of codes on the Fugaku supercomputer rather than a single n-body co-evolution program, in their simulation of the first Vlasov simulation of relic massive neutrinos and N-body simulation of CDM. The integration of the Vlasov equation within the 6-dimensional phase space was accomplished using a high-order Vlasov solver enhanced by a series of numerical methods and special instructions for specific processors. The simulation, which involved 400 trillion Vlasov grids and 330 billion-body computations, consistently reproduced the accurate nonlinear neutrino dynamics. Compared to TianNu, the simulation required only about 2 hours for time-to-solution with higher velocity resolution, earlier massiveneutrino injection, and an equivalent realization level in contrast to TianNu's 52 hours.

Conclusions
This paper presents the TianNu simulation and its key parameters, including halo data and associated scientific outcomes. Currently, a growing number of scholars are using CUBEP 3 M to investigate various aspects of neutrinos, cosmic tides and magnetic fields. We aim to leverage the existing data and analysis framework for future endeavors.
By incorporating massive neutrinos into the CUBEP 3 M N-body code, we generate fruitful scientific outcomes. These outcomes, particularly concerning the absolute neutrino mass and hierarchical open issues on the large-scale structure effects in our Universe, expand our understanding.
Our findings suggest that 89.71% of halo samples display a decreasing trend in the angular momentum modulus when massive neutrinos are added. Halos located in denser local environments with a diameter of approximately 10Mpc exhibit a more prominent descending trend, implying that halos situated within cosmic filaments display smaller moduli with the presence of massive neutrinos. Meanwhile, the mean value of the relative variation of mass-reduced modulus (∆j/j z ) for 200 compact halos in different dense cubes has a large standard deviation, suggesting that the decreasing effect from massive neutrinos is more statistical, facing the complicated and intensive interplay in over-dense regions. At last, we review some remarkable work related to TianNu.
We hope this paper will be beneficial for future N-body simulations when browsing the TianNu project, and even to further improve the code.