Abstract
The standard model of cosmology predicts the existence of cosmic neutrino background in the present Universe. To detect cosmic relic neutrinos in the vicinity of the Earth, it is necessary to evaluate the gravitational clustering effects on relic neutrinos in the Milky Way. Here we introduce a reweighting technique in the Nonebody simulation method, so that a single simulation can yield neutrino density profiles for different neutrino masses and phase space distributions. In light of current experimental results that favor small neutrino masses, the neutrino number density contrast around the Earth is found to be almost proportional to the square of neutrino mass. The density contrastmass relation and the reweighting technique are useful for studying the phenomenology associated with the future detection of the cosmic neutrino background.
Introduction
The standard model of cosmology predicts that neutrinos were decoupled from the thermal bath when the temperature of the Universe was about 1 MeV. These relic neutrinos constitute the current cosmic neutrino background. Detecting cosmic relic neutrinos^{1}, as proposed in the upcoming PTOLEMY experiment^{2}, is thus a direct test of the standard model of cosmology, and can push our understanding of the Universe to its age of about one second. In order to detect them in the neighborhood of the Earth, a prerequisite would be to figure out the number density of relic neutrinos at our local environment. Although the standard model of cosmology does predict that the average number density of relic neutrinos in the current Universe is about 56 cm^{−3} for each flavor^{3}, more relic neutrinos can be accreted around the Earth, due to the fact that massive neutrinos suffer from the gravitational potential of both dark matter (DM) and baryonic matter in the Milky Way (MW). Investigating the gravitational clustering of relic neutrinos is thus a necessary step towards interpreting the results from the future detection of cosmic neutrino background.
Gravitational clustering effects are often studied numerically with the Nbody simulation method. However, to reach a resolution of ~8 kpc, the distance from the Earth to the galactic center of the MW, the Nbody simulation turns out to be computationally expensive^{4}. In 2004, a restricted but effective method called Nonebody simulation was proposed to evaluate the gravitational clustering effects of relic neutrinos^{5}. In contrast with the Nbody simulation, where all the interactions among particles are included, relic neutrinos in the Nonebody simulation are assumed to evolve under the gravitational potential of both DM and baryonic matter. The back reaction, i.e., the gravitational effects of neutrinos on the clustering of DM and baryonic matter, and the gravitational interactions among neutrinos are both considered to be negligible^{5,6}. This assumption works for the evolution of the Universe at a late stage (z ≲ 3 with z being the redshift), when the energy density of neutrinos is much smaller than that of DM^{3}. To implement the Nonebody simulation, one first divides the initial phase space of neutrinos into Nindependent chunks, and then evolves each chunk following a onebody motion in the gravitational potential generated by DM and baryonic matter. Assembling all the N chunks with their corresponding weights after the evolution yields the final phase space distribution of neutrinos.
In this work we introduce a reweighting technique in the Nonebody simulation, so that a single Nonebody simulation is sufficient to yield neutrino density profiles for different neutrino masses and phase space distributions. For small neutrino masses, we find that the neutrino number density contrast is almost proportional to the square of neutrino mass. The dependence of gravitational clustering effects on the phase space distribution is also investigated, followed by the implications of gravitational clustering effects on interpreting the results from the future detection of cosmic neutrino background.
Results
Normalized evolution equations
Here we adopt a generalized Navarro–Frenk–White (NFW) profile^{7} for the DM distribution in the MW, while for the baryonic matter distribution a spherically symmetric profile is also assumed for simplicity^{8}. See the Methods section for the details about the matter density profiles used in the numerical simulation.
Within the spherical gravitational potential ϕ(r), the onebody motion of a test particle with the mass m_{ ν } is confined to a plane, and obeys the following Hamiltonian equations^{9}
where a = 1/(1 + z) is the scale factor of the Universe, τ is the conformal time defined as dτ = dt/a(t), and \(p_r = am_\nu \dot r\) and \(\ell = am_\nu r^2\dot \theta\) are the canonical momenta conjugate to r and θ, respectively. Here the dot denotes the derivative with respect to τ, and (r, θ) are the polar coordinates in the comoving frame. The gravitational potential ϕ(r, τ) also evolves and its evolution is assumed to be independent of relic neutrinos in the Nonebody simulation. Because of spherical symmetry, \(\ell\) is a conserved quantity and we may ignore the motion in the θ direction. A key observation in developing the reweighting technique is to identify that the evolution equations in Eq. (1) can be written in a form independent of the neutrino mass m_{ ν }. Namely, with the normalized quantities \(u_r \equiv p_r{\mathrm{/}}m_\nu = a\dot r\) and \(u_\theta \equiv \ell {\mathrm{/}}m_\nu = ar^2\dot \theta\), the following normalized evolution equations can be derived^{5}
These normalized evolution equations are the ones implemented in our Nonebody simulation.
Reweighting technique
The essence of the reweighting technique is to let a test particle represent all the particles within a fixed interval (r_{ a }, u_{r,a}, u_{θ,a}) → (r_{ b }, u_{r,b}, u_{θ,b}) in the space spanned by (r, u_{ r }, u_{ θ }). The size of the interval does not depend on the mass or the phase space distribution of relic neutrinos. The dependences on these quantities arise when associating weight to the fixed interval (r_{ a }, u_{r,a}, u_{θ,a}) → (r_{ b }, u_{r,b}, u_{θ,b}). To illustrate that, we first introduce another variable set of (r, y, ψ), with y = p/T_{ν,0}. Here p denotes the magnitude of the canonical momentum, ψ is the direction of momentum with respect to the positive radial direction, and T_{ν,0} is the neutrino temperature at the present time. The transformations between (u_{ r }, u_{ θ }) and (y, ψ) are given by \(\psi = {\mathrm{tan}}^{  1}\left( {ru_r{\mathrm{/}}u_\theta } \right)\) and y = m_{ ν }u_{ r }/(cosψT_{ν,0}). Therefore, the fixed interval (r_{ a }, u_{r,a}, u_{θ,a}) → (r_{ b }, u_{r,b}, u_{θ,b}) corresponds to an interval (r_{ a }, y_{ a }, ψ_{ a }) → (r_{ b }, y_{ b }, ψ_{ b }), which will have varying lower limit (y_{ a }) and upper limit (y_{ b }) depending on the neutrino mass m_{ ν }. For the phase space interval (r_{ a }, y_{ a }, ψ_{ a }) → (r_{ b }, y_{ b }, ψ_{ b }), we obtain its associated weight dw as follows^{5}
where spherical symmetry is applied, and f(y) is the phase space distribution function. In the case of thermal relic neutrinos, f(y) follows the Fermi–Dirac form
The effect of neutrino masses reflects then in the lower and upper limits of y for a fixed interval in terms of (r, u_{ r }, u_{ θ }), while for different phase space distributions one simply uses the corresponding forms of f(y).
In practice, one still needs to perform a benchmark simulation with definite neutrino mass and phase space distribution. This benchmark simulation serves two purposes. First, the onebody evolutions of N test particles can be obtained. Second, in the benchmark simulation the initial phase space of relic neutrinos is discretized, and such a discretization would fix the interval (r_{ a }, u_{r,a}, u_{θ,a}) → (r_{ b }, u_{r,b}, u_{θ,b}) for each evolved sample. When switching to another neutrino mass or a different phase space distribution, on one hand we can reuse those onebody evolution results from the benchmark simulation, and on the other hand we can associate a new weight to each evolved sample according to Eq. (3). With this reweighting technique, we then do not need to rerun the Nonebody simulation for different neutrino masses and phase space distributions, so that lots of computing time can be saved.
Thermal relic neutrinos
In this work we consider a benchmark Nonebody simulation with m_{ ν } = 0.15 eV and thermal Fermi–Dirac phase space distribution. See the Methods section for the details about this benchmark simulation. The reweighting technique then enables us to obtain neutrino densitiy profiles n_{ ν }(r) for other neutrino masses and phase space distributions. Consider the thermal phase space distribution first. In Fig. 1a, we show the neutrino density contrast δ_{ ν }, defined as \(\delta _\nu \equiv n_\nu {\mathrm{/}}\bar n_\nu  1\) with \(\bar n_\nu\) being the average neutrino number density (δ_{ ν } corresponds to f_{c} − 1 in ref. ^{8}), as a function of the distance r to the galactic center of MW. The results of four different neutrino masses are shown, and the neutrino halos can extend up to a few mega parsecs. At the location of the Earth (r_{⊕} = 8 kpc) the neutrino density is enhanced by about 10% (115%) for the case of m_{ ν } = 0.05 (0.15) eV, due to the gravitational clustering effects.
For a fixed distance r, the relationship between the neutrino density contrast δ_{ ν } and m_{ ν } is displayed in Fig. 1b. We observe that for the three different distances of r, all the scatter points obtained from the Nonebody simulation can be well fitted by a powerlaw function of \(\delta _\nu \propto m_\nu ^\gamma\). The obtained exponents γ are around two for all cases, indicating that the linear approximation^{5,6} is appropriate in light of current cosmological constraints that favor small neutrino masses^{10}. Recall that in the Vlasov equation^{9} for the phase space distribution function there exists a term involving both the gravitational potential ϕ and the distribution function f. In the linear approximation, one approximates the distribution function f in that term with the corresponding distribution function f_{0} without the presence of gravitational potential, so that the modified Vlasov equation becomes linear in both ϕ and f. The underlying requirement for making the linear approximation is that the perturbed distribution function f should be close to the unperturbed one f_{0}, or the neutrino density contrast does not exceed greatly over order unity. For the neutrino masses considered in this work, according to Fig. 1 we find that the gravitational clustering effects are moderate so that the linear approximation works well here. However, if larger neutrino masses or heavier halo masses were considered, because of more enhanced gravitational clustering effects, the linear approximation would no longer be applicable, and the resulting powerlaw indices could have large deviations from two^{5}.
At the location of the Earth, the fitted powerlaw function for thermal relic neutrinos is given by
A similar powerlaw function^{11} was obtained for higher neutrino masses m_{ ν } ∈ [0.15, 0.6] eV, by fitting previous Nonebody simulation results^{5}. In the recent literature two benchmark values of m_{ ν } = 0.06 eV and 0.15 eV are also studied^{8}. Since the adopted DM and baryonic matter profiles in this work are the same as those in ref. ^{8}, we can directly compare our results with those in ref. ^{8}. After taking into account the uncertainties from discrete sampling, we find that both results agree with each other, validating the reweighting technique. With the above fitted relation we can obtain neutrino number densities for all neutrinos within the mass range of [0.04, 0.15] eV. For m_{ ν } < 0.04 eV, the clustering effects due to the gravitational potential of baryonic matter and DM in the MW are insignificant, namely, the neutrino density contrast is less than about 0.05 at the location of the Earth. As a result, other astrophysical uncertainties, such as the contribution from the Virgo cluster^{8,12}, may play more important role in predicting the local neutrino number densities. Furthermore, when m_{ ν } < 0.04 eV, the required energy resolution Δ (full width at half maximum of the Gaussian distribution) is estimated to be Δ ≃ 0.7m_{ ν } < 0.03 eV^{13}, which is beyond the reach of the current proposal of the PTOLEMY experiment^{2}. For these reasons, we choose not to consider the neutrino masses below 0.04 eV here.
New physics scenarios with nonthermal relic neutrinos
It is also interesting to consider new physics (NP) scenarios, in which some chiral states of neutrinos in the early Universe are nonthermal and possess a phase space distribution that is significantly deviated from the thermal Fermi–Dirac distribution. For illustration, we consider a fully degenerate phase space distribution^{14},
where y_{0} = 1.76 ensures the same average neutrino density as the thermal case. With the reweighting technique, we can also obtain the neutrino density profiles for this nonthermal case from the benchmark simulation. From Fig. 2, we observe that the neutrino contrast δ_{ ν } is about twice of that in the thermal case when m_{ ν } ≲ 0.1 eV. This is due to the fact that in the fully degenerate case more relic neutrinos reside in the low momentum states. For a given distance r, the relationship between δ_{ ν } and m_{ ν } can also be well described by a powerlaw function. The obtained exponents γ are also around two, so that the linear approximation^{5,6} can be applied to this fully degenerate case as well.
Implications on the direct detection of relic neutrinos
Finally, we discuss the impact of gravitational clustering effects on the direct detection of relic neutrinos in the vicinity of the Earth. For concreteness, we take the PTOLEMY proposal as an example. In the PTOLEMY proposal, relic neutrinos are captured by the betadecaying tritium nuclei ^{3}H via the process ν_{e} + ^{3}H → ^{3}He + e^{−}, and in the electron energy spetrum the signal events appear at a location that is 2m_{ ν } away from the betadecay end point (assuming three degenerate neutrino masses). Although there exist many experimental issues^{13,15}, such as the limited energy resolution and the large background from beta decay, we here only discuss the overall capture rate, assuming a total mass of tritium of 100 grams as given in the PTOLEMY proposal (see the Methods section for the details on calculating the capture rate). Since the capture rate depends on the nature of neutrinos and the underlying assumptions on the thermal history of neutrinos, we here focus on the case of Dirac neutrinos and discuss three typical physics scenarios. Note that in the following discussions the left(right)handed chiral states of neutrinos also mean the right(left)handed chiral states of antineutrinos.
The Standard Case refers to the scenario predicted by the standard model of cosmology, i.e., only the lefthanded chiral states of neutrinos were thermally produced in the early Universe, while for the righthanded chiral states of neutrinos their interactions with the particles in the thermal bath were so weak that their abundances can be safely neglected.
In the first NP scenario, i.e., NP Case I, the lefthanded chiral states of neutrinos have the same thermal history as Standard Case. However, because of some NP effects, the righthanded chiral states of neutrinos can also be thermally produced in the early Universe^{15,16}, except for a higher decoupling temperature depending on the current cosmological constraints. Considering the current constraint on the extra effective number of neutrino species ΔN_{eff} < 0.53 at the 95% confidence level (C.L.) by combining the Planck TT + lowP + BAO data sets with the helium abundance measurements^{10}, the number density of the righthanded chiral states is found to be at most 28% of that of the lefthanded chiral states^{15}. For illustration, we here set the number density of the righthanded chiral states to be 28% of that of the lefthanded chiral states, when neutrinos began freestreaming.
The second NP scenario, i.e., NP Case II, is almost the same as NP Case I, except that the righthanded chiral states are nonthermal and their phase space distribution is assumed to be the fullydegenerate form^{14}. To satisfy the same cosmological constraint as NP Case I, their number density can be 52% of that of the lefthanded chiral states^{17}, when neutrinos began freestreaming.
In Fig. 3a we show the capture rate Γ as a function of the total neutrino mass \({\sum} {\kern 1pt} m_\nu\) ≡ m_{1} + m_{2} + m_{3} in Standard Case. Here m_{ i }’s (for i = 1, 2, 3) are the three neutrino masses. The bestfit values of neutrino oscillation parameters from the latest global fit results^{18} are adopted, and two neutrino mass hierarchies, normal hierarchy (NH) (m_{1} < m_{2} < m_{3}) and inverted hierarchy (IH) (m_{3} < m_{1} < m_{2}), are distinguished. Because of gravitational clustering effects, we find that if taking \({\sum} {\kern 1pt} m_\nu\) = 0.23 eV, the 95% C.L. upper limit allowed by the Planck TT + lowP + lensing + BAO + JLA + H_{0} data sets^{10}, the capture rate can be enhanced by about 23% (31%) in the case of NH (IH), compared to the scenario without clustering. Because of the low target mass of tritium nuclei, distinguishing the neutrino mass hierarchies via the direct detection of relic neutrinos, given the current specifications of the PTOLEMY proposal^{2}, will be difficult. The capture rates in the NP Case I and II scenarios are given in Fig. 3b. Taking \({\sum} {\kern 1pt} m_\nu\) = 0.23 eV, the enhancement of capture rate due to gravitational clustering effects in NP Case I is the same as that in Standard Case, since in the two scenarios both the lefthanded and righthanded chiral states of neutrinos possess the same thermal Fermi–Dirac phase space distribution. However, in NP Case II we find that the rate of enhancement turns out to be higher, 33% (45%) in NH (IH), due to the fact that larger clustering effects are observed in the fullydegenerate phase space distribution. Moreover, due to the clustering effects the differences between the capture rates in NP Case I and II (solid curves vs. dashed curves) become more prominent for larger values of \({\sum} {\kern 1pt} m_\nu\). Therefore, gravitational clustering effects are useful for distinguishing NP scenarios involving different phase space distributions of relic neutrinos.
Discussion
In this work, when performing the Nonebody simulation, we did not consider the impact of other nearby galaxies (or clusters) on the clustering of relic neutrinos in the MW. In fact, it was estimated that the Virgo cluster may create a neutrino overdensity of the same order as that generated by the MW^{8,12}. Therefore, we also need to include the gravitational potential of the Virgo cluster into simulation. With both the MW and the Virgo cluster, the shape of the gravitational potential is no longer spherically symmetric. However, since each relic neutrino still follows a onebody motion under the total gravitational potential, the Nonebody simulation and the reweighting technique are still applicable, except that the normalized evolution equations in Eq. (2) and the associated weight for each test particle in Eq. (3) would be in more general forms without the spherical symmetry. Moreover, the currently used matter density profile of the MW still has some uncertainties^{8}, and thus the reweighting technique would also become useful, were we plan to update the Nonebody simulation with more accurate matter density profile in the future. In short, the reweighting technique simplifies the task of using the Nonebody simulation method to evaluate the gravitational clustering effects on relic neutrinos, and is useful for studying the phenomenology associated with the future detection of the cosmic neutrino background.
Methods
DM and baryonic matter density profiles in the MW
The DM density profile is taken to be a generalized NFW form^{8},
The parameter η is kept to be redshiftindependent, while r_{s}(z) and \({\cal N}(z)\) are evolved with redshift z. Here we adopt the bestfit values of \((\eta ,r_{\mathrm{s}}(0),{\cal N}(0))\) = (0.53, 20.29 kpc, 0.73) from a χ^{2}fit^{8} to data^{19}. The evolutions of r_{s}(z) and \({\cal N}(z)\) are dictated by the evolutions of the viral quantities Δ_{vir}(z) and c_{vir}(M_{vir}, z), which are defined as
Here M_{vir} and r_{vir}(z) are the virial mass and radius, respectively, and they are related by the following equation
Note that M_{vir} is redshiftindependent. The evolution of the critical density ρ_{crit}(z) is
where G is the gravitational constant, and H_{0} and Ω_{m,0} are the present values of the Hubble parameter and the matter density fraction, respectively. Here we adopt the bestfit values from the Planck data (H_{0}, Ω_{m,0}) = (67.27 km s^{−1} Mpc^{−1}, 0.3156)^{10}.
Given the three equations Eqs. (8), (9) and (10) and the quantities M_{vir}, Δ_{vir}(z) and c_{vir}(M_{vir}, z), we can obtain the redshift evolution for \(\left( {r_{\mathrm{s}},r_{{\mathrm{vir}}},{\cal N}} \right)\). The evolution of Δ_{vir}(z) is taken to be^{20}
where λ(z) = Ω_{m}(z) − 1 with Ω_{m}(z) being the matter density fraction at redshift z. With Δ_{vir}(0) and Eqs. (8), (10) and (11), we obtain M_{vir} = 3.76 × 10^{12}M_{⊙}, being M_{⊙} the solar mass. Subsequently, the evolution of r_{vir}(z) is also found from Eqs. (8) and (11). The evolution of r_{s}(z) is related to that of r_{vir}(z) via the concentration parameter c_{vir}(M_{vir}, z). The evolution of c_{vir}(M_{vir}, z) is assumed to be c_{vir}(M_{vir}, z) = \(\beta c_{{\mathrm{vir}}}^{{\mathrm{avg}}}\left( {M_{{\mathrm{vir}}},z} \right)\)^{8}, where \(c_{{\mathrm{vir}}}^{{\mathrm{avg}}}\left( {M_{{\mathrm{vir}}},z} \right)\) is taken from Nbody simulations^{20},
with the functions of A(z) and B(z) given by^{20}
The value of β = 2.09 is then obtained from \(c_{{\mathrm{vir}}}^{{\mathrm{avg}}}\left( {M_{{\mathrm{vir}}},0} \right)\) and c_{vir}(M_{vir}, 0) = r_{vir}(0)/r_{s}(0). Finally, with the obtained r_{s}(z), the evolution of \({\cal N}(z)\) can be found from Eq. (10).
For the baryonic matter density profile, in reality it has a bulge shape for the central region of radius ~5 kpc^{19} and extends to be a disc in the outer region. However, since the baryonic matter density is much smaller than the DM density in the disc region, we choose to ignore the disc part, and assume the overall profile of the baryonic matter density to be spherically symmetric. The actual baryonic matter density profile as a function of the radius at the present time is adopted from Fig. 2 in ref. ^{8}. The redshift dependence of baryonic matter profile is modeled with a normalization factor \({\cal N}_{\mathrm{b}}(z)\)^{8}, and the evolution of \({\cal N}_{\mathrm{b}}(z)\) is found by averaging over the evolution of the total stellar mass as obtained in eight MWsized simulated halos^{21}.
Benchmark Nonebody simulation
In the benchmark Nonebody simulation, we set m_{ ν } = 0.15 eV and consider relic neutrinos with a Fermi–Dirac phase space distribution. The normalized evolution equations in Eq. (2) are used, and the scale factor is evolved as
We utilize an iteration procedure to discretize the initial phase space into finer parts^{5,8}, so that a smooth neutrino density profile can be achieved up to ~ 5 kpc. Because of spherical symmetry, the set of variables to be discretized is (r, y, ψ). In general, we have r ∈ [0, ∞), y ∈ [0, ∞) and ψ ∈ [0, π). However, for obtaining the neutrino overdensity around the Earth, we can restrict the ranges for (r, y, ψ). Specifically, we consider two regions of r, an inner region R_{I}: r ∈ [0, 14.3 Mpc), and a outer region R_{O}: r ∈ [14.3 Mpc, r_{max}]. The value of r_{max} is taken to be the maximal radial distance that a neutrino with mass m_{ ν } can travel, given that it must arrive at the Earth at z = 0 and with a momentum p_{max}. Here p_{max} depends on the form of f(y), e.g., for the thermal case we can truncate y to y_{max} = 10, as the contribution from y > 10 is negligible. Although for m_{ ν } = 0.15 eV setting r_{max} = 150 Mpc is sufficient, we here take r_{max} = 350 Mpc in order to cover smaller neutrino masses.
In the inner region R_{I} we consider y ∈ [0, 10] and ψ ∈ [0, π). In the coarse scan we set the resolution of r to be Δr = 14 kpc, and divide the whole ranges of y and ψ into 120 parts evenly. After the coarse scan, we select the chunks that end up in the sphere of r < r_{7} = 7 Mpc, and further divide the ranges of r, y and ψ into another 2, 4 and 4 parts, respectively, for the fine scan.
For the outer region R_{O}, we first notice from the results of ref. ^{5} that when m_{ ν } ≲ 0.3 eV and M_{m} ≲ 10^{13}M_{⊙}, with M_{m} being the total mass of the gravitational source, the clustering effect of neutrinos is rather weak beyond the distance r_{7}. As a result, the initial direction of the momentum can be confined to ψ(r) ∈ [π–r_{7}/r), as the test particles with ψ(r) < π–r_{7}/r would never reach the location of the Earth. The range of y in R_{O} is also bounded as y(r) ∈ [y_{min}(r), y_{max}(r)]. The existence of y_{min} is due to the fact that a minimum velocity is always needed in order to reach the location of the Earth, while for y_{max} it is mainly because an escape velocity exists for a given gravitational potential. In the actual simulation, we obtain the dependence of y_{min} and y_{max} on r through a quick numerical scan. Specifically, in the quick scan we vary y within a wide range while having ψ ∈ [π–r_{7}/r), and then find out the allowed range of y that leads to the final value of r below r_{7}. The lower and upper limits of such an allowed range of y yield y_{min} and y_{max}, respectively. We repeat this kind of numerical scan for a few values of r, and then perform a parametric fit to obtain the dependence of y_{min} and y_{max} on r. Having found the ranges of y and ψ, we then perform a coarse scan with a resolution of Δr = 140 kpc and divisions of 120 parts for both y and ψ. After the coarse scan, we again select the chunks that end up within the sphere of r < r_{7}, and then perform two rounds of fine scan. In the first round of fine scan, we further divide the ranges of y and ψ into another five parts for the selected chunks from the coarse scan. In the second round of fine scan we focus on the chunks that end up within the sphere of r < 140 kpc from the first round of fine scan, and another divisions of 20 parts are applied to both y and ψ.
The last step of Nonebody simulation is to reconstruct the neutrino density profile from the discrete evolved samples. Adopting the kernel method^{22} and with spherical symmetry, the neutrino number density profile can be estimated as
with the kernel function K(r, r_{ i }, d) given by
Here i is the sample index, w_{ i } (r_{ i }) denotes the weight (final radial distance) of the sample i, and N is the total number of samples. The window size d plays a similar role as the bin size in the regular histogram. Therefore, too small values of d would lead to large shot noise, while too large values would smooth out the detailed structure of interest. In this work we set d = 8 kpc, and find that varying d by 3 kpc would modify the neutrino density n_{ ν }(r_{⊕}) by around 10%. Note that a similar level of uncertainty around 8% due to discrete sampling is also observed in ref. ^{8}, and within the uncertainties our simulated results for the case of m_{ ν } = 0.15 eV agree with those in ref. ^{8}.
Capture rate in the PTOLEMY experiment
The capture rate of relic neutrinos in the PTOLEMY experiment can be calculated as^{13}
where s_{ ν } denotes the helicity of neutrinos, j labels the neutrino mass eigenstate, and n_{ j }(s_{ ν }) is the local number density of neutrinos with mass m_{ j } and helicity s_{ ν }. In addition, the total number of tritium nuclei is represented by N_{T}, and \(\bar \sigma\) = 3.834 × 10^{−45} cm^{2} is the crosssection for the capture of relic neutrino on tritium nuclei. Since only the electron neutrinos can be captured by the tritium nuclei, we need to project each neutrino mass eigenstate into the electron neutrino flavor state, resulting in the factor \( {U_{{\mathrm{e}}j}} ^2\) in the above formula. The lepton mixing parameters \( {U_{{\mathrm{e}}j}} ^2\) are found by adopting the bestfit values of neutrino oscillation parameters^{18}, i.e., \(\left( {\left {U_{{\mathrm{e}}1}} \right^2,\left {U_{{\mathrm{e}}2}} \right^2,\left {U_{{\mathrm{e}}3}} \right^2} \right)\) = (0.665, 0.314, 0.022). The spindependent factor C(s_{ ν }) is given by C(s_{ ν }) = 1 − \(2s_\nu v_{\nu _j}\) with \(v_{\nu _j}\) being the velocity of neutrino with mass m_{ j }. For the total neutrino mass \({\sum} {\kern 1pt} m_\nu\) considered in Fig. 3, all neutrino mass eigenstates are nonrelativistic at the present time, so that we can set \(v_{\nu _j}\) ≃ 0 and C(s_{ ν }) ≃ 1. As a result, the last quantity that we need to specify in Eq. (18) is the local number density of neutrinos n_{ j }(s_{ ν }), which depends on the underlying physics scenarios and the gravitational clustering effects.
In all the three physics scenarios considered in this work, the nature of Dirac neutrinos is assumed. Therefore, lepton number is conserved and thus only neutrinos, not antineutrinos, can be captured. Regarding the neutrino helicity states, in the early Universe both lefthanded and righthanded chiral states are relativistic, and therefore they are lefthanded and righthanded helical states as well. In the evolution of the Universe, the helicites of neutrinos are preserved, so that at the present time the lefthanded and righthanded helical neutrino states correspond to the lefthanded and righthanded chiral states in the early Universe, respectively. In Standard Case, the standard model of cosmology predicts the average number density of lefthanded helical neutrino states to be n_{0} = 56 cm^{−3} for each mass eigenstate at the present time, while almost no existence of righthanded helical states. For a given value of the total neutrino mass \({\sum} {\kern 1pt} m_\nu\), we can obtain the three individual neutrino masses with the two masssquared differences from neutrino oscillation experiments. Here we adopt \({\mathrm{\Delta }}m_{21}^2 \equiv m_2^2  m_1^2\) = 7.56 × 10^{−5} eV^{2} and \(\left {{\mathrm{\Delta }}m_{31}^2} \right \equiv \left {m_3^2  m_1^2} \right\) = 2.55 × 10^{−3} eV^{2} from ref. ^{18}. Taking \({\sum} {\kern 1pt} m_\nu\) = 0.23 eV as an example, we obtain (m_{1}, m_{2}, m_{3}) = (0.0711 eV, 0.0717 eV, 0.087 eV) in NH. From the fitted function in Eq. (5) we find that the neutrino density contrasts δ_{ ν } are (0.22, 0.23, 0.35) for the three mass eigenstates, respectively, and therefore the corresponding number densities of lefthanded helical states around the Earth are 68.44, 68.64 and 75.53 cm^{−3}. Finally, from Eq. (18) we obtain the capture rate Γ = 4.97 yr^{−1} in this case. The calculation of capture rates in NP Case I and II can be performed similarly, except that in NP Case I and II there are additional contributions from the righthanded helical states of neutrinos, whose average number densities at the present time are taken to be 0.28n_{0} and 0.52n_{0} for each mass eigenstate in NP Case I and II, respectively. Gravitational clustering effects on the righthanded helical states are calculated in the same way as the lefthanded helical states, except that in NP Case II the fitted relation for the fullydegenerate phase space distribution should be used.
Data availability
The data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request.
References
 1.
