Cosmic-Enu: An emulator for the non-linear neutrino power spectrum

Cosmology is poised to measure the neutrino mass sum $M_\nu$ and has identified several smaller-scale observables sensitive to neutrinos, necessitating accurate predictions of neutrino clustering over a wide range of length scales. The FlowsForTheMasses non-linear perturbation theory for the massive neutrino power spectrum, $\Delta^2_\nu(k)$, agrees with its companion N-body simulation at the $10\%-15\%$ level for $k \leq 1~h/$Mpc. Building upon the Mira-Titan IV emulator for the cold matter, we use FlowsForTheMasses to construct an emulator for $\Delta^2_\nu(k)$ covering a large range of cosmological parameters and neutrino fractions $\Omega_{\nu,0} h^2 \leq 0.01$, which corresponds to $M_\nu \leq 0.93$~eV. Consistent with FlowsForTheMasses at the $3.5\%$ level, it returns a power spectrum in milliseconds. Ranking the neutrinos by initial momenta, we also emulate the power spectra of momentum deciles, providing information about their perturbed distribution function. Comparing a $M_\nu=0.15$~eV model to a wide range of N-body simulation methods, we find agreement to $3\%$ for $k \leq 3 k_\mathrm{FS} = 0.17~h/$Mpc and to $19\%$ for $k \leq 0.4~h/$Mpc. We find that the enhancement factor, the ratio of $\Delta^2_\nu(k)$ to its linear-response equivalent, is most strongly correlated with $\Omega_{\nu,0} h^2$, and also with the clustering amplitude $\sigma_8$. Furthermore, non-linearities enhance the free-streaming-limit scaling $\partial \log(\Delta^2_\nu / \Delta^2_{\rm m}) / \partial \log(M_\nu)$ beyond its linear value of 4, increasing the $M_\nu$-sensitivity of the small-scale neutrino density.


INTRODUCTION
Cosmology will measure the neutrino mass sum M ν = m ν , one of the final unmeasured parameters of the Standard Model of particle physics, over the next several years, assuming that the dark energy is a cosmological constant (Audren et al. 2013;Chudaykin & Ivanov 2019).Future space-based experiments will provide completely independent bounds on M ν (Petracca et al. 2016;Lin et al. 2022).However, both forecasts and analyses of current data are consistent with a weakening of the M ν bound by a factor of ≈ 3 when the dark energy equation of state is allowed to vary with time (Font-Ribera et al. 2014;Upadhye 2019;Di Valentino et al. 2020).Additionally, neutrinos and dark energy may play roles in the resolution of persistent tensions in measurements of the cosmic expansion and clustering amplitude (Leauthaud et al. 2017;Böhringer & Chon 2016;Poulin et al. 2018;Gogoi et al. 2021;Di Valentino & Melchiorri 2022;McCarthy et al. 2018McCarthy et al. , 2023)).
On the particle physics side, persistent anomalies in the neutrino sector motivate models containing additional "sterile" neutrinos as well as non-standard neutrino interactions (Denton 2022;Aguilar-Arevalo et al. 2022;Alvarez-Ruso & Saul-Sala 2021).Furthermore, other hot dark matter (HDM) species such as axions could mimic the cosmological effects of massive neutrinos (Giarè et al. 2022;D'Eramo et al. 2022;Di Valentino et al. 2023).A quantitative understanding of neutrino clustering in the non-linear regime will prove invaluable for breaking these degeneracies.
In cosmology, several neutrino clustering signatures deep in the non-linear regime have been identified and quantified through Nbody simulations.These include "wakes" of neutrinos streaming coherently past halos (Zhu et al. 2016;Inman et al. 2015); an odd-parity contribution to the angular momentum field of galaxies (Yu et al. 2019); neutrino-dark-matter relative velocities (Zhu et al. 2014;Inman et al. 2017;Zhu & Castorina 2020;Zhou et al. 2022); modifications to the halo mass function (Costanzi et al. 2013;Yu et al. 2017;Biswas et al. 2019;Bocquet et al. 2020;Ryu & Lee 2022); a neutrino contribution to the scale-dependent bias of dark matter halos (LoVerde & Zaldarriaga 2014;LoVerde 2014;Chiang et al. 2018Chiang et al. , 2019;;Banerjee et al. 2020); and inhomogeneities in the cosmic neutrino background detectable in laboratory searches (Baracchini et al. 2018;Betti et al. 2019;Aker et al. 2022).Since the systematic bias associated with these phenomena are substantially different from those of joint analyses of large-scale cosmic surveys, some of them may be decisive to a convincing detection of massive neutrinos.
Accurate theoretical modeling of non-linear neutrino clustering effects, possibly through an extension of the halo model to neutrino clustering, will require fast and reliable calculations of the neutrino power spectrum.A recent code comparison by the Euclid collaboration, Adamek et al. (2023), tested a wide range of simulation methods (specifically, Schaller et al. 2023a;Teyssier 2002;Mauland et al. 2023;Adamek et al. 2016bAdamek et al. ,a, 2017;;Beck et al. 2016;Marin-Gilabert et al. 2022;Springel 2005;Springel et al. 2008Springel et al. , 2021;;Dakin et al. 2019Dakin et al. , 2022) ) and found agreement at the 30% − 40% level over the range k ≤ 1 h/Mpc of wave numbers.However, simulations are computationally expensive, motivating an exploration of alternative methods with comparable accuracies.
Linear perturbative calculations of massive neutrino clustering can be carried out to high precision by tracking the evolution of the neutrino distribution function (Ma & Bertschinger 1995) in both position and momentum.An alternative approach introduced by Dupuy & Bernardeau (2014Bernardeau ( , 2015a,b) ,b) discretizes the Fermi-Dirac distribution of initial neutrino velocities.Each neutrino fluid, defined by its initial velocity, can then be treated as a separate fluid individually obeying the continuity and Euler equations of fluid dynamics.This approach is Lagrangian in momentum space, since neutrinos cannot move from one initial-velocity bin to another.Since non-linear cosmological perturbation theory begins with the continuity and Euler equations, this momentum-Lagrangian method is a natural starting point for non-linear neutrino perturbation theories.
FlowsForTheMasses, the first non-linear perturbative power spectrum calculation for free-streaming particles such as massive neutrinos, began with precisely this approach (Chen et al. 2023b).Since a fluid with non-zero initial velocity ⃗ v has a preferred direction v, Chen et al. (2023b) began by extending the Time-Renormalization Group perturbation theory of Pietroni (2008); Lesgourgues et al. (2009) to fluids with homogeneous initial velocities.Their Fourierspace clustering depends not only upon the magnitude of the Fourier vector ⃗ k, but also its angle with v, whose cosine is µ = k • v.
Expanding the density contrast and velocity divergence in Legendre polynomials in µ, Chen et al. (2023b) showed that the modecoupling integrals of non-linear perturbation theory couple different Legendre moments, drastically increasing the computational cost.However, by applying Fast Fourier Transform (FFT) techniques introduced by McEwen et al. (2016), Schmittfull et al. (2016), and Fang et al. (2017) to its mode-coupling integrals, Chen et al. (2023b) was able to accelerate them by more than two orders of magnitude.The resulting FlowsForTheMasses perturbation theory can compute a non-linear massive neutrino power spectrum with reasonable accuracy settings on a standard desktop computer. 1hough the computational cost of FlowsForTheMasses is much lower than that of N-body neutrino simulations, it remains somewhat high; the fifty-flow production runs of Chen et al. (2023b) each took about a day on a 32-core machine.A machine learning technique known as emulation, introduced into cosmology by Heitmann et al. (2010), Heitmann et al. (2009), and Lawrence et al. (2010), is ideal for quickly approximating expensive functions.Emulation begins with a training set of evaluations of an expensive function (in our case the FlowsForTheMasses neutrino power spectrum) the size of which is determined by computational budget and required level of accuracy.To mitigate the computational expense associated with emulating multivariate data, the spectra are represented via a principal component (PC) decomposition.Gaussian process (GP) models are then used to model the cosmology-dependent PC weights, enabling fast prediction at new cosmologies.
Our goal in this work is an emulator of the FlowsForTheMasses non-linear neutrino power spectrum.Furthermore, since FlowsForTheMasses already divides neutrinos by their initial momenta and tracks the evolution of each one separately, our training set includes individual neutrino momentum deciles at no extra computational cost.For both the z ≤ 2 CDM+baryon power spectra and the emulator design points, we use the Mira-Titan IV (MT4) emulator of Moran et al. (2023).We demonstrate that our emulator, Cosmic-Eν, precisely reproduces FlowsForTheMasses to < 3.5% for 10 −3 h/Mpc ≤ k ≤ 1 h/Mpc and 0 ≤ z ≤ 3. 2Next, we compare Cosmic-Eν to a range of N-body neutrino simulations.We begin with the Euclid code comparison of Adamek et al. (2023) for M ν = 0.15 eV, which runs simulations with a va-riety of resolutions, box sizes, and massive neutrino implementations.Compared with their highest-resolution SWIFT simulation of Schaller et al. (2023b), Cosmic-Eν is accurate to 3% up to k = 3k FS = 0.17 h/Mpc, 19% to k = 0.4 h/Mpc, and 49% to k = 1 h/Mpc, all of which are somewhat larger than but comparable to the scatter among very different simulation methods.This qualitative picture is unchanged when M ν is raised to 0.3 eV and 0.6 eV in the Adamek et al. (2023) simulations, and when the dark energy is allowed to vary rapidly in one of our own simulations (described below).
Finally, we employ Cosmic-Eν to study enhancements to the ratios ∆ 2 ν /∆ 2 ν [LR] and ∆ 2 ν /∆ 2 m due to the non-linear clustering of massive neutrinos.The first of these, the non-linear enhancement of neutrino clustering relative to linear response (LR), was considered in Chen et al. (2023b).After confirming the accuracy of Cosmic-Eν for this quantity, we quantify the sensitivity of ∆ 2 ν /∆ 2 ν [LR] to each of the eight cosmological parameters, showing that Ω ν,0 h 2 is by far the most significant for determining the neutrino clustering enhancement, followed by Ω b,0 h 2 , Ω m,0 h 2 , and σ 8 .The second ratio, ∆ 2 ν /∆ 2 m , was shown by Ringwald & Wong (2004) and Wong (2008) to scale as the fourth power of M ν , hence Ω ν,0 h 2 , for LR neutrinos in the free-streaming limit.We confirm this result in the linear case, then show that non-linear corrections enhance this scaling relation.

