Local clustering of relic neutrinos with kinetic field theory

The density of relic neutrinos is expected to be enhanced due to clustering in our local neighbourhood at Earth. We introduce a novel analytical technique to calculate the neutrino overdensity, based on kinetic field theory. Kinetic field theory is a particle-based theory for cosmic structure formation and in this work we apply it for the first time to massive neutrinos. The gravitational interaction is expanded in a perturbation series and we take into account the first-order contribution to the local density of relic neutrinos. For neutrino masses that are consistent with cosmological neutrino mass bounds, our results are in excellent agreement with state-of-the-art calculations.


I. INTRODUCTION
While at large scales the formation of structures in our Universe can be described within the well-established framework of cosmological linear perturbation theory, on smaller scales our predictions almost entirely rely on Nbody simulations [1]. In general, N -body simulations reproduce observed structures in remarkable detail and allow to derive testable features. Nonetheless, such simulations, unlike analytical approaches, are very computationally demanding and may lack physical interpretability. The most desirable state would be an interplay between N -body simulations and (semi-)analytical methods. Unfortunately, analytical methods allow reliable predictions up to scales only slightly beyond the linear scale [2]. An interesting and relatively recent approach to structure formation is provided by kinetic field theory (KFT) [3][4][5][6][7][8][9]. The main effort of KFT in cosmology has so far been the calculation of the matter power spectrum 1 , showing promising results. The formalism, however, might seem relatively complicated and-being based on a path integral formulation of classical particles-is very different from standard analytical methods. As such, this recent development seems to have obtained relatively little attention. The idea of this work is to take one step back and test the reliability of KFT on a considerably less complicated problem, namely the calculation of the local overdensity of relic neutrinos. As our final results * ebholm@phys.au.dk † isabel.oldengott@uclouvain.be ‡ szentarra@ethz.ch 1 Beyond that, an application to baryonic accoustic oscillations can be found in [10] and velocity power spectra have been studied in [11].
look promising for realistic neutrino masses, this work can be seen as introducing a novel technique to approach this problem.
The relic or cosmic neutrino background (CνB) can be seen as the neutrino analogue of the cosmic microwave background (CMB). Its presence has so far been confirmed only indirectly, by measurements of the number of relativistic degrees of freedom N eff through the CMB anisotropy spectrum (e.g. [12]) and through light element abundances from big bang nucleosynthesis (BBN) (e.g. [13,14]). A direct detection would however still be the ultimate dream of every neutrino cosmologist. Originating back to temperatures T ∼ 1 MeV, it would provide the earliest direct probe of the Universe that we have. Due to the extremely low energies of the relic neutrinos today, a direct detection is unfortunately extremely challenging -if even feasible at all [15]. Nevertheless, the ambitious PTOLEMY project [16] has taken on the challenge of a first direct detection of the CνB. A crucial parameter for PTOLEMY (and for any other future attempt to measure the CνB) is the local neutrino density which directly impacts the expected event rate (by its impact on the Tritium β-decay rate induced by neutrino capture). While the homogeneous and isotropic background neutrino density is expected to be around 56 cm −3 per flavour (and as many for anti-neutrinos), the local neutrino density is in general expected to be enhanced because neutrinos start clustering when they become non-relativistic. The current strongest experimental constraint on the local neutrino overdensity from the KATRIN experiment is < 1.1 × 10 11 (95% CL) [15]. A comparable or somewhat stronger bound is expected from future experiments for cosmogenic neutrinos [17].
From the theoretical side, this problem has first been approached in [18] where it was pointed out that calculating the local neutrino clustering does not require arXiv:2305.13379v2 [hep-ph] 2 Aug 2023 full N -body simulations but can be simplified significantly by what we here dub the N -one-body approximation. Within this approximation, the relic neutrinos are treated as test particles in an external gravitational potential governed by cold dark matter (CDM) and baryons. In other words, neutrino clustering neither impacts the clustering of CDM, baryons nor itself. The validity of the N -one-body approximation was explicitly confirmed by [19] and is the basis for our work and all methods introduced in this section. Based on this assumption, [18] used the linearized Vlasov equation to solve for the local neutrino density. [20] then showed that the linearized Vlasov equation is only justified for sufficiently small neutrino overdensities by performing N -one-body simulations. Instead of simultaneously following the evolution of N particles (neutrinos and CDM) as in a conventional N -body simulation, the application of the N -one-body approximation allows the tracking of one particle at a time through N independent simulations, which reduces the computation time significantly. The same method-with various improvements on the modelling of the gravitational potentialhas been adopted in [21,22]. [20][21][22] all assumed a spherically symmetric gravitational potential to keep the computational effort in a reasonable scale. Most recently, [23] has extended the computation of the local overdensity to a non-spherical gravitational potential and applied the back-tracking technique. Here, the trajectories of the neutrinos ending up here and now are tracked back by reversing the time direction in the equations of motion. The local overdensity is then derived by assigning a statistical weight to each trajectory. Long story short, the final conclusion of [23] is that the local relic neutrino overdensity is [0.5%, 12%, 50%, 500%] for neutrino masses [10, 50, 100, 300] meV. For neutrino masses not in tension with cosmological neutrino mass bounds [12] the overdensity is hence expected to be < 50%.
While we thus already have a theoretical prediction for the local neutrino overdensity, the problem at hand offers the perfect conditions to test the reliability of KFT: Due the application of the N -one-body approximation, the KFT formalism reduces to a formalism of one-particle quantities. Therefore, we avoid two major complications arising in the context of calculations of the matter power spectrum with KFT, namely correlated initial conditions [4,24,25] and inter-particle (gravitational) interactions [4,8,9]. Furthermore, as just summarized, the problem has been well explored with other methods which offers the opportunity to directly compare our results with state-of-the art methods.