Shvartsman, B. F., Braginsky, V. B., Gershtein, S. S., Zeldovich, Y. B. & Khlopov, M. Y. Possibility of detecting relic massive neutrinos. JETP Lett. 36, 277–279 (1982).
 2.
Betts, S. et al. Development of a relic neutrino detection experiment at PTOLEMY: Princeton Tritium Observatory for Light, EarlyUniverse, MassiveNeutrino Yield. Preprint at: http://arXiv.org/abs/1307.4738 (2013).
 3.
Lesgourgues, J. & Pastor, S. Massive neutrinos and cosmology. Phys. Rep. 429, 307–379 (2006).
 4.
Emberson, J. D. et al. Cosmological neutrino simulations at extreme scale. Res. Astron. Astrophys. 17, 085 (2017).
 5.
Ringwald, A. & Wong, Y. Y. Y. Gravitational clustering of relic neutrinos and implications for their detection. JCAP 0412, 005 (2004).
 6.
Singh, S. & Ma, C. P. Neutrino clustering in cold dark matter halos: Implications for ultrahighenergy cosmic rays. Phys. Rev. D 67, 023506 (2003).
 7.
Navarro, J. F., Frenk, C. S. & White, S. D. M. The structure of cold dark matter halos. Astrophys. J. 462, 563–575 (1996).
 8.