For example, at k
ν .This study is organized as follows.Section 2 briefly describes the emulation procedure and the FlowsForTheMasses perturbation theory.Our emulator training set is assembled in Sec. 3 after improving the high-M ν numerical stability of FlowsForTheMasses.Section 4 constructs the Cosmic-Eν emulator and quantifies its accuracy with respect to FlowsForTheMasses.Cosmic-Eν is then compared with a wide variety of N-body simulation methods in Sec. 5. Finally, Sec. 6 quantifies the non-linear enhancements to the ∆ 2 ν /∆ 2 ν [LR] and ∆ 2 ν /∆ 2 m ratios, and Sec.7 concludes.

Emulation
A thorough discussion of emulation in cosmology may be found in Heitmann et al. (2009).Here, we briefly summarize their procedure, with slight differences in notation.Suppose that we wish to approximate a dimensionless function P(k, z, ⃗ C) of the wave number k, redshift z, and cosmological parameters ⃗ C.This may be proportional to the power spectrum itself, or a function of the power spectrum chosen to reduce its dynamic range.At each of m ∈ [0, N M ) cosmological models defined by parameters ⃗ C * m , we are given We seek an approximation of the form where µ * i is the mean input P * im across cosmologies, the ϕ j are a set of N PC orthogonal basis functions to be defined below, and w j ( ⃗ C ) are the corresponding basis weights.If the number of bases N PC is chosen equal to the number of training models N M minus one3 , then the approximation of Eq. ( 1) can be made exact for the N M input models, though we will typically choose N PC smaller than this.
Let D * im = (P * im − µ * i )/σ * , so D * im for fixed i has zero mean by construction, and D * im is an N kz × N M matrix.By means of a compact singular value decomposition, we may write where U is an N kz × N M orthogonal matrix, D an N M × N M diagonal matrix, and V an N M × N M orthonormal matrix.In terms of these three matrices, we may write the PC basis functions ϕ j (k i , z i ) and the weights of the input data w * jm as where the j index may be truncated to j < N PC for any chosen N PC < N M .Since the functions ϕ j (k, z) may be interpolated from the ϕ i j using standard methods, our remaining task is to model the weight functions w j ( ⃗ C ) using w * jm .We model each w j ( ⃗ C ) using a GP.A GP is an infinite dimensional generalization of a multivariate Gaussian distribution, in which any finite set of random variables is defined to follow a multivariate Gaussian distribution specified by a mean function and a covariance function (Williams & Rasmussen 2006).We define the GP over each w j ( ⃗ C ) to have a mean function of 0 and a Gaussian correlation function R j ( ⃗ C, ⃗ C ′ , ⃗ β j ).This correlation function is specified by a set of correlation hyperparameters β jℓ , one for each of the N C cosmological parameters of ⃗ C: The input data weights w * jm are now assumed to arise from the following hyperparameter-dependent probability distribution: where δ (K) mn is the Kronecker delta.Heitmann et al. (2009) has included an additional set of hyperparameters, a scaling term λ U, j and a "nugget" term λ W, j .The former scales the correlation function into a covariance function, while the latter accommodates slight numerical fluctuations in the computation of P * im .For each principal component j, we thus have one λ U, j , one λ W, j , and N C different β jℓ , for a total of N C + 2 hyperparameters.
Since we are given w * jm but not the hyperparameter values, our next step is to find the hyperparameter values most consistent with w * jm .We do so using a 50, 000-step Markov chain Monte Carlo sampling through the SEPIA code of Gattiker et al. (2020) 4 , using the default hyperparameter priors and bounds defined in SEPIA.Let λU, j , λW, j , and ⃗ β j be the posterior mean values of these hyperparameters.
We have now arrived at our goal, a predictive model for the weights w j ( ⃗ C ) in Eq. ( 1) for a given cosmology ⃗ C. Let Θ = { ⃗ β j , λU, j , λW, j }.Using these optimal hyperparameter values, we specify the conditional Gaussian probability distribution of each of the weights as follows: Here, ⃗ w * j is the N M × 1 vector of observed weights, R * j is the N M × N M matrix of realizations of Eq. ( 9) evaluated at the optimal values of the hyperparameters, ⃗ r * j ( ⃗ C ) is an N M × 1 vector having mth entry 10) indicates that w j is drawn from a normalized Gaussian distribution of mean W j ( ⃗ C ) and standard deviation Σ j ( ⃗ C ).As in Heitmann et al. (2009), we use W j ( ⃗ C ) as our emulated prediction of w j ( ⃗ C ).We see from Eq. ( 11) that this mean value is a weighted average of the observed weights ⃗ w * j , with weights determined by the covariance between the observed and predictive cosmological parameters.
Thus far we have assumed the initial data P * im to be given.Choosing the N M input models used to train the emulator is a separate problem known as emulator design.The goal of emulator design is to cover the given parameter space to a specified accuracy using the smallest number N M of input cosmological models, since each model is computationally expensive.A simple grid in parameter space is one of the least efficient designs for this purpose.Heitmann et al. (2009Heitmann et al. ( , 2016)); Lawrence et al. (2017); Moran et al. (2023) employ efficient space-filling Latin hypercube designs or similar nested, space-filling lattices; the latter is useful when runs are to be done in "batches" so that interim analyses may be performed as partial results are available.The training sets we construct, described in Sec.3.1, are built on the design choices made by these authors.
The batched MT4 emulator design allows for the addition of more design points to improve its accuracy.We leave this possibility open for future work.However, we will see in Sec.4.2 that even with the existing MT4 design, Cosmic-Eν attains an accuracy of ≈ 3.5%, which is subdominant to the 14% error in FlowsForTheMasses itself, as measured by Chen et al. (2023b).Thus we do not pursue here the possibility of including additional design points.