II. THE FORMALISM
The formalism of KFT for an ensemble of N particles has been developed thoroughly in [5] (plus references within) and [6][7][8][9]. It is intuitively clear (and can be shown rigorously) that the use of the N -one-body ap-proximation justifies the formulation of the problem in terms of one-particle quantities. While the formalism presented in this section looks significantly simpler than in most of the literature on KFT, we note that all equations can be derived from the formalism as introduced in [8].
In a nutshell, KFT is a path integral formulation of classical mechanics. Its central object is the generating functional, averaged over initial conditions, Here, the first line depicts the averaging over initial comoving positions ⃗ x (i) and initial conjugate momenta ⃗ p (i) . For the initial phase-space distribution we assume a relativistic, homogeneous and isotropic Fermi-Dirac distribution, with T 0 denoting the neutrino temperature today. We denote the norm of the momentum by p (i) . ⃗ x s is the solution to the Hamiltonian equations of motion and ⃗ J is the so-called source function 2 . Note that we explicitly choose redshift z as our time coordinate. Physical observables can be derived from the generating functional by means of functional differentiation with respect to the source function. As such, the number density in Fourier space can be derived according to In order to evaluate eq. (2) let us derive the equations of motion in the following. The Lagrangian with redshift z as a time coordinate is given by We assume a two-component Universe for the Hubble expansion rate H(z) = H 0 Ω m,0 (1 + z) 3 + Ω Λ , which is justified as neutrino clustering happens at sufficiently late times.
Here Ω m,0 and Ω Λ are the matter density and cosmological constant density parameters today, and H 0 is the Hubble rate today. The effective gravitational potential φ on the RHS of eq. (3) fulfills the comoving Poisson equation (see e.g. the appendix of [3]), where ρ is the matter density perturbation. The conjugate momentum to the comoving position is which leads to the Hamiltonian and hence the equations of motion If the equations of motion (7)-(8) had an analytical solution, we could now derive an exact analytical expression for the number density (2) in terms of an integral over initial phase-space. As this is not the case, we follow the perturbative treatment of the generating functional outlined in [8]. For that purpose, we write down the formal solution of eqs. (7)- (8), with the Green's function where 2 F 1 refers to the Gaussian hypergeometric function. Following the derivation in [8] (for the case of an N -particle ensemble), one can then write the generating functional in eq. (1) as ≡ exp iŜ I Z 0 ⃗ J .
Here, we identified eq. (14) (including the integration over initial conditions in eq. (12)) as the free generating functional Z 0 and eq. (13) as an interaction operator exp(iŜ I ) acting on it. Applying the series definition of the exponential exp(iŜ I ) leads to a perturbation series for the generating functional (1) and in turn to a perturbation series for the number density, In this work, we restrict our analysis to KFT in firstorder perturbation theory. As discussed in sec. 4.4 of [8], the particle trajectories from the first-order generating functional correspond to trajectories within the Born approximation. To be more explicit, those are trajectories in eq. (9) where the gradient of the gravitational potential is evaluated along the free trajectory. At zeroth order, applying the functional derivative in eq. (2) on the free generating functional directly gives the standard expression for the homogeneous and isotropic background number density, i.e.
The derivation in appendix A shows that the first order contribution to the local relic neutrino number density is The above expression has a transparent interpretation: At first-order, the contribution to the relic neutrino overdensity is obtained by integrating the external matter density along the free trajectory back in time from the position ⃗ x.