de Salas, P. F., Gariazzo, S., Lesgourgues, J. & Pastor, S. Calculation of the local density of relic neutrinos. JCAP 1709, 034 (2017).
 9.
Bertschinger, E. Cosmological dynamics: course 1. Preprint at: http://arXiv.org/abs/astroph/9503125 (1995).
 10.
Ade, P. A. R. et al. Planck 2015 results: XIII. Cosmological parameters. Astron. Astrophys. 594, A13 (2016).
 11.
Li, Y. F. & Xing, Z. Z. Captures of hot and warm sterile antineutrino dark matter on ECdecaying Ho163 nuclei. JCAP 1108, 006 (2011).
 12.
VillaescusaNavarro, F., MiraldaEscudé, J., PeñaGaray, C. & Quilis, V. Neutrino halos in clusters of galaxies and their weak lensing signature. JCAP 1106, 027 (2011).
 13.
Long, A. J., Lunardini, C. & Sabancilar, E. Detecting nonrelativistic cosmic neutrinos by capture on tritium: phenomenology and physics potential. JCAP 1408, 038 (2014).
 14.
Chen, M. C., Ratz, M. & Trautner, A. Nonthermal cosmic neutrino background. Phys. Rev. D 92, 123006 (2015).
 15.
Zhang, J. & Zhou, S. Relic righthanded Dirac neutrinos and implications for detection of cosmic neutrino background. Nucl. Phys. B 903, 211–225 (2016).
 16.