Non-linear perturbation theory for neutrinos
We are particularly interested in the non-linear growth of the neutrino density perturbations, which occurs at late times and at scales well within the Hubble horizon.Thus, to excellent approximation, we assume that all matter in the Universe obeys the scalarized nonrelativistic continuity, Euler, and Poisson equations in a box expanding uniformly at a rate given by the time-dependent Hubble parameter.General Relativistic clustering including vector and tensor perturbations, as well as multiple fluids, has been considered previously (e.g., Hwang & Noh 2006b,a, 2007, 2013b,a;Hwang et al. 2016;Jeong et al. 2011;Gong et al. 2017;Yoo & Zaldarriaga 2014;Yoo 2014;Magi & Yoo 2022;Adamek et al. 2014Adamek et al. , 2016aAdamek et al. , 2017;;Fidler et al. 2015Fidler et al. , 2016Fidler et al. , 2017Fidler et al. , 2019)).Generalizing these results to multiple fluids with different non-zero initial velocities is beyond the scope of the present study, as well as unnecessary to our goal of accuracy in the non-linear regime.Furthermore, following (Moran et al. 2023), we restrict our consideration to spatially-flat cosmologies.
The chief difficulty in applying the continuity and Euler fluid equations to neutrinos is that, in the Eulerian fluid description, neutrinos are not fluids.At each spatial point, neutrinos have a velocity dispersion arising from their initial Fermi-Dirac distribution; the number of degrees of freedom required to describe neutrino perturbations is therefore technically infinite.Führer & Wong (2015) devised a Eulerian fluid-like perturbation theory for the neutrino bispectrum by absorbing the infinite degrees of freedom into temporally non-local couplings, while Garny & Taule (2021, 2022) combined a linear treatment of free streaming with non-linear corrections to the density and velocity monopole perturbations.
On the other hand, the approach of Dupuy & Bernardeau (2014, 2015a,b) and Chen et al. (2021aChen et al. ( , 2023b)), which we describe in below, is to formulate a neutrino perturbation theory that is Lagrangian in momentum space.Let ⃗ τ be the lower-index threemomentum in the limit of an unperturbed universe, P (0)  i , which is time-independent.The Fermi-Dirac distribution may be binned by the magnitude τ = |⃗ τ| of this momentum, allowing us to approximate the neutrino population using N τ momenta, τ α , for α ∈ [0, N τ ).
For a given momentum vector ⃗ τ α in bin α, the set of neutrinos with initial momentum ⃗ τ α has no thermal velocity dispersion.Thus it behaves as a fluid obeying the continuity and Euler equations.Spatial isotropy implies that these equations depend upon the direction of ⃗ τ α only through its angle with respect to the Fourier vector ⃗ k, that is, through µ α := ⃗ k •⃗ τ α /(kτ α ).Chen et al. (2021a) demonstrates that this µ α -dependence can be expanded in N µ Legendre polynomials P ℓ (µ α ), using an appropriate boundary term at ℓ = N µ − 1, with an error of order N −2 µ .Furthermore, for given τ α , all neutrino fluids with initial momenta ⃗ τ such that |⃗ τ| = τ α obey the same fluid equations.We use the term "flow" for this entire set of fluids.
The subhorizon non-relativistic linear theory of Chen et al. (2021a) began with dimensionless scalar perturbations to the density, δ α (⃗ x) = (ρ α (x)− ρα )/ ρα , and the momentum divergence, θ α (⃗ x) = − ⃗ ∇ • ⃗ P/(m ν aH), with m ν the neutrino mass and H the conformal Hubble expansion rate.Since linear theory allows us to choose an arbitrary normalization for the perturbation variables, we normalize them to the square roots of their corresponding power spectra, where P ℓ (x) is the Legendre polynomial of order ℓ.Above we have employed two conventions we will use henceforth.
(i) Power spectrum indices b and c in P αbc take the value 0 for δ and 1 for θ, so, for example, P α00 = P αδδ .
(ii) Wave number superscripts denote a functional dependence, so that P ⃗ k αbc = P αbc ( ⃗ k) and δ k αℓ = δ αℓ (k).The perfect correlation between the random variables corresponding to δ α ( ⃗ k) and θ α ( ⃗ k) in linear theory breaks down beyond the linear order, leading the FlowsForTheMasses perturbation theory of Chen et al. (2023b) to introduce the quantities Since P α01 = P α10 , Eqs. (14-17) completely specify the power spec- αbc in terms of the perturbation variables.
In terms of the bispectrum integrals I k α,acd,be f,ℓ of Chen et al. (2023b), the evolution of δ αℓ , θ αℓ , and χ αℓ is given by where v α := τ α /(m ν a) is the flow velocity, and primes denote derivatives with respect to η := log(a/a in ) for given initial scale factor a in .
The gravitational potential Φ is given by the Poisson equation where δ cb is the density contrast of the CDM and baryons, treated as a single fluid, labelled cb; Ω cb (η) = Ω cb,0 H 2 0 /(aH 2 ) is the timedependent density fraction of this cb fluid; Ω cb,0 its density fraction today; ) the time-dependent density fraction of neutrino flow α; and Ω α,0 its density fraction today.
The bispectrum integrals I k α,acd,be f,ℓ are defined and thoroughly studied in Chen et al. (2023b).For our purposes, we may define them by their evolution equations with initial conditions I k α,acd,be f,ℓ = 2A k α,acd,be f,ℓ at η = 0. Here, the mode-coupling integrals are given by (25) =: with all other γ abc vanishing.In Eqs.(22,25), we have assumed implicit summation over the indices g and h for compactness.Numerical computation of these mode-coupling integrals A k α,acd,be f,ℓ is the main computational expense of FlowsForTheMasses perturbation theory.Chen et al. (2023b) shows that they may be reduced to Fast Fourier Transforms (FFTs) and then computed using the methods of Hamilton ( 2000 3 TRAINING DATA SET

Emulator design
We begin by constructing the training set P * im upon which the emulator is built.Emulator design is described broadly in Heitmann et al. (2009), and the particular design of the Mira-Titan IV (MT4) emulator upon which Cosmic-Eν is built is described in Heitmann et al. (2016), Lawrence et al. (2017), andMoran et al. (2023).The 111 design points in cosmological parameter space, ⃗ C * m , are chosen to strike a balance between broad parameter coverage and a high density of points, necessary for achieving high accuracy.Out of these, N M = 101 design points have non-zero neutrino masses, with physical density fractions Ω ν,0 h 2 ranging from 0.00017 to 0.01, corresponding to M ν from 0.0158 eV to 0.931 eV.
Our strategy is to build upon the MT4 emulator design.Since massless neutrinos are already described well by linear theory, Cosmic-Eν uses only these N M = 101 massive-neutrino points.We build our emulator upon MT4 for two reasons.Firstly, the MT4 design has already been optimized and thoroughly tested for a balance between breadth and accuracy.Secondly, the MT4 CDM+baryon power spectrum is used as an input to the FlowsForTheMasses neutrino perturbation theory.Since MT4 is most accurate at its own design points, choosing this same design for Cosmic-Eν avoids the compounding of errors that would result from using the outputs of one emulator as the inputs for another.Another potential error is the backreaction of enhanced neutrino clustering on the CDM+baryon power.Chen et al. (2021a) quantified the linear response backreaction to be 0.05% for Ω ν,0 h 2 = 0.01, the largest value considered here, while Section 6.3 argues that non-linear clustering only increases this by a factor of ∼ 3. Thus this backreaction is a negligible source of error for Cosmic-Eν.
Table 1 lists the parameter ranges covered by the Cosmic-Eν emulator.The allowed range of Ω ν,0 h 2 values is determined by the set of 101 massive-neutrino MT4 models.Its lower bound of 0.00017 is over three times smaller than the lower bound imposed by laboratory oscillation experiments (de Salas et al. 2018;Capozzi et al. 2018;Esteban et al. 2020).In Sec.4.2 we will quantify Cosmic-Eν errors near this low-Ω ν,0 h 2 boundary.For still smaller values, the nonrelativistic-neutrino approximation made by FlowsForTheMasses becomes increasingly inaccurate, and we instead recommend the use of relativistic linear perturbation theories such as CLASS (Lesgourgues 2011; Blas et al. 2011;Lesgourgues & Tram 2011) and CAMB (Lewis et al. 2000;Lewis & Bridle 2002).The upper bound on Ω ν,0 h 2 , consistent with M ν = 0.931 eV, allows for a broad exploration of the parameter space, which may be useful for finding solutions to the Hubble and σ 8 tensions (McCarthy et al. 2018;Di Valentino & Melchiorri 2022;McCarthy et al. 2023).As with the MT4 emulator, we have assumed three degenerate-mass neutrinos.
Allowed ranges on the remaining seven parameters are taken directly from the MT4 emulator.We parameterize the dark energy equation of state as w(z) = w 0 + w a z/(1 + z), following Chevallier & Polarski (2001); Linder (2003).The final range in Table 1 implies that if w 0 = −1, then −1.77 ≤ w a ≤ 0.99.Allowing significant variation in the dark energy equation of state is important for cosmological constraints, as Upadhye (2019) showed that this variation weakens the neutrino mass bound by a factor of ≈ 3.

Stabilizing perturbation theory
Our next challenge is the high-k numerical instability of the FlowsForTheMasses perturbation theory.The mode-coupling integral A k α,acd,be f,ℓ at low k and large ℓ rises sharply with k, increasing its dynamic range.Since FFTs spread errors across the entire range, small numerical errors near the peak of A α,acd,be f,ℓ lead to large fractional errors where A α,acd,be f,ℓ is small.These lead to instabilities, preventing integration of the equations of motion at high k.Additionally, the computational cost of the full set of mode-coupling integrals A k α,acd,be f,ℓ of Eqs.(25-26) was shown by Chen et al. (2023b) to scale as N 6 µ , further motivating a truncation in the range of ℓ.In practice, FlowsForTheMasses integrates 128 values of the wave number, logarithmically distributed between 10 −4 h/Mpc and 10 h/Mpc, and sets a stability threshold k st to 10 h/Mpc at the beginning of integration.The integration step size in η is chosen dynamically.Each time that high-k numerical instabilities drive this step size below 10 −6 , the integrator discards the highest k value, effectively lowering k st by 9.1%, and then resumes integration for k ≤ k st .Since non-linear physics causes power to flow from low to high k, we may safely discard wave numbers above this stability threshold without affecting those below it.We will see below that no significant noise or discontinuities affect the power spectrum in the range k ≤ k st .As an added precaution, we will require stability up to a threshold k st that is 20% larger than our largest wave number of interest, 1 h/Mpc.
In order to stabilize FlowsForTheMasses for its fiducial model, with Ω ν,0 h 2 = 0.005, Chen et al. (2023b) truncated the power spectra P k αbcℓ used to compute A ⃗ k α,acd,be f to ℓ < N µ,NL , while allowing N µ,NL ≤ ℓ < N µ elsewhere.They demonstrated that N µ,NL of 6, 7, and 8, respectively agreed with their N-body neutrino simulations to 14%, 12%, and 10% for k ≤ 1 h/Mpc.However, N µ,NL of 9 suffered from severe numerical instabilities preventing the equations of motion from being integrated to the present time.Henceforth we fix N µ,NL = 6, since its substantial reduction in computational expense relative to 7 and 8 only modestly decreases its accuracy.
Although the ℓ < N µ,NL truncation of Chen et al. (2023b) sufficed to stabilize their Ω ν,0 h 2 = 0.005 model over the range k ≤ 3 h/Mpc, we find that increasing the neutrino density fraction tends to exacerbate the numerical instabilities in FlowsForTheMasses.The MT4 emulator of Moran et al. (2023) allows Ω ν,0 h 2 to be twice as high as the fiducial model of Chen et al. (2023b).Requiring stability up to k = 1.2 h/Mpc, so as to allow for a buffer around our desired range k ≤ 1 h/Mpc, we find that 19 of the 101 design models are numerically unstable (that is, have k st < 1.2 h/Mpc).The mean and minimum Ω ν,0 h 2 values for these models are 0.0083 and 0.0063, respectively, so this instability is a high-Ω ν,0 h 2 problem.
As Ω ν,0 h 2 , hence the neutrino mass sum, increases, flow velocities v α decrease, gradually decoupling the neutrino density and velocity monopoles from higher Legendre moments ℓ.Thus we employ a more aggressive high-ℓ truncation in order to stabilize these high-Ω ν,0 h 2 models.In addition to truncating the power spectrum used Figure 1.Accuracy of truncation ℓ < N µ,AI in mode-coupling and bispectrum integrals.(Top) Ω ν,0 h 2 = 0.005.All N µ,AI are stable up to k = 2 h/Mpc, so the power spectrum for each has been divided by that for N µ,AI = 11, the maximum value.(Bottom) MT4 design model with Ω ν,0 h 2 = 0.0091.Each power spectrum has been divided by that for N µ,AI = 5, the largest value for which the calculation is stable up to k = 1.2 h/Mpc.
to compute mode-coupling integrals, we truncate the ℓ expansions of the mode-coupling integrals A k α,acd,be f,ℓ and bispectrum integrals I k α,acd,be f,ℓ themselves, ℓ < N µ,AI .The power truncation ℓ < N µ,NL itself implies an N µ,AI of 2N µ,NL − 1, or 11 for N µ,NL = 6, so we allow N µ,AI to be reduced below this number.
As an upper bound on the error associated with this truncation, consider the Ω ν,0 h 2 = 0.005 fiducial model of Chen et al. (2023b).Figure 1 (Top) compares neutrino power spectra for several N µ,AI values to the maximum value 2N µ,NL − 1 = 11.In the range k ≤ 1 h/Mpc, N µ,AI = 4 is accurate to 4%, while the higher N µ,AI considered are accurate to 1.3% or better.
Figure 1 (Bottom) makes a similar comparison for one of the MT4 design points with Ω ν,0 h 2 = 0.0091.Since power spectra with N µ,AI > 5 cannot be stably integrated to z = 0 all the way to k = 1.2 h/Mpc, we use the N µ,AI = 5 power spectrum for comparison.Encouragingly, the N µ,AI = 4 power spectrum agrees with this to better than 1% at all k ≤ 1 h/Mpc.Also, even those higher-N µ,AI power spectra with k st < 1 h/Mpc agree with the N µ,AI = 5 power to better than 1% across their entire stable ranges k ≤ k st .Evidently from the figure, our stabilization procedure discards higher k early enough to prevent them from contaminating k ≤ k st at even the percent level.
Out of the nineteen MT4 design points requiring a mode-coupling truncation N µ,AI < 2N µ,NL − 1, fifteen reach k = 1.2 h/Mpc with N µ,AI = 5, and the remaining four reach that wave number with N µ,AI = 4.These four N µ,AI = 4 models all have Ω ν,0 h 2 > 0.008, with a mean Ω ν,0 h 2 value of 0.0092.All have stability thresholds k st > 0.7 h/Mpc when run with N µ,AI = 5.By comparison with Fig. 1 (Bottom), we may estimate their truncation error as ∼ 1%.Another estimate of the truncation error, for all nineteen stabilized models, is the difference between the N µ,AI = 4 and N µ,AI = 5 power spectra, either to k = 1 h/Mpc or to k st if this is less than one.By this measure, we find that one of these nineteen has a 2.5% truncation error, while all others have ≤ 2% and nine of them have ≤ 1%.Since these nineteen represent nearly half of the forty design models with Ω ν,0 h 2 ≥ 0.0063, we estimate that truncation leads to a ≈ 1% power underestimate in that range.

Reduced power spectrum
The power spectrum of massive neutrinos declines sharply below the free-streaming scale, k ≫ k FS , giving ∆ 2 ν (k) a large dynamic range, which in turn makes it difficult to emulate.Our strategy is to divide the neutrino power spectra, the total power as well as the singledecile power spectra, by quickly-calculable linear approximations.Further reduction in the dynamic range is achieved by taking the natural logarithm of this power spectrum ratio, for each of the MT4 models, resulting in the training set P * im .As a starting point for a fast linear approximation to the neutrino power spectra, we take the matter power spectra of Eisenstein and Hu (Eisenstein & Hu 1998;Hu & Eisenstein 1998;Eisenstein & Hu 1997).In the clustering limit, k ≪ k FS , neutrinos and cold matter cluster very similarly.In the free-streaming limit, k ≫ k FS , the ratio of the neutrino and total matter density contrasts scales as k 2 FS /k 2 , as shown by Ringwald & Wong (2004) and Wong (2008).Those studies developed and tested an interpolation function, δ ν /δ m ≈ (1 + k/k FS ) −2 , that applies to linearly-clustering neutrinos, even when the CDM and baryons cluster non-linearly.Chen et al. (2021a,b) generalized the free-streaming scale of Ringwald & Wong (2004)  where, in the non-relativistic approximation, v(a) = τ/(am ν ), implying that k FS (a, v) ∝ √ a.Thus our Eisenstein-Hu-like approximations to the total and decile neutrino power spectra are respectively where ⟨v α ⟩ L in the final line denotes an average over all flows α making up decile L, and ζ(x) is the Riemann zeta function.
Next, we consider the range of scales over which we emulate the neutrino power spectra.We set the upper limit of our range of wave numbers to 1 h/Mpc, since beyond that we expect nonperturbative effects, such as the capture of neutrinos by halos, to dominate the power spectrum.Meanwhile, Fig. 2 demonstrates that our minimum wave number should not be much larger than ∼ 0.001 h/Mpc.Below that value, all power spectra shown agree with CAMB to ≤ 5%, with percent-level agreement for nearly all models below k = 0.002 h/Mpc.Above k = 0.001 h/Mpc, fractional differences rise rapidly, reaching ∼ 20% by k = 0.01 h/Mpc.Thus we choose to emulate log(∆ 2 ν /∆ 2 EH,ν ) and log(∆ 2 L /∆ 2 EH,L ) over the range 0.001 h/Mpc ≤ k ≤ 1 h/Mpc.

Non-linear enhancement: A first look
The chief goal of this study is to quantify the non-linear clustering of neutrinos across a wide range of parameters.One tool we will use for this is the non-linear enhancement ratio, that is, the ratio of the FlowsForTheMasses power spectrum to that from the Multi-Fluid Linear Response code MuFLR of Chen et al. (2021a): 5 MuFLR allows the CDM+baryon fluid to cluster non-linearly while limiting neutrino clustering to the linear terms.That is, the modecoupling integrals A α,acd,be f , hence also the bispectrum integrals I α,acd,be f and the non-linear correlations ξ αℓ , are set to zero, leaving only Eqs. (18,19,21) to be solved for the neutrinos.The enhancement ratio is instructive because it allows us to isolate directly the effects of non-linearity in the neutrino sector.High-quality simulations using neutrino linear response include Ali-Haimoud & Bird (2012), McCarthy et al. (2017), andLiu et al. (2018).R ν (k, z) quantifies the amount by which these underestimate neutrino clustering.Furthermore, Chen et al. (2023a,b) showed that FlowsForTheMasses itself becomes inaccurate in the regime R ν − 1 ≳ 1, indicating the scales on which a particle neutrino simulation is necessary for accurate predictions of neutrino clustering.Section 6 will use emulation to isolate the impact upon R ν of individual cosmological parameters.
Figure 3 provides a glimpse of the parameter-dependence of R ν (k, 0).From Fig. 3 (a) we immediately see that Ω ν,0 h 2 ∝ M ν has a significant effect upon R ν .Furthermore, non-linear corrections to the lightest neutrinos reduce R ν slightly below unity.This is not surprising, as non-linear corrections in the weakly non-linear regime are known to suppress even the CDM clustering (Bernardeau et al. 2002).Evidently from Fig. 3, the clustering amplitude σ 8 also has a discernible effect upon R ν , while the impact of the total matter and baryon fractions are less obvious.We will revisit the impact of different cosmological parameters upon R ν (k, 0) in Sec.6.1.

Emulation using SEPIA
For each neutrino momentum decile L, we use the SEPIA code of Gattiker et al. (2020) to determine the values of the principal component basis functions ϕ (L) j (k i , z i ) as well as sample from the posteriors of the hyperparameters β (L)  jℓ , λ (L) U, j , and λ (L) W, j .We use N PC = 50 5 MuFLR is publicly available at github.com/upadhye/MuFLR .principal components for each L, comparable to the 45 used for the MT4 emulator.Hyperparameters are optimized in SEPIA using N MCMC = 50000 Markov chain Monte Carlo steps.Since the error on the total neutrino power spectrum for N MCMC of 10000, 20000, and 50000 is, respectively, 3.92%, 3.86%, and 3.48%, measured against ten randomly-chosen models outside of the MT4 sample, we regard N MCMC = 50000 to have converged.
SEPIA is designed to draw weights from the probability distribution of Eq. ( 10).We instead prefer a deterministic emulator such as that described in Heitmann et al. (2009), which uses the mean weights W(L) j ( ⃗ C ) of Eq. ( 11) as the emulator prediction.These depend on the hyperparameters' posterior means β(L) jℓ , λ(L) U, j , and λ(L) W, j .Appendix A shows how to obtain these quantities from SEPIA.
Our goal is now in sight: the mean weight W(L) j ( ⃗ C ) of Eq. ( 11), for each decile L and principal component j, for a given cosmological parameter vector ⃗ C. The final ingredient needed to emulate W(L) With fixed L and j, the quantity R * (L) Since the Kriging basis is independent of ⃗ C, we compute it once and save the result.Now, given ⃗ C, we may readily compute r * jm ( ⃗ C ) of Eq. ( 13) for each m = 0, . . ., N M − 1.The dot product of the vector ⃗ r * j ( ⃗ C ) with the Kriging basis gives W(L) j ( ⃗ C ), as in Eq. ( 11).The emulated reduced neutrino power ) for decile L is given by: and the total neutrino power spectrum ∆ 2 ν (k i , z i ) is formed by averaging over the individual decile powers ∆ 2 L (k i , z i ): As an alternative, we could have emulated ∆ 2 ν (k i , z i ) separately.However, we find averaging over decile powers to be more accurate.Further, averaging ensures consistency between ∆ 2 ν (k i , z i ) and the individual ∆ 2 L (k i , z i ), which should obey Eq. (34).

