Neutrino clustering in the Milky Way and beyond
P. Mertsch, G. Parimbelli, P.F. de Salas, S. Gariazzo, J. Lesgourgues, S. Pastor
TTTK-19-44
Prepared for submission to JCAP
Neutrino clustering in the Milky Wayand beyond
P. Mertsch, a G. Parimbelli, b,c,d
P.F. de Salas, e S. Gariazzo, f J.Lesgourgues a and S. Pastor f a Institute for Theoretical Particle Physics and Cosmology (TTK)RWTH Aachen University, D-52056 Aachen, Germany b SISSA - International School for Advanced Studies, Via Bonomea 265, 34136 Trieste, Italy c INFN - National Institute for Nuclear Physics, Via Valerio 2, 34127 Trieste, Italy d IFPU - Institute for Fundamental Physics of the Universe, Via Beirut 2, 34151 Trieste,Italy e The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, StockholmUniversity, SE-106 91 Stockholm, Sweden f Instituto de F´ısica Corpuscular (CSIC-Universitat de Val`encia)Parc Cient´ıfic UV, C/ Catedr´atico Jos´e Beltr´an, 2, E-46980 Paterna (Valencia), SpainE-mail: [email protected], [email protected],[email protected], gariazzo@ific.uv.es, [email protected],pastor@ific.uv.es
Abstract.
The standard cosmological model predicts the existence of a Cosmic NeutrinoBackground, which has not yet been observed directly. Some experiments aiming at itsdetection are currently under development, despite the tiny kinetic energy of the cosmologicalrelic neutrinos, which makes this task incredibly challenging. Since massive neutrinos areattracted by the gravitational potential of our Galaxy, they can cluster locally. Neutrinosshould be more abundant at the Earth position than at an average point in the Universe.This fact may enhance the expected event rate in any future experiment. Past calculationsof the local neutrino clustering factor only considered a spherical distribution of matter inthe Milky Way and neglected the influence of other nearby objects like the Virgo cluster,although recent N -body simulations suggest that the latter may actually be important. Inthis paper, we adopt a back-tracking technique, well established in the calculation of cosmicrays fluxes, to perform the first three-dimensional calculation of the number density of relicneutrinos at the Solar System, taking into account not only the matter composition of theMilky Way, but also the contribution of the Andromeda galaxy and the Virgo cluster. Theeffect of Virgo is indeed found to be relevant and to depend non-trivially on the value of theneutrino mass. Our results show that the local neutrino density is enhanced by 0.53% for aneutrino mass of 10 meV, 12% for 50 meV, 50% for 100 meV or 500% for 300 meV. a r X i v : . [ a s t r o - ph . C O ] O c t ontents N -one-body method 44 Computing neutrino clustering with back-tracking 55 Density profiles and gravitational potential 6 A.1 Spherical symmetry: the dark matter halo 18A.2 Spherical symmetry: the bulge 18A.3 Cylindrical symmetry: exponential profile 19
The existence of a background of relic neutrinos, produced in the early Universe in a similarway as the photons that constitute the Cosmic Microwave Background (CMB), is one ofthe yet unconfirmed predictions of the standard model of cosmology. While we have goodindirect indications, such as the number of relativistic species in the early Universe, whichis constrained [1] to be very close to the expected theoretical value N eff = 3 .
045 [2–4], orthe imprint of relativistic species on the CMB spectrum, which is compatible with those offree-streaming relics (see e.g. [5]), a direct probe of the existence of a background of relicneutrinos would be a major discovery and a confirmation of what we know about cosmologyand neutrinos. In particular, it would discard the possibility that ordinary neutrinos havedecayed at some stage during the evolution of the universe (see e.g. [6, 7]) or have beenproduced with unexpectedly low abundance (e.g. in low reheating scenarios [8]), while anotherform of dark radiation would contribute to N eff (cid:39) β -decaying nuclei [9], in particular tritium [10].Using this method one would look for a (small) peak in the electron energy spectrum oftritium due to the capture of relic neutrinos, just above the endpoint of β decay. Althoughsuch determination is a real challenge due to the required energy resolution (that must becomparable with the absolute value of the neutrino mass) and the high number of backgroundevents, coming from β decay, that the experiment has to distinguish from the signal, a projectnamed PTOLEMY [11] is nowadays starting to test innovative technology that could lead,for favorable values of the neutrino masses, to the first direct observation of relic neutrinos– 1 –12]. The successful observation of relic neutrinos by PTOLEMY would also offer the uniqueopportunity to study for the first time the interactions of non-relativistic neutrinos .Since the number of events that a direct detection experiment could measure dependslinearly on the local number density of relic neutrinos, it is important to have a preciseknowledge of how many relic neutrinos are present today at Earth. The average numberdensity of neutrinos predicted by the standard cosmological model is 56 cm − per family andper degree of freedom. Relic neutrinos, however, have a very small mean energy today (of theorder of 10 − eV), therefore their average number density can be enhanced because of thegravitational attraction of the matter content of the Galaxy, as well as other neighbouringgalaxies and galaxy clusters, provided that their masses are large enough to let them cluster atsmall scales. The calculation of the clustered number density of relic neutrinos was proposedfor the first time in [16] using a method based on the collisionless Boltzmann equation, and in[17] using a method called N -one-body simulations. The latter case consists in computing thetrajectories of several ( N ) independent test particles (one-body) in the evolving gravitationalpotential of the Galaxy, starting from some high redshift until today, and then reconstructingthe profile of the neutrino halo according to the final positions of all the test particles. Thesame method has been adopted later in [18, 19], where an updated treatment of the darkmatter (DM) and baryonic content of the Milky Way was considered.In this paper, we improve the calculation presented in [18]. For the first time, we takeinto account not only the contribution of the Milky Way, but also the contribution of relevant,nearby astrophysical objects, such as the Andromeda galaxy and the Virgo cluster. In orderto perform this task, we need to relax the assumption of spherical symmetry that has beenused in previous works. In sections 2 to 4 we discuss the theoretical aspects of the calculationand the practical implementations in our code. The treatment of the matter content of thetwo galaxies and the Virgo cluster is presented in section 5. In the two final sections 6 and7 we present and discuss our results on the local number density of relic neutrinos. The motion of the N test particles must be computed in an evolving matter background. Ongalaxy scales and at recent times, the behavior of at least the two heaviest neutrino states iswell captured by Newton’s theory, in which the motion equations can be obtained from thefollowing Lagrangian: L = a (cid:18) m ν v − m ν Φ( (cid:126)x, t ) (cid:19) , (2.1)where a is the scale factor, m ν is the mass of the test neutrino, v its velocity and Φ the grav-itational potential. Let us now write the corresponding Hamiltonian, expressed in Cartesiancoordinates x, y, z and the corresponding conjugate momenta p x , p y , p z : H = 12 am ν (cid:0) p x + p y + p z (cid:1) + am ν Φ( (cid:126)x, t ) . (2.2)We can now compute the equations of motion. Denoting with a dot the derivative withrespect to conformal time, we find p i = am ν ˙ x i , ˙ p i = − am ν ∂ Φ( (cid:126)x, t ) ∂x i , with x i = x, y, z . (2.3) Given the values of the mass splittings provided by neutrino oscillation experiments, the second-to-lightestneutrino mass eigenstate must be heavier than at least 8 meV, see e.g. [13–15], while the mean energy of relicneutrinos is of the order of 10 − eV. – 2 –or a spherically symmetric potential, the equations of motion simplify significantly due tothe conservation of angular momentum and are best expressed in spherical coordinates. Notethat previous works always considered a spherically symmetric Milky Way.The most efficient way to do the calculations, however, does actually not involve thesolution of the above equations. It is more convenient to rescale the momenta of the testparticles in order to eliminate the neutrino mass from the equations: u i = p i /m ν , (2.4)thus replacing p i → u i and m ν → u i using a single neutrino mass will allow to obtain the results for different neutrino masses,simply rescaling the parameter space volume appropriately (see [17, 19]).In order to solve the equations of motion, we need the gravitational potential Φ of theGalaxy as well as those of other nearby objects like Andromeda and Virgo. We will makeuse of the Poisson equation to obtain the contribution to the total gravitational potential ofeach component described by its energy density ρ ,4 πGa ρ = ∇ Φ( (cid:126)x, t ) , (2.5)where the Laplacian operator is in comoving coordinates. Since the Poisson equation is linear,it is always possible to solve it separately for the different constituents of the total matterdensity. When assuming spherical symmetry, the potential depends only on the distancefrom the center of the halo, r , so that it is possible to have an analytic expression for thederivative of the potential from equation (2.5): ∂ Φ ∂r ( r, z ) = 4 πGa r (cid:90) r ρ halo ( x, z ) x d x = Gar M halo ( r, z ) . (2.6)When the density is instead not symmetric, the Poisson equation must be solved numerically.The most convenient way is to use Fourier transforms, as we discuss in the appendix A. OnceΦ is obtained, one has to compute the partial derivatives that enter in the Hamilton equationsabove and that we discuss in details in section 5.The N -one-body simulation method requires the solution of the equations of motionof many test particles with different initial conditions. When dealing with the sphericallysymmetric case, one has to sample different values for the parameter space of only threequantities: the initial distance from the center of the halo, the initial momentum of theparticle, and the initial angle between the initial momentum vector and the radial direction.Because of spherical symmetry, the motion of each test particle will always be contained ina plane. The calculation of the number density profile of the relic neutrino halo (and of itsparticular value at Earth) takes into account the final position of all test particles weightedby their initial phase space, see [17]. If one wants to relax the assumption of sphericalsymmetry, however, the final position will have to be computed as a function of six inputvariables (three for the position and three for the momentum) instead of just three. Since thenumber of test particles that are required in order to obtain a sufficiently precise result scalesexponentially with the number of dimensions, repeating the calculation without sphericalsymmetry would require an unreasonably high number of computational hours. Moreover,many of the simulated test particles will end up very far from the position of the Earth, andwill give very little or no contribution to the local density of relic neutrinos.– 3 –ortunately, a simple way to overcome this problem has been known for many yearsin the context of cosmic ray propagation. Instead of forward-tracking the particles startingfrom homogeneous and isotropic initial conditions at high redshift, it is more efficient toconsider only those particles that are at Earth today. This is done by inverting the arrow oftime in the equations, and back-tracking the particles from our position today. Afterwards,we can attribute an initial phase-space volume and an appropriate statistical weight to eachtrajectory. The main advantage of this method is that one only needs to sample over the3-momentum of the neutrinos reaching the Earth today, since their position is fixed byassumption, regardless of the assumed symmetries of the astrophysical environment. Thecomputational time will thus remain comparable to that of previous works assuming sphericalsymmetry, while allowing us to introduce a complex distribution of matter with many objects.The drawback, however, is that with the back-tracking method one does not obtainthe full shape of the neutrino halo around the Earth, but only the local number density.To estimate the shape of the neutrino profile, multiple simulations at different positions arerequired. More details on the back-tracking method and on our specific implementation arediscussed in the next sections. N -one-body method Some of us used the forward-tracking technique in a previous publication [18], followingAppendix A.3 in [17] and using the kernel method of [20]. In this approach, the numberdensity is reconstructed from a set of N representative particles of the phase-space interval( r a , p r,a , p T,a ) i → ( r b , p r,b , p T,b ) i . Each trajectory is given a weight w i ( i = 1 , · · · , N ) w i = (cid:90) ( r b ,p r,b ,p T,b ) i ( r a ,p r,a ,p T,a ) i (cid:90) θ,φ,ϕ f ( p ) d r d p, (3.1)where f ( p ) is assumed to be the homogeneous and isotropic Fermi-Dirac distribution (neglect-ing small linear perturbations far from the Milky Way), while we used d r = r sin θ d θ d φ d r and d p = p T d p T d p r d ϕ , p T being the transverse momentum and p r the radial momentum.The final number density at radius r is then given by n ( r ) = N (cid:88) i =1 w i ξ K ( r, r i , ξ ) , (3.2)where the Gaussian kernel K ( r, r i , ξ ) = ξ (2 π ) / rr i exp (cid:18) − r + r i ξ (cid:19) sinh (cid:18) rr i ξ (cid:19) (3.3)plays the role of smearing the particles around the surface of a sphere in order to get a profilethat is spherically symmetric. The parameter ξ is the window width [20] and its value canbe optimized for each step in the simulation.When switching to the back-tracking method, we take the opposite perspective. Wedraw trajectories from samples of the neutrino momentum today, at the position of theEarth. At that time and location, the phase-space distribution of neutrinos is no longer closeto the Fermi-Dirac distribution of the average neutrino background, due to the non-lineardynamics inside the halo. Fortunately, we can make use of Liouville’s theorem to compute– 4 –he statistical weight of each phase-space volume element around the Earth to that at theother end of the trajectory, where neutrinos still obey the average homogeneous and isotropicFermi-Dirac distribution.Liouville’s theorem [21] implies the conservation of phase-space density along the solu-tions of the equations of motions,˙ f + ˙ (cid:126)x · (cid:126) ∇ (cid:126)x f + ˙ (cid:126)p · (cid:126) ∇ (cid:126)p f = 0 . (3.4)After tracking a particle from redshift z = 0, the Earth position (cid:126)x ⊕ and some arbitrarymomentum (cid:126)p j (0) back to redshift z back , position (cid:126)x j ( z back ) and momentum (cid:126)p j ( z back ), we cancompute the phase-space distribution today and in the right direction by applying f ( (cid:126)x ⊕ , (cid:126)p j (0) ,
0) = f back ( (cid:126)x j ( z back ) , (cid:126)p j ( z back ) , z back ) , (3.5)where f back can be identified with the Fermi-Dirac distribution of the average neutrino back-ground. This gives us f today in any direction. The final number density is then obtained byintegrating over the observed momentum (cid:126)p j (0), without any need for Gaussian smoothing.In this work we have chosen z back = 4, but we verified that the value of z back has no signifi-cant impact on the final number density, see section 6. Note that this Liouville mapping isroutinely used when back-tracking cosmic rays, see e.g. Ref. [22].In both methods, after obtaining the local number density n ν i ( (cid:126)x ⊕ ) for each mass eigen-state i , one can compute the clustering factor: f i ≡ n ν i /n ν, , (3.6)where n ν, = 112 cm − is the cosmological average number density for one family of neutrinosplus anti-neutrinos. For solving the equations of motion, we have used a symplectic ODE solver that also conservesphase-space volume, the symplectic_rkn_sb3a_mclachlan solver, that is the symmetricB3A method of the Runge-Kutta-Nystr¨om scheme of sixth order [23] from the odeint packageof the
Boost libraries [24].The symplectic solvers of the odeint package require the equations of motion to be sep-arable, that is the time-derivatives of the coordinates are functions of the conjugate momentaonly and vice versa , and autonomous, that is all right-hand sides must not depend on time t explicitly. The latter requirement poses a problem since both the expanding space-timeand the redshift evolution of the gravitational potential introduce a time dependence in theHamiltonian, cf. eq. (2.2). A common fix consists of treating time as an extra variable to beintegrated on top of u i ( t ) and x i ( t ), with a trivial derivative ˙ t = 1 (cf. e.g. [25]). With suchan addition, the system is formally autonomous and still separable. Finally, we note that ifwe substitute time for the new variable s , s ( z ) = − (cid:90) z d z ˙ a , (4.1) – 5 –he equations of motion (2.3) take on the even simpler formd x i d s = u i , d u i d s = − a ∂φ∂x i , (4.2)allowing to further speed up the computation.Although back-tracking dramatically reduces the number of particles to be simulated,we still have to deal with a large number of trajectories, obtained by solving eqs. (4.2) forseveral initial conditions u i ( z = 0). This requirement is most efficiently fulfilled in a “SingleInstruction, Multiple Data” (SIMD) architecture, modern graphic processing units (GPUs)being an example. We have used the CUDA framework via the Thrust library which canbe interfaced with odeint ’s solvers. In order to increase speed, we have pre-computed thebaryonic contributions to the gravitational potential (see next section) at redshift z = 0on a grid in cylindrical coordinates R and z , loaded them as textures onto the GPU, bi-linearly interpolated them between grid points, and finally scaled the results up to higherredshifts z . For the results below, we have isotropically sampled the arrival directions ofneutrinos (20 points for polar angle, 20 points for azimuth) and logarithmically sampled inmomentum (100 points over 3 decades), which leads to a grid of 4 × velocities. We havechecked that this is sufficient for getting well-converged clustering factors even in the non-axisymmetric case (Milky Way dark matter plus baryons plus Andromeda and Virgo). Allthe computations were performed on an Nvidia Quadro P6000. Depending on the numberof different contributions to the gravitational potential, back-tracking the 4 × particlesfrom redshift z = 0 to z = 4 required between 120 and 500 minutes. In this Section we describe how we implement the gravitational potential of the objects thatwe include in our analysis. For the Milky Way, we consider a spherical dark matter haloplus a number of baryonic components, which are described in axial symmetry. Beyond theMilky Way, we consider spherical dark matter halos for the Andromeda galaxy and the VirgoCluster, which are the largest objects relatively close to Earth that can have an impact onthe local density of relic neutrinos. Finally, we report technical details on the discretizationof the grid that we adopt in the numerical calculation for the interpolation of the derivatives.
For the dark matter in our Galaxy, we consider two distinct cases: a Navarro-Frenk-White(NFW, [26]) and an Einasto [27] profile, which read, respectively, ρ NFW ( r ) = ρ rR s (cid:18) rR s (cid:19) for r < R vir , (5.1) ρ Ein ( r ) = ρ exp[ − ( r/R s ) α ] , (5.2)where ρ is the normalization, R s is the scale radius, R vir is the virial radius of the NFWprofile (which is related to R s through the concentration parameter c = R vir /R s ), and α is http://developer.nvidia.com/cuda-zone http://thrust.github.io – 6 –n additional parameter for the Einasto profile that controls the change in the slope of thedensity.Concerning baryons, we follow the treatment of [28] and adopt five separate compo-nents: stars, warm and cold dust, atomic HI and molecular H gas. The density of stars isparameterized using a disk plus a bulge. The bulge of the Milky Way has been shown to havea triaxial shape [29]. However, since we are mainly interested in the neutrino clustering atthe Earth position, which is located at distances (8 . ± . n = 4 (i.e. a de Vaucouleurs profile): ρ DeVac ( r ) = ρ exp (cid:34) − A (cid:18) rR b (cid:19) / (cid:35) (cid:18) rR b (cid:19) − / , (5.3)where A = 2 n − / ≈ .
67 and we take R b = 0 .
74 kpc. The other baryonic components areassumed to be distributed according to a double exponential disk profile ρ exp ( R, z ) = ρ e − R/R s e −| z | /z s . (5.4)The present day values of the parameters of the different profiles are obtained as follows.The parameters related to the dark matter component of the Milky Way at z = 0 are obtainedby fitting the dark matter contribution to the rotation curve data as reported in [31], followingthe same procedure already adopted in [18] . Better estimates of the Galactic rotation curveare nowadays accessible thanks to the second data release of the ESA/Gaia mission [32] (seee.g. [33] for an analysis of the dark matter contribution to the rotation curve data presentedin [34]). However, given the limited radial extent of the data, instead of fixing the total darkmatter mass to the values predicted either in [31] or [33], we use an estimate based on orbitingMilky Way satellites up to ∼
300 kpc from the center of the Galaxy [35]. Concerning thebaryon components, we take the warm dust, cold dust, H and HI profile parameters from[28], as well as the scale parameters of the bulge and the disk. For the central value of thedensity of the bulge we use [30], which also provides an estimate for the total stellar mass inthe Galaxy (5 × M (cid:12) ). From this number we can derive the total mass of the stellar diskby subtracting the total mass of the bulge. The values of the parameters that we adopt arelisted in Tables 1 and 2.Regarding the HI density profile, observations (e.g. [36, 37]) have shown that the distri-bution of neutral hydrogen in the outskirts of the Galaxy follows an exponential profile, aswe assume in this work; conversely, the central 2.75 kpc [37] seem to be devoid of it. Thisfeature would in principle spoil the analyticity of our potentials. However, we found thatneglecting the central hole in the hydrogen distribution, i.e. extrapolating the exponentialprofile until the origin of the coordinates, would just cause an increase of the total HI massof 1%, which is in turn an overestimate of order 0.01% on the total mass of the Milky Way.We then feel safe to ignore such a feature in the HI profile, and we consider it to be also adouble exponential disk, following eq. (5.4).In order to compute the clustering factor today, we also need the time evolution ofthe density profiles. As we comment in section 6, most of the clustering happens at smallredshifts, so there is no need to compute the density profiles very precisely at all times. Notice that to switch from our parameterization of the Einasto profile in Eq. (5.2) to the one used by [18],one has to substitute ρ → ρ exp ( − /α ) and R s → R s (2 /α ) /α . – 7 – FW Einasto M vir [ M (cid:12) ] 2 . × . × ρ [ M (cid:12) /kpc ] 1 . × . × R s [kpc] 19.9 0.737 R vir [kpc] 333.5 (cid:55) α (cid:55) Table 1 . Dark matter density parameters for the Milky Way at z = 0, obtained by fitting the datafrom [31], following the same procedure as in [18]. ρ [ M (cid:12) /kpc ] R s [kpc] z s [kpc] M tot [ M (cid:12) ]Bulge 1 . × (cid:55) . × Disk 3 . × . × Warm dust 1 . × . × Cold dust 2 . × . × H . × . × HI 7 . × . × Table 2 . Density profile parameters for the baryonic components at z = 0. We also provide thetotal mass for each component. All the components have a profile described by equation (5.4), exceptfor the bulge, which follows a de Vaucouleurs profile (equation (5.3)). The scale radii and heightsare taken from [28], as specified in the main text. The redshift evolution of the total mass is foundfollowing the N -body simulation results of [38], while we assume that R s and z s do not evolve in time. The evolution in redshift of the density profiles is accounted for in the following way. Weassume the total virial mass of the dark matter halo to be constant in redshift, while theconcentration parameter changes according to [39]log βc vir ( z ) = a ( z ) + b ( z ) log (cid:18) M vir h − M (cid:12) (cid:19) , (5.5)where β is a parameter (considered constant in time) which denotes the offset of the MilkyWay concentration with respect to the average one. The functions a ( z ) and b ( z ) are differentfor the NFW and Einasto profiles. For the NFW they correspond to a ( z ) = 0 .
537 + (1 . − . (cid:2) − . z . (cid:3) , b ( z ) = − .
097 + 0 . z , while for the Einasto profile a ( z ) =0 .
459 + (0 . − . (cid:2) − . z . (cid:3) , b ( z ) = − .
130 + 0 . z .The evolution of the virial radius is obtained from M vir = 4 πa (cid:90) R vir ( z )0 ρ ( x, z ) x d x, (5.6) R vir ( z ) = (cid:18) M vir π ∆ vir ( z ) ρ crit ( z ) (cid:19) / , (5.7)where ρ crit = 3 H / (8 πG ) is the critical density of the Universe and ∆ vir ( z ) = 18 π +82 [Ω m ( z ) − −
39 [Ω m ( z ) − [40] for the NFW. For the Einasto profile it is instead fixed to∆ vir = 200, since this was the approach followed by [39] in order to obtain the numerical val-ues of the corresponding a ( z ) and b ( z ) equations. Combining equations (5.5), (5.6) and (5.7)allows us to find the scale radius as a function of redshift. The cosmology used in this work– 8 –as h = 0 . m = 0 . N -body simulations obtained by [38]. Inparticular we assume that the fraction of each component with respect to the total baryonmass is conserved.For the equations of motion, we need to obtain the derivatives of the gravitationalpotentials. A detailed description of the method we employ to compute the potentials andtheir derivatives for all the matter components of the Galaxy can be found in Appendix A.The derivative of the total Milky Way potential, split in all its matter components, is givenby ∂ Φ tot ∂x i ( x ) = ∂ Φ DM ∂x i ( r, ρ DM , R DM , R vir , DM ) dark matter (eq. 5.1/5.2)+ ∂ Φ b ∂x i ( r, ρ b , R b ) stellar bulge (eq. 5.3)+ ∂ Φ d ∂x i ( R, z, ρ d , R d , z d ) stellar disk (eq. 5.4)+ ∂ Φ w ∂x i ( R, z, ρ w , R w , z w ) warm dust (eq. 5.4)+ ∂ Φ c ∂x i ( R, z, ρ c , R c , z c ) cold dust (eq. 5.4)+ ∂ Φ H ∂x i ( R, z, ρ H , R H , z H ) H (eq. 5.4)+ ∂ Φ HI ∂x i ( R, z, ρ HI , R HI .z HI ) HI (eq. 5.4). (5.8) We also wish to incorporate in our system nearby objects whose presence may have a signif-icant impact on the clustering factor of neutrinos in the Milky Way. Results from N -bodysimulations in [41] (see their Figure 2) show that the neutrino halo of Virgo-like clusters mayextend up to distances comparable to the one between the Milky Way and the Virgo clusteritself. The neutrino overdensity caused by the Virgo halo at the Milky Way distance is ex-pected to be of a few percent, even for the minimum masses allowed by neutrino oscillations(Σ m ν = 60 meV). At the location of the Earth, we thus expect the Virgo effect to be almostof the same order of magnitude as the Milky Way effect.We assume that the dark matter halo of the Virgo cluster follows a NFW profile with amass of 6 . × M (cid:12) [42]. Its distance and position in the sky in Galactic coordinates aretaken from the NASA Extragalactic Database : D Virgo ≈ . ± .
0) Mpclatitude
Virgo = 74 . ◦ longitude Virgo = 283 . ◦ = ⇒ x Virgo = 1 .
056 Mpc y Virgo = − .
299 Mpc z Virgo = 15 .
895 Mpc . (5.9)We also include in our work the Andromeda galaxy, which is much lighter than Virgo (bya factor ∼ ∼
20) to the Milky Way. Also for Andromeda https://ned.ipac.caltech.edu/ – 9 – irgo cluster Andromeda M vir [ M (cid:12) ] 6 . × . × ρ [ M (cid:12) /kpc ] 8 . × . × R s [kpc] 399.1 21.8 R vir [kpc] 2328.8 244.7 Table 3 . Dark matter density parameters for the Andromeda galaxy and the Virgo cluster at z = 0.The parameters for Virgo are taken from [42], and for Andromeda from [43]. we consider a NFW profile, but we neglect its baryon content. The galactic latitude andlongitude of Andromeda are taken from the Vizier database , while its distance, mass anddensity profile parameters are given by [43], leading to D And ≈ . ± .
120 Mpclatitude
And = − . ◦ longitude And = 121 . ◦ = ⇒ x And = − .
377 Mpc y And = 0 .
623 Mpc z And = − .
288 Mpc . (5.10)The density parameters at z = 0 are listed in Table 3 for both Andromeda and theVirgo cluster. The evolution of the density profile parameters for these objects is governedby the same equations as for the Milky Way halo (see Section 5.1).The complete astrophysical setup we consider, with the Milky Way, Andromeda and theVirgo cluster, is shown in Figure 1. The size of the dots corresponds to the virial radius of theNFW halos. We can note that Virgo is much bigger and distant than Andromeda. However,as N -body simulations show [41], the neutrino halo of each object is always much moreextended than the one of dark matter, due to the high neutrino thermal velocities. Thus,despite its high distance from the Milky Way, we expect that Virgo will also contribute tothe neutrino overdensity at the Earth location. Solving the Hamiltonian equations of motion requires the derivative of the gravitationalpotentials listed in eq. (5.8). It is convenient to provide these derivatives explicitly to thecode in order to benefit from the use of textures in the GPU calculations.First of all, we safely assume that outside the virial radius of each dark matter halo, thepotential is just given by Kepler’s formula. In this way we do not have to build very broadgrids.Inside the halos, the choice of the grid size depends on how much we want to characterizethe halo itself. For us, the most interesting structure is of course the Milky Way. We wantour grid to be much finer than the distance between the Earth and the Galactic center ( ≈ ∼
50 kpc more than the maximum virial radius at z = 4. On the other hand,despite the fact that the Virgo cluster is much more extended than the Milky Way (its virial http://vizier.u-strasbg.fr/viz-bin/VizieR-S?NGC%20224 – 10 – [ M p c ] − − y [ M p c ] − − z [ M p c ] Figure 1 . Relative position of the Milky Way, Andromeda Galaxy and the Virgo Cluster. The sizeof the dots matches the virial radius of the object. The grey shaded plane represents the plane of theMilky Way. radius reaches up to 3 Mpc), we do not need a very narrow binning there, since we are notinterested in what happens on very small scales. We use a 1 kpc bin size in radius.After computing these derivatives in spherical coordinates as a function of the radius,we get the derivatives in Cartesian coordinates by means of the chain rule (see Appendix A).The baryonic components only have cylindrical symmetry, leading to a more subtlesituation. Their 2-D grid in R and z must extend at least up to a point where we can safelyapproximate the potential generated by a disk-like profile with the one generated by a sphereof the same mass. This depends of course on the ratio of scale radius and scale height: thelarger the ratio, the further the grid needs to extend before we approach a Keplerian law.Looking at Table 2, we see that in the Milky Way the maximum ratio between the scaleradius and the height of the disk is 50 (for cold dust). For this configuration, we compute thepotential of an exponential profile as well as its Keplerian counterpart (i.e. a point-like objectwith the same mass) to check where the two potentials start to be equal to each other. InFigure 2, we plot the different iso-potential lines for these two configurations: it is easy to seethat at distances of R ≈ R s , the red and black isocontours, which refer to the cylindricaland spherical case respectively, differ approximately by just 1%. We therefore extend thegrid on which we calculate the derivative of the potential to at least 30 times the largest scaleradius among all the components. All in all, for the Milky Way, we compute the derivativeof the potential up to ≈
550 kpc from the Galactic center.– 11 –
R/R s z / R s - - - - - - - - - - - - CylindricalSpherical, same mass − − − − [ ( k m / s ) ] Figure 2 . The colormap shows the potential generated by an exponential disk. The red lines denoteisocontours for this potential, while the black ones denote the isocontours for the potential generatedby a point-like source with the same mass. At
R/R s ∼ z/R s ∼
25 the difference between the sphericaland cylindrical potentials is smaller than 1%.
Milky Way Andromeda Virgo r / R . −
550 kpc 0 . −
350 kpc 1 − r / ∆ R z − −
550 kpc∆ log ( z ) 0.0337 Table 4 . Characteristics of the grid used for the derivative of the contributions to the potential.
The bin sizes must be chosen carefully, especially along the direction z orthogonal tothe baryonic disks. For dark matter, a bin size of 0.1 kpc would be sufficient, but some of thebaryonic components have a disk much thinner than that. We therefore opt for a logarithmicgrid in z that spans from 10 − to 550 kpc.All the above choices are summarized in Table 4. In Figure 3, we show the clustering factor, i.e. n ν i /n ν, , at the Earth’s position for a givenneutrino mass eigenstate as a function of m ν i , both for the case with an NFW distributionand an Einasto distribution for the dark matter in the Milky Way. For Virgo and Andromeda,– 12 – m ν [meV]1 + 10 − − n ν / n ν , NFW
NFWNFW + baryonsNFW + baryons + VirgoNFW + baryons + Virgo + AndromedaNFW (de Salas et al. )NFW + baryons (de Salas et al. )NFWhalo (Ringwald & Wong)MWnow (Ringwald & Wong)NFW + baryons (Zhang & Zhang) m ν [meV]1 + 10 − − n ν / n ν , Einasto
EinastoEinasto + baryonsEinasto + baryons + VirgoEinasto + baryons + Virgo + AndromedaEinasto (de Salas et al. )Einasto + baryons (de Salas et al. ) Figure 3 . For each neutrino mass state, we plot the ratio n ν /n ν, at the Earth’s position as a functionof the neutrino mass m ν . We consider contributions to the gravitational potential from the Galacticdark matter halo ( top panel: NFW profile, bottom panel:
Einasto profile), from baryons in the Galaxy,from the Virgo cluster and from the Andromeda galaxy. We also compare with earlier studies [17–19]. – 13 –e only consider dark matter with an NFW profile. We also compare our results with thoseof previous studies [17–19].As expected, regardless of our assumptions on the gravitational potential, the clusteringfactor increases with the neutrino mass. The impact of baryons in our Galaxy is significantfor any value of the mass. In contrast, adding the Virgo contribution leads to an enhance-ment at small neutrino mass, but can actually lead to less clustering at masses larger thanapproximately 200 meV. In the forward-tracking picture, this is easily explained as some ofthe neutrinos that would have clustered at the Earth’s position in the absence of Virgo arenow clustering in the Virgo potential well instead. In the back-tracking picture, a fraction ofthe particles sent out from the Earth that would have lost energy by leaving the Milky Way’sgravitational potential have fallen into Virgo’s gravitational potential instead. This leads toan increase in momentum of these particles with increasing redshift. Thus the phase-spacedensity is sampled only at large momenta for these particles (instead of all momenta), andthe clustering is overall less pronounced.We can also see in Figure 3 that both the effect of Andromeda and the difference betweenan NFW and an Einasto profile for the Milky Way’s dark matter are negligible. Assuming m ν = 50 meV, the overdensity is ( n ν /n ν, − (cid:39) z back . In the back-tracking approach, z back controls the time at which we stop integratingthe neutrino trajectories – while with forward-tracking it would be the initial redshift. Inboth cases, z back gives the time at which we assume a perfect homogeneous and isotropicFermi-Dirac distribution for the neutrinos. Figure 4 shows the reconstructed value of theclustering factor today when z back is floated – rather than being fixed to our baseline case ofaveraging z back ∈ [3 . , z back is floated in therange from 0 to 0.5, and a gradual convergence towards an asymptotic value for z back > z < .
5. Second, this convergence test proves that it is sufficientto assume a perfect homogeneous and isotropic Fermi-Dirac distribution for the neutrinosat z back . Indeed, in principle, one should either push the simulation up to z back → ∞ ,or introduce some small phase-space density fluctuations δf ( t back , (cid:126)x, (cid:126)p ) accounting for theamount of clustering that took place between the onset of structure formation and z back . Ifgravitational potential wells at z back were so large that such fluctuations should be taken intoaccount, neglecting them would introduce a bias in the results, that would depend on z back .A non-observation of this dependence shows that the clustering between z → ∞ and z back can be safely neglected.As one can see from Figure 4, for masses below 100 meV, the convergence of the clus-– 14 –ering factor is achieved for z back >
2. Instead, when the neutrino mass grows, we note thatthe solution is slightly less converged, due to the existence of trapped orbits for some of theneutrinos around the Milky Way and Virgo halos, which originate well before z = 4. Inthese cases, the value of z back can have an impact on the results, but the magnitude of theoscillations seen in Figure 4 shows that this is at most a 10% effect for n ν /n ν, −
1. Since thiseffect is smaller than the uncertainties coming from the assumptions on the dark matter andbaryon composition of the Galaxy, and that neutrino masses above 100 meV are disfavoredby cosmological measurements, we simply present the results (Figure 3) at high masses as anaverage of the values n ν /n ν, obtained considering z back ∈ [3 . , The cosmic neutrino background, one of the predictions of the standard cosmological model,has never been directly detected. Due to their extremely small energy, this is remarkablychallenging, although there are ongoing efforts to build the first relic neutrino detector. Inparticular, the PTOLEMY project [11], under development at Gran Sasso National Labora-tory in Italy, will use the capture of relic neutrinos by tritium [9, 10] as the process meantto unveil the existence of such elusive particles [12]. The event rate for the neutrino cap-ture process depends on the density of relic neutrinos at the location of the Earth, whichis expected to be larger than the average cosmological density ( n (cid:39)
56 cm − ) due to thegravitational attraction of our Galaxy.Past works studied the clustering of relic neutrinos in the Milky Way using a methodcalled “ N -one-body simulations” [17–19], which consists in computing the trajectories of N independent test particles, sampled assuming homogeneous and isotropic initial conditions athigh redshift, and in using the obtained information to reconstruct the shape of the neutrinohalo of the Galaxy. This method was found to be tractable only when assuming a sphericallysymmetric distribution of matter in our Galaxy and around: relaxing this assumption wouldrequire an enormous number of trajectory calculations to properly sample the six-dimensionalinitial phase space. In the context of the propagation of cosmic rays in the magnetic field of theGalaxy, a similar problem has been addressed since many years with a back-tracking method:instead of evolving the trajectories from all the possible initial positions and momenta, onecan compute them backwards from the location of the Earth. In that way, one only considersthe trajectories that are relevant for the calculation of the local cosmic rays flux.In this paper, we used the back-tracking technique to expand the reach of the original N -one-body method beyond the spherically symmetric case. We included in our calculationa more realistic (cylindrical) description of the baryonic components of our Galaxy, as well asthe contribution of two close-by astrophysical objects: the Andromeda galaxy and the Virgocluster. We found that the main contribution comes from the dark matter halo of the MilkyWay, especially for the largest considered neutrino masses. However, the Virgo cluster mustbe taken into account in order to obtain the correct number density for the smallest neutrinomasses. The effect of the Virgo cluster is not trivial, as its presence may actually divert someof the neutrinos that would otherwise cluster on the Milky Way if their mass (velocity) waslarge (small) enough. On the other hand, we find that the nearest galaxy with a reasonable In the forward picture, it is easier to understand the phenomenon: since neutrinos are already clusteringaround the Milky Way and the Virgo cluster at z = 4, their momentum distribution function is not thehomogeneous and isotropic Fermi-Dirac at such redshifts. In the backward case, one has to think that theneutrinos cannot escape the Milky Way and the Virgo cluster until higher redshifts. – 15 – . . . . . . . . . z − − n ν / n ν ,
10 meV50 meV100 meV300 meV NFWNFW + baryonsNFW + baryons + Virgo
Figure 4 . Clustering factor as a function of the earliest redshift z back at which neutrino trajectoriesare integrated, for different values of the neutrino mass and different astrophysical configurations. size, Andromeda, is giving an almost negligible contribution to the overall clustering. Forthis reason, other nearby galaxies can be safely ignored.We conclude quoting our results for a few representative values of neutrino masses.The two cases m ν = 10 meV and 50 meV are particularly interesting, because they standfor plausible values of the mass of the second and third neutrino mass eigenstates in theminimal normal hierarchy scenario (that is, when the lightest neutrino is massless and themass ordering is normal). Additionally, in the minimal inverted hierarchy scenario (when the– 16 –ightest neutrino is massless and the mass ordering is inverted), the two heaviest neutrinoshave a mass close to m ν = 50 meV. We also quote our results for a mass of 300 meV, in tensionwith recent cosmological bounds, but still well below the strong and model-independent limitcurrently set by KATRIN [44].For such masses of 10 meV, 50 meV, 100 meV and 300 meV, we obtain that the localnumber density of the relic neutrinos is respectively enhanced by 0.53%, 12%, 50% and 500%.with respect to the cosmological average. We therefore find that the local number density ofrelic neutrinos is 56 . − , 63 . − , 85 cm − and 300 cm − for these cases. Acknowledgments
SG and SP were supported by the Spanish grants SEV-2014-0398 and FPA2017-85216-P(AEI/FEDER, UE), PROMETEO/2018/165 (Generalitat Valenciana) and the Red Con-solider MultiDark FPA2017-90566-REDC, as well as the European Union’s Horizon 2020research and innovation program under the Marie Sk(cid:32)lodowska-Curie individual Grant Agree-ment No. 796941 (SG). They, together with GP, thank the Institute for Theoretical ParticlePhysics and Cosmology (TTK) at RWTH Aachen University for hospitality and supportduring the development of this work. PFdS acknowledge support by the Vetenskapsr˚adet(Swedish Research Council) through contract No. 638-2013-8993 and the Oskar Klein Centrefor Cosmoparticle Physics. GP was supported by the INFN INDARK PD51 grant.
A Solving the Poisson equation
In order to correctly determine the clustering factor of neutrinos in the Earth neighborhood,an accurate modelling of the gravitational potential, not only of the Milky Way but also ofthe surrounding structures, is required. The potential is related to the matter density viathe well-known Poisson equation: ∇ Φ( x ) = 4 πGρ ( x ) , (A.1)where x is a physical (i.e. not comoving) coordinate.In case of spherical symmetry, equation (A.1) is easy to solve:Φ( r ) = − πG (cid:20) r (cid:90) r d x x ρ ( x ) + (cid:90) ∞ r d x x ρ ( x ) (cid:21) . (A.2)This solution can be applied to all the components of the Galaxy which satisfy sphericalsymmetry, namely the dark matter halo and the bulge.The other baryonic components of our Galaxy which we wish to incorporate, such as gasand stars, are distributed in a disk, so that a more general solution to the Poisson equationmust be used. In terms of Green’s functions such solution writesΦ( x ) = − G (cid:90) d r ρ ( r ) | x − r | , (A.3)but the integral becomes more and more difficult to solve with the increasing complexity ofthe system. For this reason we employ a different approach. When Fourier transforming bothsides of equation (A.1), the Laplacian operator converts into a factor −| k | , independently– 17 –f the geometry of the system. Moving this factor at right-hand side and switching back toconfiguration space givesΦ( x ) = − πG (cid:90) d k (2 π ) e i k · x | k | (cid:20)(cid:90) d r ρ ( r ) e − i k · r (cid:21) . (A.4)This is the equation we are going to use for the baryon components in the Milky Way. A.1 Spherical symmetry: the dark matter halo
For dark matter halos, we consider two distinct cases of density profiles, namely the NFWand Einasto profiles, written in equations (5.1) and (5.2), respectively. We solve the Poissonequation by simply working out the integrals in equations (A.2):Φ
NFW ( r ) = − πGρ R s ln (cid:16) mR s (cid:17) r/R s − R vir /M R vir R s (A.5)Φ Ein ( r ) = − πGρ R s α ×× (cid:20) Y − /α Γ (cid:18) α , , Y (cid:19) + Γ (cid:18) α , Y, ∞ (cid:19)(cid:21) , (A.6)where m = min( r, R vir ), M = max( r, R vir ), Y = ( r/R s ) α and we have defined the incompleteΓ function: Γ( a, L, U ) = (cid:90) UL d t t a − e − t . (A.7)To compute the derivative with respect to any axis x i we use the chain rule: ∂ Φ ∂x i = dΦd r ∂r∂x i = dΦd r x i r . (A.8) A.2 Spherical symmetry: the bulge
As we mentioned in the paper, the bulge of the Milky Way is ellipsoidal rather than spherical,with a ratio of semi-axes of ∼ .
6. However, for our purposes we can safely neglect thedifference between the semi-axes, and consider the density profile given in equation (5.3).Again, the gravitational potential can be computed using equation (A.2):Φ
DeVac ( r ) = − Gπρ R s A / [ F in ( X ) + F out ( X )] , (A.9)– 18 –here X = A (cid:16) rR s (cid:17) / and we have defined F in ( X ) = 116 X (cid:40) √ π erf (cid:16) X / (cid:17) − e − X X / ×× (cid:16) X + 960 X + 6240 X + 34320 X ++ 154440 X + 540540 X + 1351350 X + 2027025 (cid:17)(cid:41) (A.10) F out ( X ) = (cid:40) √ π erfc (cid:16) X / (cid:17) ++ 2 e − X X / (cid:0) X + 28 X + 70 X + 105 (cid:1) (cid:41) . (A.11)Again, computing the derivatives of the potential (A.9) is quite straightforward usingthe chain rule: ∂ Φ DeVac ∂x i = dΦ DeVac d X d X d r ∂r∂x i = dΦ DeVac d X X r x i r . (A.12) A.3 Cylindrical symmetry: exponential profile
Not all the components in a halo have spherical symmetry: the dust, gas and stellar disk ofthe Milky Way have an axial symmetry with respect to the Galactic plane. We can thereforeexploit cylindrical coordinates. We assume that all the matter component which satisfy axialsymmetry display an exponential profile, given by equation (5.4).In this particular case the solution of the Poisson equation given by eq. (A.4) is (almost)analytical: Φ exp ( R, z ) = − πGρ z s R s ×× (cid:90) ∞ d k k k e − k | z | − z s e −| z | /z s (cid:104) kR s ) (cid:105) / [1 − k z s ] J ( kR ) . (A.13)The integral appearing in equation (A.13) is a Hankel transform of order 0: we solve itusing the FFTlog code .The derivatives of this potential are also analytical. For the x and y directions, weexploit the chain rule combined with the property of Bessel functions d J α ( x )d x = − J α +1 ( x ),while the derivative with respect to the z -axis is obtained by just deriving the integrand Specifically we use the python package
FFTLog ( https://github.com/prisae/fftlog ) by D. Werthm¨uller,based on the Fortran FFTLog code by A. Hamilton [45, 46]. – 19 – − R [kpc]10 − − − ρ = 7 . × M (cid:12) / kpc R s = 5 . z s = 0 . z/z s = 0 . → . − Φ( R, z ) [(km / s) ] − F R ( R, z ) [(km / s) kpc − ] − Φ Kepler − F R, Kepler
Figure 5 . The colored solid lines represent the potential generated by an exponential disk with thesame parameters as the cold dust profile in the Milky Way. Dashed lines represent the derivative withrespect to R of this potential, i.e. the force per unit mass in the radial direction. The different colorlabel different z : the blue curves are the closest to the Galactic plane, the red ones are the farthestones. For comparison, also the Keplerian potential and force are plotted (solid and dashed black linesrespectively): for radii much larger than R s (dotted vertical line) both the potential and its derivativeconverge to this value. function: ∂ Φ exp ∂ ( x, y ) = 4 πGρ z s R s ( x, y ) R ×× (cid:90) ∞ d k k k (cid:0) k e − k | z | − z s e −| z | /z s (cid:1)(cid:104) kR s ) (cid:105) / [1 − k z s ] J ( kR ) , (A.14) ∂ Φ exp ∂z = 4 πGρ z s R s z | z | ×× (cid:90) ∞ d k k e − k | z | − e −| z | /z s (cid:104) kR s ) (cid:105) / [1 − k z s ] J ( kR ) . (A.15)The behavior of these last equations can be seen in Figure 5, where the example of thecold dust disk in the Milky Way is shown. The solid lines represent the gravitational potentialas a function of radius for different values of z , indicating with blue (red) the closest (farthest)to the Galactic plane. Dashed lines in turn are the derivative of this potential with respect to– 20 – . As can be noticed, both the potential and its derivative converge to the Keplerian limit(solid and dashed black) when going sufficiently far outside the scale radius (dotted verticalline). References [1]
Planck
Collaboration, N. Aghanim et al.,
Planck 2018 results. VI. Cosmological parameters , arXiv:1807.06209 .[2] G. Mangano, G. Miele, S. Pastor, T. Pinto, O. Pisanti, and P. D. Serpico, Relic neutrinodecoupling including flavor oscillations , Nucl. Phys. B (2005) 221–234, [ hep-ph/0506164 ].[3] P. F. de Salas and S. Pastor,
Relic neutrino decoupling with flavour oscillations revisited , JCAP (2016) 051, [ arXiv:1606.06986 ].[4] S. Gariazzo, P. F. de Salas, and S. Pastor, Thermalisation of sterile neutrinos in the earlyUniverse in the 3+1 scheme with full mixing matrix , JCAP (2019) 014,[ arXiv:1905.11290 ].[5] B. Audren et al.,
Robustness of cosmic neutrino background detection in the cosmic microwavebackground , JCAP (2015) 036, [ arXiv:1412.5948 ].[6] J. F. Beacom, N. F. Bell, and S. Dodelson, Neutrinoless universe , Phys. Rev. Lett. (2004)121302, [ astro-ph/0404585 ].[7] Z. Chacko, A. Dev, P. Du, V. Poulin, and Y. Tsai, Cosmological Limits on the Neutrino Massand Lifetime , arXiv:1909.05275 .[8] P. F. de Salas, M. Lattanzi, G. Mangano, G. Miele, S. Pastor, and O. Pisanti, Bounds on verylow reheating scenarios after Planck , Phys. Rev. D (2015) 123534, [ arXiv:1511.00672 ].[9] S. Weinberg, Universal Neutrino Degeneracy , Phys. Rev. (1962) 1457–1473.[10] A. G. Cocco, G. Mangano, and M. Messina,
Probing low energy neutrino backgrounds withneutrino capture on beta decaying nuclei , JCAP (2007) 015, [ hep-ph/0703075 ].[11] PTOLEMY
Collaboration, E. Baracchini et al.,
PTOLEMY: A Proposal for Thermal RelicDetection of Massive Neutrinos and Directional Detection of MeV Dark Matter , arXiv:1808.01892 .[12] PTOLEMY
Collaboration, M. G. Betti et al.,
Neutrino Physics with the PTOLEMY project , JCAP (2019) 047, [ arXiv:1902.05508 ].[13] F. Capozzi, E. Lisi, A. Marrone, and A. Palazzo, Current unknowns in the three neutrinoframework , Prog. Part. Nucl. Phys. (2018) 48–72, [ arXiv:1804.09678 ].[14] P. F. de Salas, D. V. Forero, C. A. Ternes, M. T´ortola, and J. W. F. Valle,
Status of neutrinooscillations 2018: 3 σ hint for normal mass ordering and improved CP sensitivity , Phys. Lett.B (2018) 633–640, [ arXiv:1708.01186 ].[15] I. Esteban, M. C. Gonzalez-Garcia, A. Hernandez-Cabezudo, M. Maltoni, and T. Schwetz,
Global analysis of three-flavour neutrino oscillations: synergies and tensions in thedetermination of θ , δ CP , and the mass ordering , JHEP (2019) 106, [ arXiv:1811.05487 ].[16] S. Singh and C.-P. Ma, Neutrino clustering in cold dark matter halos : Implications forultrahigh-energy cosmic rays , Phys. Rev. D (2003) 023506, [ astro-ph/0208419 ].[17] A. Ringwald and Y. Y. Y. Wong, Gravitational clustering of relic neutrinos and implicationsfor their detection , JCAP (2004) 005, [ hep-ph/0408241 ].[18] P. F. de Salas, S. Gariazzo, J. Lesgourgues, and S. Pastor, Calculation of the local density ofrelic neutrinos , JCAP (2017) 034, [ arXiv:1706.09850 ]. – 21 –
19] J. Zhang and X. Zhang,
Gravitational Clustering of Cosmic Relic Neutrinos in the Milky Way , Nature Commun. (2018) 1833, [ arXiv:1712.01153 ].[20] D. Merritt and B. Tremblay, NonParametric Estimation of Density Profiles , Astrophys. J. (1994) 514.[21] H. Goldstein, C. Poole, and J. Safko,
Classical mechanics . Addison Wesley, 2002.[22] P. Mertsch,
Test particle simulations of cosmic rays , arXiv:1910.01172 .[23] R. I. McLachlan, On the Numerical Integration of Ordinary Differential Equations bySymmetric Composition Methods , SIAM J. Sci. Comput. (1995) 151–168.[24] B. Sch¨aling, The Boost C++ Libraries . XML Press, 2014.[25] S. Blanes and P. C. Moan,
Splitting Methods for Non-autonomous Hamiltonian Equations , J.Comp. Phys. (2001) 205–230.[26] J. F. Navarro, C. S. Frenk, and S. D. M. White,
The Structure of cold dark matter halos , Astrophys. J. (1996) 563–575.[27] J. Einasto,
On the Construction of a Composite Model for the Galaxy and on theDetermination of the System of Galactic Parameters , Trudy Astrofizicheskogo InstitutaAlma-Ata (1965) 87–100.[28] A. Misiriotis, E. M. Xilouris, J. Papamastorakis, P. Boumis, and C. D. Goudis, The distributionof the ISM in the Milky Way. A three-dimensional large-scale model , Astron. Astrophys. (2006) 113–123, [ astro-ph/0607638 ].[29] M. Portail, C. Wegg, O. Gerhard, and I. Martinez-Valpuesta,
Made-to-Measure models of theGalactic Box/Peanut bulge: stellar and total mass in the bulge region , Mon. Not. Roy. Astron.Soc. (2015) 713–731, [ arXiv:1502.00633 ].[30] J. Bland-Hawthorn and O. Gerhard,
The Galaxy in Context: Structural, Kinematic, andIntegrated Properties , Ann. Rev. Astron. Astrophys. (2016) 529–596, [ arXiv:1602.07702 ].[31] M. Pato and F. Iocco, The Dark Matter Profile of the Milky Way: a Non-parametricReconstruction , Astrophys. J. (2015) L3, [ arXiv:1504.03317 ].[32]
Gaia
Collaboration, A. G. A. Brown et al.,
Gaia Data Release 2 , Astron. Astrophys. (2018) A1, [ arXiv:1804.09365 ].[33] P. F. de Salas, K. Malhan, K. Freese, K. Hattori, and M. Valluri,
On the estimation of theLocal Dark Matter Density using the rotation curve of the Milky Way , JCAP (2019) 037,[ arXiv:1906.06133 ].[34] A.-C. Eilers, D. W. Hogg, H.-W. Rix, and M. Ness, The Circular Velocity Curve of the MilkyWay from to kpc , Astrophys. J. (2019) 120, [ arXiv:1810.09466 ].[35] L. L. Watkins, N. W. Evans, and J. H. An,
The masses of the Milky Way and Andromedagalaxies , Mon. Not. Roy. Astron. Soc. (2010) 264–278, [ arXiv:1002.4565 ].[36] P. M. W. Kalberla, W. B. Burton, D. Hartmann, E. M. Arnal, E. Bajaja, R. Morras, andW. G. L. P¨oppel,
The Leiden/Argentine/Bonn (LAB) Survey of Galactic HI. Final data releaseof the combined LDS and IAR surveys with improved stray-radiation corrections , Astron.Astrophys. (2005) 775–782, [ astro-ph/0504140 ].[37] P. J. McMillan,
The mass distribution and gravitational potential of the Milky Way , Mon. Not.Roy. Astron. Soc. (2017) 76–94, [ arXiv:1608.00971 ].[38] F. Marinacci, R. Pakmor, and V. Springel,
The formation of disc galaxies in high resolutionmoving-mesh cosmological simulations , Mon. Not. Roy. Astron. Soc. (2014) 1750–1775,[ arXiv:1305.5360 ]. – 22 –
39] A. A. Dutton and A. V. Macci`o,
Cold dark matter haloes in the Planck era: evolution ofstructural parameters for Einasto and NFW profiles , Mon. Not. Roy. Astron. Soc. (2014)3359–3374, [ arXiv:1402.7073 ].[40] G. L. Bryan and M. L. Norman,
Statistical Properties of X-Ray Clusters: Analytic andNumerical Comparisons , Astrophys. J. (1998) 80–99, [ astro-ph/9710107 ].[41] F. Villaescusa-Navarro, J. Miralda-Escud´e, C. Pe˜na-Garay, and V. Quilis,
Neutrino halos inclusters of galaxies and their weak lensing signature , JCAP (2011) 027,[ arXiv:1104.4770 ].[42] P. Fouque, J. M. Solanes, T. Sanchis, and C. Balkowski,
Structure, mass and distance of thevirgo cluster from a tolman-bondi model , Astron. Astrophys. (2001) 770,[ astro-ph/0106261 ].[43] P. R. Kafle, S. Sharma, G. F. Lewis, A. S. G. Robotham, and S. P. Driver,
The need for speed:escape velocity and dynamical mass measurements of the Andromeda galaxy , Mon. Not. Roy.Astron. Soc. (2018) 4043–4054, [ arXiv:1801.03949 ].[44]
KATRIN
Collaboration, M. Aker et al.,
An improved upper limit on the neutrino mass from adirect kinematic method by KATRIN , arXiv:1909.06048 .[45] J. D. Talman, Numerical Fourier and Bessel transforms in logarithmic variables , J. Comput.Phys. (1978) 35–48.[46] A. J. S. Hamilton, Uncorrelated modes of the nonlinear power spectrum , Mon. Not. Roy.Astron. Soc. (2000) 257–284.(2000) 257–284.