Horvat, R., Trampetic, J. & You, J. Spacetime deformation effect on the early Universe and the PTOLEMY experiment. Phys. Lett. B 772, 130–135 (2017).
 17.
Huang, G. Y. & Zhou, S. Discriminating between thermal and nonthermal cosmic relic neutrinos through an annual modulation at PTOLEMY. Phys. Rev. D 94, 116009 (2016).
 18.
de Salas, P. F., Forero, D. V., Ternes, C. A., Tortola, M. & Valle, J. W. F. Status of neutrino oscillations 2017. Preprint at: http://arXiv.org/abs/1708.01186 (2017).
 19.
Pato, M. & Iocco, F. The dark matter profile of the Milky Way: a nonparametric reconstruction. Astrophys. J. 803, L3 (2015).
 20.
Dutton, A. A. & Macciò, A. V. Cold dark matter haloes in the Planck era: evolution of structural parameters for Einasto and NFW profiles. Mon. Not. Roy. Astron. Soc. 441, 3359–3374 (2014).
 21.
Marinacci, F., Pakmor, R. & Springel, V. The formation of disc galaxies in high resolution movingmesh cosmological simulations. Mon. Not. Roy. Astron. Soc. 437, 1750–1775 (2014).
 22.
Merritt, D. & Tremblay, B. Nonparametric estimation of density profiles. Astron. J. 108, 514–537 (1994).
Acknowledgements
We are indebted to Zhizhong Xing for suggesting the topic and useful discussions. We also thank Ran Ding, Guoyuan Huang, Yandong Liu and Shun Zhou for useful discussions. This work was supported by the National Natural Science Foundation of China (Grants No. 11522540 and No. 11690021), the National Program for Support of Topnotch Young Professionals, the Provincial Department of Education of Liaoning (Grant No. L2012087), and the China Postdoctoral Science Foundation (Grant No. 2017M610008).
Author information
Affiliations
Contributions
J.Z. and X.Z. jointly initiated the project. Both authors have contributed to this work.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Additional information
Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Electronic supplementary material
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Zhang, J., Zhang, X. Gravitational clustering of cosmic relic neutrinos in the Milky Way. Nat Commun 9, 1833 (2018). https://doi.org/10.1038/s4146701804264y
Received:
Accepted:
Published:
Further reading

Exploring neutrino mass and mass hierarchy in interacting dark energy models
Science China Physics, Mechanics & Astronomy (2020)

Dark energy versus modified gravity: Impacts on measuring neutrino mass
Science China Physics, Mechanics & Astronomy (2020)
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.