Tests of Cosmic-Eν
Now that our Cosmic-Eν emulator is complete, we quantify its accuracy.We first consider the total neutrino power spectrum ∆ 2 ν (k, z), the main goal of this study.Power spectra of individual momentum deciles, ∆ 2 L (k, z), are covered at the end of this subsection.We begin with the most accurate error measurement, which compares emulator predictions to FlowsForTheMasses computations for a set of ten test models outside of the MT4 design set.Table 2 lists the cosmological parameters of these out-of-sample models, E001 through E010.They cover a large range of Ω ν,0 h 2 from 0.0007,  2. For each model, lines show the redshifts 3, 1, and 0, in order of increasing thickness.The gray shaded region is the mean plus-or-minus one standard deviation, maximized over all emulated redshifts.Across all k and z, this error is less than 3.5%.Lines show the holdout-to-FlowsForTheMasses power spectrum ratios at z = 0, with colors corresponding to Ω ν,0 h 2 .The dark gray shaded region is the mean plusor-minus one standard deviation, maximized over all emulated redshifts; the maximum 1σ error across all k and z is 5.7%.
for E006, to 0.009, for E005, and allow for a substantial variation in the dark energy equation of state.Figure 4 compares the total neutrino power spectra of Cosmic-Eν and FlowsForTheMasses for the ten out-of-sample models of Table 2.For each of the 27 redshifts emulated and at each k, the mean and standard deviation of the Cosmic-Eν-to-FlowsForTheMasses ratio are computed.The 1σ error interval is the region within one standard deviation of the mean ratio at each k and z emulated, i.e., the standard deviation across cosmologies for a single redshift and wave number.At each k, the gray shaded region shows the 1σ error interval maximized over redshift.The 1σ error, the difference between this region and unity, is everywhere less than 3.5%.The maximum error across all ten models, over all k and z, is 6.9%.The largest error is for model E006, which has a small neutrino density.
Another emulator error estimate, the leave-one-out holdout test, is performed directly through SEPIA.For each model m of the N M models in the training set, SEPIA builds an emulator with all N M − 1 models excluding m, then compares the resulting emulator predic-  Holdout tests for the two lowest-neutrino-mass models, which also have the largest errors in Fig. 5, color-coded by the scale factor a. Solid (dashed) curves correspond to Ω ν,0 h 2 = 0.00017 (Ω ν,0 h 2 = 0.00019).
holdout tests tend to overestimate the emulator error, particularly in the context of a space-filling design.
Figure 5 shows the results of leave-one-out holdout tests of Cosmic-Eν.The maximum 1σ error is now 5.7%, about 60% higher than the out-of-sample test.Nevertheless, at all but the largest scales, k < 0.002 h/Mpc, and the smallest scales, k > 0.8 h/Mpc, the error is 4% or less.Meanwhile, the largest holdout test error across all k and z, and all 101 models, is 25.4%.
The two models with the largest holdout test errors are further studied in Fig. 6.Out of the models in the MT4 training set, these two have the smallest Ω ν,0 h 2 ; they are the only models allowing Ω ν,0 h 2 to fall below half of its lower bound ≈ 0.00064 from neutrino oscillation experiments.Evidently, the 25% errors noted above are due to sharp falls in the Cosmic-Eν predictions at the earliest times and largest scales.Errors are under 20% for all k ≥ 0.0012 h/Mpc and all z.Since the holdout test errors of Fig. 5 overestimate the more accurate out-of-sample errors of Fig. 4 by about 60%, we estimate Cosmic-Eν errors of 10% − 12% at the lowest neutrino masses, in the range k ≥ 0.0012 h/Mpc.
Finally, we compare individual neutrino momentum deciles between Cosmic-Eν and FlowsForTheMasses.In multi-fluid perturbation theories such as ours of Sec.2.2, the neutrino flow speed v α behaves similarly to a sound speed, and the resulting power spectrum exhibits oscillatory behavior on sufficiently small scales.Averaging over a large number of flows eliminates these oscillations.However, each decile averages over only 5 flows, rather than 50 for the total neutrino power, making the small-scale decile powers noisier and more difficult to emulate.This is especially true for the higher-momentum deciles for which individual-flow oscillations are more prominent.Thus we expect less accuracy in ∆ 2 L (k, z), particularly for high L, than in the ∆ 2 ν (k, z) shown above.Figure 7 combines out-of-sample tests and leave-one-out holdout tests for each of the ten momentum deciles.The 1σ errors shown are in line with our expectations above.For deciles L = 0 and 1, errors from out-of-sample tests are < 7%, or about twice the maximum error on ∆ 2 ν (k, z).For higher L, errors remain at this level for k ≤ 0.1 h/Mpc but grow substantially at large k, rising to nearly four times the ∆ 2 ν (k, z) emulator error at k = 1 h/Mpc.