III. RESULTS
Eq. (18) and the numerical calculation thereof form the main result of this work. We compare it to the most state-of-the art technique to calculate the local neutrino overdensity, i.e., the back-tracking method applied by [23]. While [23] showed that in general non-spherical contributions to the gravitational potential (in particular the Virgo cluster and baryonic components of the Milky way) have a significant impact on the local neutrino overdensity, the aim of this work is to provide a test of KFT and a comparison between the two methods. For this purpose, the modelling of the non-spherical potential as described in [23] would present an unnecessary complication 3 . We therefore restrict our work to the case of a spherically symmetric NFW potential [26] representing the dark matter halo of the Milky Way and model it exactly as described in [23]. We truncate the NFW profile at its virial radius R vir (z), where ρ 0 (z) denotes the central density and R s (z) the scale radius. We refer the reader to [23] for the time evolution of ρ 0 (z), R s (z) and R vir (z) as well as all numerical values that are needed to fully characterize the density profile. As in [23], we assume the Earth to be located at a distance of 8.2 kpc from the galactic center [27] and assume the cosmological parameters Ω m,0 = 0.3111, Ω Λ = 0.6889 and H 0 = 67.77 km s −1 Mpc. Numerically, we evaluate the integral (18) as a fourdimensional iterated integral. However, by choosing a coordinate system in which ⃗ x = x⃗ e z , one can show that the magnitude of the vector ⃗ x−g(z, z ′′ )⃗ p (i) is independent of the azimuthal angle of the initial momentum ⃗ p (i) , reducing the dimensionality of the integral by one. We take the redshift integral as the outermost integral and evaluate the hypergeometric function, occurring in the Green's function (11), with a (12, 12) minimax rational approximation [28]. The redshift and polar integrals are evaluated using adaptive quadrature with a 15 point Gauss-Kronrod formula [29] and the dp (i) integral is performed using Gauss-Laguerre quadrature [28] since it is weighted by the initial Fermi-Dirac distribution. Without parallellization, the integral only takes on the order of 10 ms to evaluate 4 . This is many orders of magnitude faster than the back-tracking method, which takes about 120 minutes, according to [23]. However, our assumption of spherical symmetry reduces the dimensionality of the integral (18) compared to the computation of [23], which of course makes a one-to-one comparison of the computation times difficult. Figure 1 shows the estimates of the local relic neutrino overdensity n/n − 1 from eq. (16) and eq. (18) compared to the state-of-the-art results of [23] as a function of the neutrino mass m ν (the blue line corresponds to the line dubbed as "NFW" in fig. 3 of [23]). The lower panel shows the relative difference between the results of both methods. The two computations are seen to agree at the percent level for masses up to ≈ 0.1 eV, while the firstorder KFT result underestimates the neutrino clustering at larger masses. At m ν = 0.3 eV the difference between both methods is around 47%.
The fact that the result from first-order KFT underestimates the overdensity can be well understood by keeping in mind that it corresponds to Born trajectories that underestimate the deviation from straight lines (or free 4 Our code is publicly available at https://github.com/EBHolm/KFT-Neutrinos on the branch with the arXiv ID of this paper. trajectories). This can also be seen as an underestimation of the gravitational force and it therefore comes as no surprise that the neutrino clustering is underestimated as well. For low neutrino masses this effect is less severe because the true trajectories are less curved (due to the higher velocities) and the Born approximation performs better. Furthermore, for high neutrino masses one can in general expect neutrinos to follow bound trajectories, which KFT is unable to reproduce at first order.

IV. CONCLUSION
In this work, we have for the first time applied KFT to calculate the local overdensity of relic neutrinos. For neutrino masses consistent with cosmological mass bounds, first-order KFT shows excellent agreement with the state-of-the art calculation by [23] that applies the backtracking method. For larger neutrino masses (m ν ≳ 0.1 eV) the overdensity is underestimated (up to around 47% at m ν = 0.3 eV). The method introduced in this work is a priori not restricted to relic neutrinos but can be applied to other light particles, under the assumption that the N -one-body approximation holds.
We focused on a spherically symmetric NFW potential governed by dark matter. This is justified as our work should be rather seen as a proof-of-principle for the reliability of KFT-and not so much as a precision calculation of the local clustering of relic neutrinos which would require accounting for baryonic contributions from the Milky Way and the Virgo cluster (see [23]). Though the calculation of the local relic neutrino density is not a cosmological problem limited by computation time, it is remarkable that the numerical evaluation of eq. (18) only takes about 10 ms (per neutrino mass). To put this into context, the back-tracking method takes around 120 minutes according to [23]. However, caution has to be allowed with this direct comparison of computation times as the generalization to a non-symmetric gravitational potential like in [23] could increase the computation time in a non-trivial way. However, even in comparison to previous works based on N -one-body simulations [20][21][22], which also assumed spherical symmetry, the first-order KFT computation is considerably faster. On the other hand, it was pointed out early on in [20] that for relatively low overdensities (and hence for what we today consider as realistic neutrino masses) solving the linearized Vlasov equation already provides reliable results. Indeed, the contribution from first-order KFT in eq. (18) has a form similar to the linear response expression in [20] (or [30] in the context of N -body simulations including massive neutrinos), where neutrinos respond linearly to an external gravitational potential. Note however that in case of the linearized Vlasov equation the perturbation parameter is the phase-space density while in KFT it is the gravitational interaction. It is a priori not clear how both of these approaches compare. In a future work we will compare KFT to the linearized Vlasov equation (and its numerical performance). We believe that such a comparison is also relevant for a general comparison between KFT and standard cosmological perturbation theory [31]. Furthermore, we leave the investigation of higher-order contributions to the clustering for future work.
In conclusion, KFT reproduces the standard results of a homogeneous and isotropic background neutrino density and performs remarkably well at first order. For the specific calculation of the local clustering of relic neutrinos and for neutrino masses consistent with cosmological mass bounds, the method introduced in this work is competitive with the state-of-the art technique based on back-tracking [23].