Variation of the ν implementation
The previous section quantified the precision with which Cosmic-Eν reproduced its underlying FlowsForTheMasses perturbation theory.Next, we consider its accuracy relative to numerical simulations of the massive neutrino power spectrum.Adamek et al. (2023) carried out an extensive comparison of neutrino simulation methods for a spatially-flat νΛCDM model with M ν = 0.15 eV (Ω ν,0 h 2 = 0.00161), Ω m,0 h 2 = 0.1432, Ω b,0 h 2 = 0.022, A s = 2.215×10 −9 (σ 8 = 0.815), h = 0.67, and n s = 0.9619.Though all of the neutrino power spectra agree at about the percent level in the linear regime, k ≲ 0.1 h/Mpc, the different methods are discrepant at the 30% − 40% level by k = 1 h/Mpc.Even within a given simulation method, choices of initial conditions and simulation parameters have a substantial impact on the massive neutrino power spectrum.Sullivan et al. (2023) implemented the "tiling" method of Banerjee et al. (2018).They find that the number of allowed neutrino momentum directions is the most important parameter for determining convergence, and they judge their highest-direction-number run to have converged at the ∼ 10% level up to half of the neutrino Nyquist frequency, or ∼ 1 h/Mpc in their standard run.Combining this with the discrepancies among the different simulation methods, we may regard simulation uncertainty in the neutrino power at k ≲ 1 h/Mpc to be in the 30% − 50% range.
Figure 8 compares Cosmic-Eν and MuFLR linear response to the results of Adamek et al. (2023).Each power spectrum is compared to that computed using the SWIFT simulation of Schaller et al. (2023a), with their ratio smoothed using a centered 10-point moving average.Compared with SWIFT, the CONCEPT code of Dakin et al. (2019) predicts 22% more power at k = 1 h/Mpc, while gevolution (Adamek et al. 2016a) predicts ≈ 15% − 20% less.
Cosmic-Eν agrees closely with all of the simulated ∆ 2 ν at low k.Its predicted power spectrum falls below that of SWIFT by ≤ 3% up to k = 0.17 h/Mpc, or three times the free-streaming wave number for  15 eV, at z = 0, computed using a variety of methods by Adamek et al. (2023).Solid (dashed) lines represent power spectra that are greater (less) than that of SWIFT.Ratios have been smoothed using a centered 10-point moving average.this neutrino mass.This power deficit grows rapidly, rising to 19% at k = 0.4 h/Mpc and 49% at k = 1 h/Mpc, somewhat greater than, but comparable to, the scatter between different N-body simulation methods.Thus FlowsForTheMasses and the Cosmic-Eν emulator appear to provide accurate computations of the non-linear neutrino power, given the current level of simulation uncertainty.Most of the scatter among the N-body methods at wave numbers k ≈ 0.5 h/Mpc is due to a systematic difference between particlebased methods, such as SWIFT, and CONCEPT, which integrates the massive neutrino fluid equations on a grid.A priori we have no reason to consider one of these more accurate.However, if we exclude CONCEPT as an outlier among the N-body simulations, then the range spanned by the remaining simulations drops significantly, and the Cosmic-Eν power deficit exceeds this range by a factor of a few for k ≳ 0.5 h/Mpc.Thus it is not clear whether the Cosmic-Eν small-scale power deficit is due to systematic uncertainties among the different non-linear methods or to a genuine non-perturbative effect such as the capture of neutrinos by CDM+baryon halos.
Also shown in Fig. 8 is the MuFLR linear response power spectrum.As with Cosmic-Eν, its power is less than that of the SWIFT simulation used as a reference.Its power deficits relative to SWIFT are significantly larger than those of Cosmic-Eν, 6.1 times larger at k = 0.1 h/Mpc and 1.8 times larger at k = 0.4 h/Mpc.Thus Cosmic-Eν represents a significant accuracy improvement over linear response approximations.

Rapidly-evolving dark energy
Since neutrino mass bounds are dependent upon constraints on the growth factor of large-scale structure, they are degenerate with variations in the dark energy equation of state.For example, Upadhye (2019) found a factor-of-three degradation in the 95%-confidence M ν bound when w 0 and w a were allowed to vary.In recognition of this degeneracy, MT4 and Cosmic-Eν allow for substantial variations in w 0 and w 0 + w a .Here we test the accuracy of Cosmic-Eν for such a rapidly-varying equation of state by comparison to a νwCDM N-body simulation.
Figure 10 tests Cosmic-Eν for this νwCDM simulation, with a rapidly-varying equation of state, at redshifts 0 and 1. Inner and outer shaded bands show the regions within 20% and 50% of the N-body power, respectively.The accuracy of Cosmic-Eν is in line with our previous νΛCDM comparisons to the SWIFT simulations: 20% up to k = 0.3 h/Mpc − 0.4 h/Mpc, and 50% up to k ≈ 1 h/Mpc, with slightly higher accuracy at z = 1.Thus we conclude that even |w a | ∼ 1 does not diminish the accuracy of Cosmic-Eν.
This Section has quantified the accuracy of the Cosmic-Eν ∆ 2 ν (k, z) emulator across a wide range of M ν , for a cosmological constant as well as a rapidly-evolving equation of state, by comparison to N-body simulations using a few very different massive neutrino simulation methods.Its error at z = 0 is fairly consistent across a wide range of methods: a few percent up to k ≈ 0.15 h/Mpc, an ≈ 20% power underestimate at k = 0.4 h/Mpc, and an ≈ 50% underestimate at k = 1 h/Mpc.Since ∆ 2 ν (k, z) for k ≫ k FS scales approximately as M 4 ν (Ringwald & Wong 2004;Wong 2008), these under-estimates at k = 0.4 h/Mpc and k = 1 h/Mpc are consistent with 5% and 13% biases in M ν , respectively.Before proceeding, we comment upon the discrepancy between the 14% error in FlowsForTheMasses reported in Chen et al. (2023b), over the entire range k ≤ 1 h/Mpc, and the larger differences with N-body simulations evident in Figs. 8, 9, and 10 for k ≈ 1 h/Mpc.There are two possibilities: errors in the hybrid simulations of Chen et al. (2023a) used to test FlowsForTheMasses, and a systematic error causing the discrepancies between differing N-body implementations of neutrinos.
Errors in Chen et al. (2023a) may be due to a finite number of neutrino flows, residual shot noise, and a finite simulation volume.Finite-flow-number errors in the linear response calculations of Chen et al. (2021a) were found to be ≈ 10%, and non-linear neutrino clustering likely increases them somewhat.Of course, increasing the number of flows, or sampling the Fermi-Dirac distribution more efficiently, will improve the accuracy of FlowsForTheMasses as well as the simulations.While the estimated simulation shot noise (kL sim ) 3 /(2π 2 N sim ), for N sim particles in a volume L 3 sim , was subtracted from ∆ 2 ν (k), residual shot noise remains.In an effort to mitigate shot noise, Chen et al. (2023a) chose a small box, L sim = 128 Mpc/h, at the cost of neglecting the contributions of larger modes to small-scale non-linear growth.
Meanwhile, the 30% − 40% spread among the different simulations in Fig. 8 at k = 1 h/Mpc suggests small-scale systematic errors in some of these methods.The gevolution power spectrum of Adamek et al. (2016a) is about 15% − 20% lower than SWIFT at k = 1 h/Mpc, meaning that correcting the 14% underestimate of Cosmic-Eν relative to Chen et al. (2023a), as well as the ≈ 10% underestimate due to a finite number of flows, would put Cosmic-Eν within 5% − 10% of gevolution.Moreover, in the k ≲ 0.2 h/Mpc range where the N-body methods of Fig. 8 agree with one another to a few percent, Cosmic-Eν also agrees with them at that level.The SWIFT-Cosmic-Eν difference rises along with the CONCEPTgevolution difference in Fig. 8. Thus we cannot conclusively attribute the discrepancy between the 14% FlowsForTheMasses error estimate of Chen et al. (2023b) and the 49% SWIFT-Cosmic-Eν difference to errors in Chen et al. (2023a).

Parameter-sensitivity of the enhancement ratio
Now that we have quantified the accuracy of Cosmic-Eν, we may use it to study the non-linear clustering of massive neutrinos.We focus here on the non-linear enhancement ratio R ν (k, z) of Eq. ( 31), that is, the ratio of the neutrino power spectra using FlowsForTheMasses and MuFLR, with the CDM+baryon treatment held fixed.We emulate R ν by taking the ratio of Cosmic-Eν to a MuFLR emulator.
Figure 11 compares perturbative (dashed) and emulated (solid) calculations of R ν at z = 0 for the out-of-sample models of Table 2.A couple trends are evident.Firstly, the total neutrino power is typically more accurate at high k than individual decile powers.At k = 1 h/Mpc, the emulated R ν agrees with the FlowsForTheMasses computation to better than 2% for eight of the ten models.For decile 0, this error rises to 3.2%, and for decile 1 to 4.8%.
Secondly, the lower deciles and higher neutrino masses tend to have smaller errors.This is due to the fact that larger L and smaller masses lead to larger average velocities, hence more prominent oscillatory behavior in the free-streaming limit, making these flows . Sensitivity of R ν (k, z) at z = 0 to variations of the eight cosmological parameters shown in Table 1.difficult to emulate.Both of these trends are consistent with the individual-decile out-of-sample and holdout tests of Sec.4.2.
Next, we consider the sensitivity of R ν to the cosmological parameters.As a fiducial model at which to test this sensitivity, we choose a ΛCDM model in which each parameter except for w 0 and w a is set to the midpoint of its range in Table 1. Figure 12 shows the derivative of log R ν with respect to each parameter about this fiducial model.We have checked that the results are qualitatively similar for Ω ν,0 h 2 = 0.002.
Above k ≈ k FS = 0.16 h/Mpc, the dominant effect is a rise in R ν with Ω ν,0 h 2 ∝ M ν .This can be understood by noting that, with all parameters in Table 1 other than Ω ν,0 h 2 held fixed, the small-scale linear and linear-response neutrino power scales as (Ω ν,0 h 2 ) 4 , while the non-linear power rises relative to the LR power.We can estimate the non-linear enhancement to the linear scaling law from Fig. 12: The final term on the right, the logarithmic derivative of R ν , is Ω ν,0 h 2 times d ln(R ν )/d(Ω ν,0 h 2 ) from the figure, or about 0.57 at k = 1 h/Mpc ≈ 6k FS .A similar analysis for Ω ν,0 h 2 = 0.002 finds this scaling enhancement to be 0.35 at k = 0.4 h/Mpc ≈ 6k FS and 0.59 at k = 1 h/Mpc ≈ 15k FS , suggesting a rise in this scaling enhancement with both Ω ν,0 h 2 and k/k FS .Section 6.2 will explore this enhancement further.At k ≳ k FS = 0.16 h/Mpc, the next most significant parameter for determining R ν is the physical baryon fraction Ω b,0 h 2 .The oscillatory nature of ∂ ln R ν /∂(Ω b,0 h 2 ) for k ≲ 0.8 h/Mpc, and its decline for k ≳ 0.9 h/Mpc, suggests that Ω b,0 h 2 affects R ν primarily through the baryon acoustic oscillations (BAO).Modifications to the BAO, in turn, are amplified by the non-linear clustering of neutrinos.
After Ω ν,0 h 2 and Ω b,0 h 2 , the next parameters to which R ν is most sensitive at k ≳ k FS are σ 8 and Ω m,0 h 2 .At first glance, the relative sensitivities of R ν to σ 8 , Ω m,0 h 2 , and Ω b,0 h 2 appear to contradict Fig. 3.However, the range of R ν associated with each parameter in that figure is the derivative in Fig. 12 times the parameter range in Table 1.This range for σ 8 is about six times larger than for Ω m,0 h 2 , and 100 times larger than for Ω b,0 h 2 .
Thus, given a particular data combination, we may define an alternative sensitivity measure for each parameter p by multiplying ∂ log R ν /∂p from Fig. 12 by the range ∆p allowed by the data.As an example, we choose the νwCDM analysis of Upadhye (2019), constrained using a combination of CMB, galaxy, and supernova data, and marginalized over a five-parameter model of scale-dependent  galaxy bias.We approximate the 95% CL interval of Ω m,0 h 2 by the corresponding one for the CDM alone, its dominant component.Table 3 shows the result at two wave numbers.While Ω ν,0 h 2 remains by far the most significant parameter for determining R ν , σ 8 is also important.In summary, while R ν is most sensitive to Ω ν,0 h 2 , σ 8 , Ω m,0 h 2 , and Ω b,0 h 2 , the first two of these are the most important given current parameter constraints.

Relative clustering in the free-streaming limit
Next, we consider further the mass scaling of the neutrino power spectrum in the free-streaming limit, raised in the previous subsection.Ringwald & Wong (2004) argues that neutrino linear response to non-linear CDM+baryon growth results in a scaling ∆ 2 ν /∆ 2 m ∝ M 4 ν , while the Tremaine-Gunn bound of Tremaine & Gunn (1979);Shu (1978Shu ( , 1987)); Kull et al. (1996) Ringwald & Wong (2004) demonstrate using N-body simulations that halos approach this latter bound, and, further, that the bound can be exceeded, especially in the case of small M ν , if one includes all neutrinos present, rather than only those captured by the halo's gravitational potential.
Thus we may expect ∂ log(∆ 2 ν /∆ 2 m )/∂ log(M ν ) to rise above four while remaining below six. Figure 13 numerically differentiates the ratio of the neutrino power spectrum to the MT4 total-matter power spectrum with respect to Ω ν,0 h 2 ∝ M ν with a step size of 1% in Ω ν,0 h 2 ∝ M ν .Solid and dashed lines respectively use non-linear (Cosmic-Eν) and linear response (MuFLR) neutrino power spectra.
Consider first the larger wave numbers.The MuFLR curves approach 4 from below but never exceed it, as expected.Meanwhile, the Cosmic-Eν non-linear logarithmic derivatives for k = 0.5 h/Mpc and k = 1 h/Mpc both exceed 4 for small Ω ν,0 h 2 , where these wave numbers are many times the free-streaming scale.The non-linear enhancement to the mass scaling in Eq. ( 35) is the difference between the solid and dashed lines.Focusing on k = 1 h/Mpc, we find this to be 0.67 for Ω ν,0 h 2 = 0.002 and 0.58 for Ω ν,0 h 2 = 0.005.Since a 3.5% emulator error implies an error of ∼ 0.1 in the logarithmic derivative, these are consistent with the results of Sec.6.1.
This emulator error means that the difference between the linear and non-linear curves for k = 0.1 h/Mpc is consistent with zero.The same is true for k = 0.2 h/Mpc for Ω ν,0 h 2 ≳ 0.003.Further, the small oscillations observed in some of the logarithmic derivatives are consistent with emulator fluctuations.Thus Sections 6.1-6.2 consistently demonstrate a non-linear enhancement of ≳ 0.5 to the power law scaling ∂ log(∆ 2 ν )/∂ log(Ω ν,0 h 2 ) at k = 1 h/Mpc for 0.002 ≲ Ω ν,0 h 2 ≲ 0.005.

Neutrino contribution to the matter power
The MT4 emulator includes fully linear neutrinos, as implemented in the CAMB code of Lewis et al. (2000); Lewis & Bridle (2002), in their CDM+baryon and total matter power spectra, as described in Saito et al. (2008); Agarwal & Feldman (2011); Upadhye et al. (2014Upadhye et al. ( , 2016)).Since the non-linear clustering of neutrinos increases their power by an order of magnitude relative to linear theory, as shown in Fig. 2, we quantify here the impact of neutrino nonlinearity on the matter power spectrum.Neutrinos will affect ∆ 2 m (k, z) in two ways: indirectly, by adding to the gravitational potential, hence enhancing CDM+baryon clustering; and directly, through their inclusion in the total matter power.
Quantifying the indirect effect precisely, by incorporating FlowsForTheMasses into an N-body simulation, is beyond the scope of this study.However, we may bound this effect.Chen et al. (2021a) carries out an N-body simulation with neutrino linear response through the MuFLR code.For the largest neutrino fraction considered here, Ω ν,0 h 2 = 0.01, that study finds an indirect enhancement of 0.05% to the CDM+baryon power spectrum.Since the nonlinear enhancement ratios R ν are less than 5 in Figs. 3, 11, the indirect enhancement is < 0.25%.A more accurate estimate directly multiplying the linear response enhancement of Chen et al. (2021a) by R ν for their model finds an indirect enhancement of 0.16%.
Figure 14 quantifies the direct effect, which is everywhere less than one percent.Even this is a slight overestimate, as FlowsForTheMasses, hence Cosmic-Eν, assume that CDM, baryon, and neutrino density-contrast monopoles are perfectly correlated.Under this approximation, the matter power spectrum is (36) Bird et al. (2018) shows that the actual neutrino-CDM correlation function drops below unity by ≤ 4% for k ≤ 0.5 h/Mpc and ≤ 8% for k ≤ 1 h/Mpc.The smallness of this deviation is due to the fact that the initially-slowest neutrinos, which contribute the most to smallscale clustering, are also closely correlated with the CDM.This correlation being slightly less than one implies that the actual direct contribution of neutrinos to ∆ 2 m is slightly smaller than in Fig. 14. m (k, z), computed using the MT4 CDM+baryon power and the Cosmic-Eν neutrino power as in Eq. ( 36), is divided by the MT4 emulated total matter power ∆ 2 m (k, z).

CONCLUSIONS
FlowsForTheMasses, the first non-linear perturbative power spectrum calculation for free-streaming particles such as massive neutrinos, provides detailed information on the clustering of neutrinos of different initial momenta.We have emulated the total nonlinear neutrino power spectrum as well as separate power spectra for the ten momentum deciles, each representing a tenth of the neutrino number density.Our emulated ∆ 2 ν (k, z) agrees precisely with FlowsForTheMasses to < 3.5% for 10 −3 h/Mpc ≤ k ≤ 1 h/Mpc and 0 ≤ z ≤ 3, as shown in Fig. 4. Individual-decile errors range from about twice as large for the lowest momenta to four times as large for the fastest-moving neutrinos with highly oscillatory density contrasts; see Fig. 7.We have released our emulator as Cosmic-Eν.
Comparing Cosmic-Eν to the highest-resolution simulations of Adamek et al. (2023) in Fig. 9, we found agreement to 3% up to k = 3k FS = 0.17 h/Mpc and 19% to k = 0.4 h/Mpc.Above this wave number, Cosmic-Eν increasingly underpredicts the simulations of Adamek et al. (2023), with this underprediction reaching nearly 50% by k = 1 h/Mpc.Even this error is not substantially larger than the 30% − 40% scatter between different simulation methods seen in Fig. 8, so we cannot definitively attribute it either to a nonperturbative effect beyond the capabilities of FlowsForTheMasses or to a systematic error in the simulations.Importantly, Cosmic-Eν provides a neutrino power spectrum in about ten milliseconds on a standard desktop machine, and we have confirmed that its accuracy is unaffected by rapid variations in the dark energy equation of state.
One strength of the emulation technique is our ability to differentiate numerically the emulated function without the result being dominated by the shot noise and sample variance affecting N-body power spectra.Section 6 took full advantage of this capability by studying the non-linear enhancement ratio R ν (k, z) of Eq. ( 31) and the neutrino-to-matter ratio ∆ 2 ν /∆ 2 m .Differentiating R ν with respect to each of the cosmological parameters, we find that it is most sensitive to the physical neutrino density Ω ν,0 h 2 , but also to Ω b,0 h 2 , σ 8 , and Ω m,0 h 2 .Furthermore, we demonstrated a non-linear enhancement of ≈ 0.5 to the free-streaming-limit scaling ∂ log(∆ 2 ν /∆ 2 m )/∂ log(M ν ) → 4, meaning that non-linear clustering makes the small-scale density of neutrinos even more sensitive to their mass.Our results demon-strate the speed and efficacy of the emulation technique in neutrino cosmology.
to individual neutrino flows through the replacement of their neutrino sound speed c ν by the flow velocity v, k FS(a, v)

Figure 2 .
Figure2.Fractional difference between FlowsForTheMasses and CAMB for ten cosmologies spanning the parameter space, listed in Table2.
length N M .Since R * (L) j and ⃗ w * (L) j are both known, we find the Kriging basis by solving the linear system: R *

Figure 4 .
Figure 4. Cosmic-Eν vs. FlowsForTheMasses for the total neutrino power spectra ∆ 2 ν (k, z) of the out-of-sample models of Table2.For each model, lines show the redshifts 3, 1, and 0, in order of increasing thickness.The gray shaded region is the mean plus-or-minus one standard deviation, maximized over all emulated redshifts.Across all k and z, this error is less than 3.5%.

Figure 5 .
Figure5.Leave-one-out holdout tests of Cosmic-Eν.Lines show the holdout-to-FlowsForTheMasses power spectrum ratios at z = 0, with colors corresponding to Ω ν,0 h 2 .The dark gray shaded region is the mean plusor-minus one standard deviation, maximized over all emulated redshifts; the maximum 1σ error across all k and z is 5.7%.

Figure 9 comparesFigure 9 .
Figure 9.Comparison of Cosmic-Eν (solid) and MuFLR (dashed) to GADGET-3 N-body neutrino power spectra of Adamek et al. (2023), at z = 0, for neutrino masses M ν ranging from 0.15 eV to 0.6 eV.For M ν = 0.15 eV, the high-resolution 1024 3 -particle simulations of that reference are used.(Top) N-body power spectra are shown as filled circles for (512 Mpc) 3 simulation volumes and open circles for the (1024 Mpc) 3 M ν = 0.15 eV simulation.(Bottom) Fractional errors in Cosmic-Eν and MuFLR compared with the (512 Mpc) 3 -box simulations at k ≥ 0.1 h/Mpc and the (1024 Mpc) 3 -box simulation below that wave number.

Table 1 .
Allowed ranges of each cosmological parameter in Cosmic-Eν.As with the MT4 emulator, we have assumed three degenerate-mass neutrinos.

Table 2 .
Cosmological parameters for the ten out-of-sample test models.Each is a spatially-flat νwCDM model with w(a) = w 0 + (1 − a)w a .
tion P im at ⃗ C * m with the excluded training data P * im for that model.Since leaving out model m creates a gap in the training data at ⃗ C * m ,