Two-loop power spectrum with full time- and scale-dependence and EFT corrections: impact of massive neutrinos and going beyond EdS

We compute the density and velocity power spectra at next-to-next-to-leading order taking into account the effect of time- and scale-dependent growth of massive neutrino perturbations as well as the departure from Einstein--de-Sitter (EdS) dynamics at late times non-linearly. We determine the impact of these effects by comparing to the commonly adopted approximate treatment where they are not included. For the bare cold dark matter (CDM)+baryon spectrum, we find percent deviations for $k\gtrsim 0.17h~\mathrm{Mpc}^{-1}$, mainly due to the departure from EdS. For the velocity and cross power spectrum the main difference arises due to time- and scale-dependence in presence of massive neutrinos yielding percent deviation above $k\simeq 0.08, 0.13, 0.16h~\mathrm{Mpc}^{-1}$ for $\sum m_{\nu} = 0.4, 0.2, 0.1~\mathrm{eV}$, respectively. We use an effective field theory (EFT) framework at two-loop valid for wavenumbers $k \gg k_{\mathrm{FS}}$, where $k_{\mathrm{FS}}$ is the neutrino free-streaming scale. Comparing to Quijote N-body simulations, we find that for the CDM+baryon density power spectrum the effect of neutrino perturbations and exact time-dependent dynamics at late times can be accounted for by a shift in the one-loop EFT counterterm, $\Delta\bar{\gamma}_1 \simeq - 0.2~\mathrm{Mpc}^2/h^2$. We find percent agreement between the perturbative and N-body results up to $k\lesssim 0.12h~\mathrm{Mpc}^{-1}$ and $k\lesssim 0.16h~\mathrm{Mpc}^{-1}$ at one- and two-loop order, respectively, for all considered neutrino masses $\sum m_{\nu} \leq 0.4~\mathrm{eV}$.

The sum of neutrino masses is known to be greater than 0.06 eV from oscillation experiments [54], and β-decay experiments sets an upper bound m β < 0.8 eV at 90% C.L. for the effective electron anti-neutrino mass [55]. Cosmological probes offer complementary constraints because of observable impacts left by the neutrinos on the cosmological evolution. In particular, mainly due to late ISW and lensing effects, CMB measurements constrain the sum of neutrino masses to M ν ≡ m ν ≤ 0.26 eV at 95% C.L. [56], improving to M ν ≤ 0.12 eV when combined with data on the scale of baryon acoustic oscillations (BAO). In addition, the large neutrino velocity dispersion after their non-relativistic transition in the late Universe slows down structure formation on scales smaller than the neutrino free-streaming length, corresponding to a wavenumber k FS ∼ 0.01h Mpc −1 [9]. This characteristic impact is expected to be detectable by the Euclid satellite, yielding a forecasted measurement of the neutrino mass sum with σ(M ν ) 0.02-0.03 eV [9][10][11][12][13][14][15][16][17][18][19]. On the other hand, the presence of the neutrino free-streaming scale introduces a scale-dependence in the clustering dynamics, making the modelling of neutrino perturbations beyond linear theory complicated.
In order to model the mildly non-linear scales in a robust and efficient way as necessary for MCMC analyses, certain approximations have to be made. Accordingly, the approximations should be scrutinized both in ΛCDM and extensions of it to determine the theoretical uncertainty. In this work we perform a precision calculation of the matter power spectrum with the aim of examining the accuracy of certain common approximations. In particular, our analysis involve the following: • computing the matter density and velocity power spectrum at next-to-next-to-leading order (NNLO) in perturbation theory, using the method that can capture scale-and time-dependent dynamics introduced in [57] (see also [58][59][60]), • including neutrino perturbations up to fifth order (i.e. at two-loop) and their impact on the cold dark matter (CDM) and baryons fluid via gravity, utilizing a hybrid two-component fluid scheme introduced in [61] and further developed in [57] 1 , • relaxing the Einstein-de-Sitter (EdS) approximation [39,57,[66][67][68][69][70] and taking into account the exact ΛCDM (+M ν ) time-dependence at two-loop order (as first done in [57], see also [70]), • extending the work of [57] by including EFT corrections to the power spectrum in the twocomponent fluid model, promoting the approach of [71] to a cosmology with scale-dependent dynamics, and further resumming the effect of large bulk flows in the IR [25][26][27][28][29][30], • comparing theory predictions to and calibrating EFT parameters using N-body data from the Quijote simulation suite [72] for three cosmologies with M ν = 0.1, 0.2 and 0.4 eV.
Our work extends previous studies [16,[73][74][75] for massive neutrinos within the EFT context in evaluating the power spectrum at NNLO (two-loop) order, and including the full scale-and timedependence of neutrinos imprinted on non-linear kernels. The EFT setup used in this work is based on [71], and we argue that it can be extended to the case of massive neutrinos in the limit k k FS 0.01h Mpc −1 , being satisfied for current and future galaxy surveys. Nevertheless, we stress that our unrenormalized results capture also the scale dependence in the transition region k ∼ k FS .
With these elements, we can perform an assessment of the accuracy of neglecting neutrino perturbations as well as the EdS approximation on the power spectrum at two-loop order: • We first compare the unrenormalized CDM+baryon power spectrum between the full solution including non-linear neutrino perturbations and departure from EdS (which we name the 2F scheme) to the simplified treatment with linear neutrinos and EdS dynamics for CDM+baryons (named 1F scheme). We find deviations beyond one percent for k 0.17h Mpc −1 at z = 0, mainly arising from the departure from EdS.
• As a proxy for redshift space distortion effects, we analyze the neutrino mass dependence of the velocity divergence power spectrum and the cross spectrum. We indeed find deviations between the full 2F and approximate 1F scheme that strongly depend on the value of the neutrino mass, and reach more than one percent for k 0.08, 0.13, 0.16h Mpc −1 for either the cross-or velocity spectrum and M ν = 0.4, 0.2, 0.1 eV, respectively. Therefore, the 1F approximation is not sufficient to access the neutrino mass information encoded in redshift space [75,76].
• After including EFT corrections, we compare both 2F and 1F results to N-body data, and determine to what extent the discrepancy of the simplified treatment can be absorbed by a readjustment of counterterms at NNLO. For the CDM+baryon density spectrum we find that this is indeed the case to high accuracy, with a shift mainly in the one-loop EFT term by ∆γ 1 Figure 1. CDM+baryon density power spectrum predictions at linear, next-to-leading-and next-to-nextto-leading order in perturbation theory, normalized to Quijote N-body results, for three neutrino masses Mν = 0.1, 0.2 and 0.4 eV. The NLO curve correspond to the 0-parameter model and the NNLO curve to the 2-parameter model defined in Sec. 5, and we use a pivot scale kmax = 0.148h Mpc −1 . The gray area indicates uncertainty from the N-body simulation.
We also address the level of degeneracy between neutrino free-streaming suppression and the overall amplitude of fluctuations [76] for the power spectrum within perturbation theory. We note that the approach presented here can be extended to the (one-loop) bispectrum in a straightforward manner. It has been shown that this additional information could be instrumental to break the degeneracy of neutrino masses with the overall power spectrum amplitude [17,77].
A main result of this work is shown in Fig. 1, where we plot the power spectrum computed in the full perturbative solution, incorporating neutrino perturbations beyond the linear order as well as the departure from EdS. The results are normalized to the N-body results, and we show the linear, NLO and NNLO perturbative predictions for three neutrino masses. It is clear that adding higher order corrections extends the wavenumber range with percent accuracy. In particular, the NLO and NNLO results deviates by more than 1% for k 0.12h Mpc −1 and k 0.16h Mpc −1 in all cosmologies, respectively. We refer to Sec. 5 for further discussion, including potential limitations related to overfitting.
Our work is structured as follows: in the next section we describe how we capture neutrino perturbations in a two-component fluid model and compare it to a simplified approach for the (unrenormalized) density and velocity power spectra in Sec. 3. We detail how we renormalize the power spectrum using an EFT approach in Sec. 4 and compare and calibrate both models to N-body simulations in Sec. 5. We conclude in Sec. 6.

Perturbative modelling with massive neutrinos
In this section we briefly review Standard Perturbation Theory (SPT) and the generic extension of it introduced in [57] that can describe cosmologies with multiple species in addition to scale-and timedependent dynamics. We will apply this framework to cosmologies with massive neutrinos, therefore we next summarize a two-component fluid setup of CDM+baryons and neutrinos introduced in [61] and further developed in [57], which can readily be captured by the extension of SPT. Finally, we compare the two-component fluid model at the linear level to the full neutrino Boltzmann hierarchy.

Standard perturbation theory
The equations of motion for the density contrast δ and velocity divergence θ = ∂ i v i (neglecting vorticity) in Fourier space reads where τ is conformal time, H = d ln a/ dτ the conformal Hubble rate and Ω m the time-dependent matter density parameter. We introduced the shorthand notations k 12 = k 1 + k 2 and k = d 3 k and used δ D to denote the Dirac delta function. The mode coupling functions are as usual. In SPT one assumes that the anisotropic stress of the fluid vanishes. We also make this assumption here initially, but will relax it when discussing an effective field theory setup in Sec. 4. The equations of motion can be written in a compact form after defining the tuple ψ = (δ, −θ/Hf ) and using η = log D, with D and f being the linear growth factor and growth rate, respectively, thus [20]: The matrix Ω ab governing the linear evolution is given by 4) and the only non-zero components of the non-linear vertex are γ 121 (k, k 1 , k 2 ) = α(k 1 , k 2 ) and γ 222 (k, k 1 , k 2 ) = β(k 1 , k 2 ). In SPT one typically adopts the Einstein-de-Sitter (EdS) approximation, in which Ω m /f 2 = 1 so that the Ω ab -matrix becomes time-independent. This approximation greatly simplifies Eq. (2.3), allowing for analytic solutions order by order. Only the decaying mode is affected by changes in the ratio Ω m /f 2 and moreover the ratio only departs significantly from one at late times, z 2 in ΛCDM (and also in moderate extensions). Consequently, the EdS-approximation has been shown to work at the percent-level for the power spectrum [57,[66][67][68]70] as well as the bispectrum [39,69]. In particular, in EFT analyses the departure from EdS can be largely degenerate with counterterms, only leading to a shift in the EFT parameters [39,68]. In this work we consider schemes with and without the EdS approximation, which we will specify accordingly.

Extension of SPT
Following [57], we extend SPT by allowing for multiple species in the fluid as well as allowing for a general time-and wavenumber-dependence. More precisely, for an N -fluid we collect the density contrast and velocity divergence for each component i into into the field vector ψ a = (. . . , δ i , θ i , . . . ) with the index a running from 1 to 2N . In addition, we permit a general dependence on time and wavenumber in the (now 2N × 2N ) matrix describing the linear evolution Ω ab = Ω ab (|k|, η). This extension can capture multiple models beyond ΛCDM in addition to effective models of clustering dynamics. It has in general no analytic solution however, hence we will mostly need to solve the dynamics numerically.
The equations of motion (2.3) can also in this case be solved perturbatively, at each order furnished by 2N kernels F (n) a labeled by the index a at order n: a (q 1 , . . . , q n ; η) δ 0 (q 1 ; η ini ) · · · δ 0 (q n ; η ini ) , (2.5) where ∆η ≡ η − η ini and δ 0 is an initial condition that we discuss shortly. Note that due to the assumed non-trivial time-dependence in the dynamics, we allow for a dependence on η in addition to the wavenumbers q for the kernels. Inserting this solution into Eq. (2.3) yields the following recursive solution at n-th order in perturbation theory: .
Here, k = i q i , and the right hand side is understood to be symmetrized with respect to all permutations exchanging momenta in the {q 1 , . . . , q m } set with momenta in the {q m+1 , . . . , q n } set and normalized to the number of permutations.
Note that setting N = 1 and using the Ω ab -matrix from Eq. (2.4) with Ω m /f 2 = 1 (EdS approximation) in Eq. (2.6), we recover in the limit η ini → −∞ the usual kernel recursion relations with the replacements F We still need to specify suitable initial conditions in order to solve the above equations. Taking η ini after recombination but long before non-linearities become important at the scales of interest, we assume that the initial conditions for each fluid component is correlated, so that which holds for adiabatic initial conditions. Furthermore, we assume that δ 0 is a Gaussian random field, so that we only need to specify the initial linear power spectrum δ 0 (k)δ 0 (k ) = δ D (k+k )P 0 (k) in order to compute correlations of ψ a 's. The initial linear kernels F a (k, η ini ) impose the relative normalization for each perturbation and wavenumber. Deep in the linear regime η ini → ∞ the higher order initial kernels can be set to zero. In practice however, using η ini after recombination, we find that those n > 1 initial conditions work poorly because they excite transient solutions that do not entirely decay by η = 0 (z = 0). We return to this issue below.
In short, at L-loop it consists of integrating over L loop momenta using Monte Carlo integration (with CUBA [78]) where at every integration point a set of (2L + 1)-order kernels needs to be evaluated using the recursion relation (2.6). In general, there is no analytic solution of Eq. (2.6), therefore we solve for the kernels numerically. We refer to [57] for further details on the algorithm.

Two-component fluid: CDM+baryons and massive neutrinos
We will employ the hybrid two-component fluid setup described in [57] to model structure formation in massive neutrino cosmologies. In the following, we repeat the main elements of this setup for convenience, but refer to [57] for details of the implementation.
In the hybrid two-component fluid model, the system is described linearly by the full Boltzmann hierarchy until some intermediate redshift z match , after which the evolution is mapped onto a twocomponent fluid, suitable for computing non-linear corrections (see also [61]). Baryons and CDM comprise jointly one fluid component, which is coupled to the second component, the neutrinos, via gravity. A fluid description of massive neutrinos is suitable at late times because the coupling to higher moments in the Boltzmann hierarchy is suppressed by powers of T ν /m ν . We may therefore follow only the lowest moments of the hierarchy for the neutrinos, with an effective sound velocity to capture free-streaming. We show below that this only leads to sub-percent differences compared to the full hierarchy. The CDM+baryon component is treated both as a perfect, pressureless fluid (SPT) as well as an imperfect fluid with corrections to the effective stress tensor (EFT).
We consider three massive neutrino cosmologies with M ν ≡ i m νi = 0.1, 0.2 and 0.4 eV. In all models we assume that the neutrinos are degenerate, each with mass M ν /3. This approximation works well because cosmological observables are very insensitive to the neutrino mass hierarchy except for the total mass, see e.g. [79]. Furthermore, we choose a matching redshift z match = 25, which for the neutrino masses of consideration is well after the non-relativistic transition z nr = 1890 m ν 1 eV , (2.10) in addition to being sufficiently earlier than the point at which non-linearities become important, z 10. The weighted density contrast for the CDM+baryon fluid component reads with f i ≡ Ω i /Ω tot , and similarly for the velocity divergence. We collect the density contrasts and (suitably rescaled) velocity divergences of both fluid components into a vector 2 The equations of motion for the two-component fluid are given by Eq.
Ωm 13) and the non-zero components of the vertex γ abc being (2.14) In the part of the evolution matrix Ω ab describing the neutrino dynamics, we introduced the effective neutrino sound velocity c 2 s,eff . This term captures the free-streaming of the neutrinos c 2 s = δP ν /δρ ν as well as the neutrino anisotropic stress σ ν : where we defined the associated effective free-streaming wavenumber k FS . We compute the terms in the bracket in linear theory, however in a manner that is informed about the complete neutrino distribution function, including all higher moments (see [57]). The equations for the two-component fluid are precisely of the form that can be solved by the extension of SPT discussed above. We expand ψ a in powers of δ 0 = δ cb (η match ) using Eq. (2.5). The accompanying kernels are solutions of Eq. (2.6), with Ω ab and γ abc from Eqs. (2.13) and (2.14). We find that we need to be cautious in choosing initial conditions for the kernel hierarchy at the matching redshift; one easily excites transient solutions that do not entirely decay away by η = 0. Ultimately, we opt for numerically finding the growing mode solution for the kernels with fixed Ω ab (k) = Ω ab (k, η match ) and using this as initial condition. This method is described more comprehensively in [57]. We check that this initial condition matches at the linear level the transfer functions from the Boltzmann solver CLASS [80] for the wavenumbers of interest, and that we obtain the same linear power spectrum at z = 0 (see below).
The total matter power spectrum is the sum of the weighted auto spectra for each component as well as the cross-spectrum, The neutrino energy fraction is f ν = 0.75, 0.15 and 0.3% (constant for z z nr ) for the three models with M ν = 0.1, 0.2 and 0.4 eV, respectively. For wavenumbers larger than the neutrino free-streaming scale k FS ∼ 0.01 h/Mpc, neutrinos are hardly captured in gravitational wells and the neutrino spectra are suppressed compared to P cb,cb . In total, the first term in Eq. (2.16) gives the dominant contribution to the matter power spectrum. Nevertheless, via the backreaction on the gravitational potential, the presence of massive neutrinos leads to a reduction of growth of the CDM+baryon fluid. At the linear level, the relative suppression of the matter power spectrum between a model with massive neutrinos and one without (with Ω cdm tuned so that Ω m is the same) is the well-known −8f ν [9].
In Fig. 2, we compare the two-fluid scheme with the full solution using the Boltzmann hierarchy (CLASS, with maximum neutrino multipole l max = 17) at the linear level. The matter power spectra have sub-permille agreement for k > 0.03 h/Mpc for all neutrino masses. For very low wavenumber, k → 10 −3 h/Mpc there are percent-level deviations. These arise because at the matching redshift z match = 25, the wavenumbers have not yet completely reached the growing mode after entering the horizon, while the two-fluid model initial condition assumes that all scales already reside in the growing mode. In addition, for scales close to Hubble, relativistic corrections become important. Nevertheless, these large scales have negligible impact on loop corrections for a realistic power spectrum.
The neutrino auto power spectrum agrees better than 1% for scales k ∼ 0.01 -0.1 h/Mpc. On smaller scales the difference increase, however in this region errors from the truncation of the Boltzmann hierarchy (choice of l max ) also become significant. We see that the deviation of the neutrino spectrum has a minimal effect on the total matter power spectrum, whose dominant contribution is the CDM+baryon component. Table 1. Approximate time in microseconds to compute one integration point (the integrand for a fixed set of integration variables) at one-and two-loop in the different schemes. Measured on a laptop using one core. 3 Effect of neutrino perturbations on density and velocity spectra

EdS
Having set up the formalism and procedure for computing non-linear corrections taking neutrino perturbations into account beyond the linear level, we proceed in this section to comparing it to simplified treatments. We consider the bare density-and velocity power spectra, and leave a comparison of the renormalized power spectra to Sec. 5. Furthermore, we assess to which extent the effect of free-streaming neutrinos is degenerate with the primordial amplitude of fluctuations at one-and two-loop on BAO scales. The schemes for describing the growth of structure in presence of massive neutrinos that we will consider in this work in increasing order of complexity are: 1. EdS-SPT scheme (1F). The impact of neutrinos is only taken into account at the linear level.
Therefore, the only species modelled non-linearly by the fluid Eqs. (2.3) is CDM+baryons (joint fluid). In addition, the ratio Ω m /f 2 is approximated to 1 so that analytical EdS-SPT kernels may be used. This approach is frequently used in the literature because of its simplicity and efficiency, the latter being particularly favorable for data analysis and parameter scans.
2. One-fluid scheme (1F+). Equivalent to the 1F-scheme, but relaxing the EdS approximation by including the time-dependence of Ω m /f 2 . This scheme is included in order to estimate to what degree the EdS approximation is the source of inaccuracies when comparing 1F and 2F.
3. Two-fluid scheme (2F). Baryons+CDM and neutrinos are both modelled beyond the linear theory in the two-component fluid setup. The time-and scale-dependence of the dynamics both arising from the departure from EdS at late times as well as neutrino free-streaming are captured by solving the kernel hierarchy (2.6) numerically, with the linear evolution and nonlinear vertices given by Eqs. (2.13) and (2.14), respectively.
For the 1F and 1F+ schemes, we use the linear power spectrum at z = 0 (from CLASS) as input to the loop correction calculations. The 2F scheme however, captures also accurately the linear evolution (see discussion above), thus we take the linear power spectrum at z match = 25 as input. The increasing complexity of the different schemes leads to longer processing times; in Table 1 we list the approximate times to compute one integration point on one core on a laptop at one-and two-loop in each scheme. A typical calculation of a loop correction for a single external wavenumber entails O(10 6 ) integral evaluations.

Comparison
To address the accuracy of treating neutrinos only linearly as well as the EdS approximation, we compare the 1F and 2F schemes for the CDM+baryon density and velocity power spectra at z = 0. We consider the power spectrum both at NLO (linear + one-loop) and NNLO (linear + one-loop + two-loop), without EFT corrections (we repeat the comparison including these in Sec. 5).
The comparison is presented in Fig. 3. In the leftmost panel we show the density auto spectrum P δ cb ,δ cb in the 1F and 1F+ schemes, normalized to the 2F result. At one-loop (dashed lines), the 1F scheme deviates at most by 0.5% from the full 2F solution, while at two-loop (solid lines) the deviation exceeds one percent at k 0.17 h Mpc −1 and grows as the two-loop correction becomes increasingly significant. A similar comparison was done for different neutrino masses in [57], and it is reassuring to see consistent results. Since the deviations are roughly independent of neutrino mass  Pθ cb,θcb Figure 3. Comparison of the 1F, 1F+ and 2F schemes for the CDM+baryon density and velocity spectra. We show the power spectrum at one-loop (dashed lines) and two-loop (solid lines) in the 1F scheme, divided by that from the 2F scheme. The dotted lines correspond to the 1F+ scheme normalized to the 2F scheme, at two-loop order. Shaded regions indicate uncertainty from the Monte Carlo integration. Left: density auto spectrum, middle: density-velocity cross spectrum and right: velocity auto spectrum. and there is insignificant difference between the 1F+ (EdS approximation relaxed) and 2F schemes at two-loop (dotted lines), we conclude that the main source of inaccuracy for the density spectrum comes from the EdS approximation. 3 The middle and right panels display the same comparison for the cross spectrum P δ cb ,θ cb and velocity spectrum P θ cb ,θ cb . In contrast to the density spectrum, they show a clear dependence on the neutrino mass in the fractional difference. Moreover, the deviations of the 1F and 1F+ schemes from the full 2F solution are very similar for P θ cb ,θ cb , indicating that the main error comes from neglecting the effect of non-linear neutrino perturbations. It is not too surprising that this effect is larger on the fluid velocity than the density contrast: for scales smaller than the free-streaming scale the growing mode for the CDM+baryon component is approximately (1, 1 − 3/5 f ν ) in the 2F scheme, while it remains (1, 1) in 1F (and 1F+). Therefore, there is an extra suppression for the velocity field compared to the density field, in addition to the change of the growth factor. Nevertheless, the deviation for the smallest neutrino mass M ν = 0.1 eV is less than a percent at two-loop for k 0.15 and 0.2 h Mpc −1 for the cross-and velocity spectra, respectively. We note that for wavenumbers k 0.2 h Mpc −1 the two-loop correction is of the order of the one-loop correction, signalling the breakdown of perturbation theory. Finally, for the velocity spectra it is curious to see that the deviations of 1F to 2F for the one-and two-loop corrections go in opposite directions, thus adding the two-loop seemingly improves the agreement.

Degeneracy with amplitude of fluctuations
The suppression of power on scales smaller than the neutrino free-streaming scale is largely degenerate with a change of the clustering amplitude A s . This may pose a major limitation on constraining the neutrino mass with large-scale structure [81][82][83][84]. The complication is somewhat mitigated by measuring at different redshifts, including higher-order statistics [17,77] and by redshift-space distortions [34]. Moreover, the degeneracy can be broken by including other probes, such as CMB data.
In this light, we check to what extent we can reproduce our results for the 2F matter power spectra in massive neutrino models using massless models by adjusting A s . To isolate the effect of neutrino suppression and overall amplitude, we use the 1F+ scheme for the massless models. Furthermore, we remove completely the double-hard limit of the two-loop corrections, which otherwise yields a large contribution at the scales of interest (after renormalization this contribution is considerably reduced). In Sec. 4 Figure 4. Comparison of the total matter power spectrum in massive neutrino models and massless models with the amplitude of fluctuations changed. We show massless models with A s = (1 − 8fν )As normalized to massive models with amplitude As for the linear (blue), NLO (orange) and NNLO (green) power spectrum.
To isolate the effect of the neutrinos, we compute the massless model in the 1F+ scheme and the massive model in the 2F scheme. Errorbars indicate uncertainty from Monte Carlo integration.
the power spectrum. For simplicity, we set A s = (1 − 8f ν )A s to mimic neutrino free-streaming suppression in the massless models. The results are shown in Fig. 4. For the lowest neutrino mass, the massive and modified massless models agree within a percent. As the characteristic free-streaming wavenumber increases for larger neutrino mass, it becomes more difficult to reconstruct the shape of the massive neutrino non-linear power spectrum with a change of A s . The shape dependence of the linear suppression is more prominent at k 0.1 -0.2h Mpc −1 and due to mode-coupling, the loop-corrections are affected beyond just a rescaling by the presence of the neutrinos.

Effective field theory setup
In order to precisely assess the performance of the simplified treatment of neutrinos (1F) and the full two-fluid solution (2F) versus N-body simulations, we need to take into account EFT corrections. In this section we elaborate on the EFT setup we employ for obtaining a renormalized CDM+baryon power spectrum at next-to-next-to-leading order.
In the effective theory approach, long wavelength fluctuations are modelled perturbatively while systematically taking into account the impact of short wavelengths [22,23]. This is done by coarsegraining the fluid, leading to a modified Euler equation with an effective stress tensor that in addition to the microscopic stress includes corrections from coarse-grained products of short fluctuations. In the notation introduced in Sec. 2, the equations of motion (2.3) in the two-fluid model are modified as is the Kronecker delta, and the effective stress terms are defined equivalently for the CDM+baryon and neutrino fluid in real space as The effective stress tensor τ ij can in the effective theory be split into a deterministic and stochastic contribution; the deterministic part can be modelled in terms of the perturbative solution, while the stochastic part is uncorrelated with the long modes and must be modelled statistically. Due to momentum conservation, the deterministic part must scale as the density contrast times k 2 as k → 0, while the stochastic power spectrum vanishes faster than k 4 in this limit [85,86]. Therefore, we neglect the stochastic contribution in this work.
In the effective theory the deterministic part of the stress tensor is not described from first principles, but written down as a sum of all possible operators allowed by symmetries, with a priori unknown EFT coefficients that must be fitted to simulations or marginalized over in data analysis. Due to Galilean invariance and the equivalence principle, the operators are restricted to gradients and products of the building blocks ∂ i ∂ j Φ and ∂ i v j , where Φ is the rescaled gravitational potential satisfying ∆Φ = δ. At leading order in gradients and number of factors of fields, the effective stress term is [86] τ where d 2 δ and d 2 θ are (time-dependent) EFT coefficients, and δ (1) and θ (1) are the linear solutions in perturbation theory. In massless neutrino cosmologies, they are related by δ (1) = −θ (1) /Hf , so that the sum above can be reduced to Hf . This is not the case for massive neutrino models, due to free-streaming, therefore we in principle need to write down in the two-fluid model with EFT coefficients d cb,i , d ν,i . This leads to a counterterm CDM+baryon density fieldψ in Fourier space. Each coefficient d i depends on a linear combination of both sets of coefficients d cb,i , d ν,i because the linear propagation after the insertion of the sources τ θ cb and τ θν mixes the contributions from the four components.
In practice, we are interested in weakly non-linear scales k 0.1h Mpc −1 , and assume a large scale separation k FS k k NL , where the non-linear scale is k NL ∼ 0.3h Mpc −1 . As we return to below, this assumption is justified for M ν = 0.1 and 0.2 eV but begins to break down for M ν = 0.4 eV. In particular, the linear CDM+baryon density contrast and (rescaled) velocity divergence can be related by a constant factor far below the free-streaming scale: we have approximately ψ 4 -terms. Note also that the sound velocity in the microscopic part of the τ ij ν stress tensor is already captured in the 2F model, as defined in Eq. (2.15). In total, the leading CDM+baryon counterterm density can be written as where γ 1 is an EFT parameter and we extracted the overall growth e 2∆η for convenience. Since we will determine γ 1 via calibration to N-body simulations, its explicit relation to the coefficients d i and kernels F i (η) is not important.

Renormalization of the one-loop correction
To assess whether the counterterms from the effective stress term in the EFT can cure the cutoffdependence of SPT (also in the 2F scheme), we study its UV sensitivity, starting with the one-loop correction. In the EdS approximation (1F scheme), the leading contribution (i.e. at leading power in the gradient expansion k 2 /k 2 NL ) to the SPT one-loop power spectrum for hard loop momentum q → ∞ is where is the displacement dispersion. P 0 is the input linear power spectrum from the initial condition, as defined above. From here on we will only consider the CDM+baryon power spectrum, and we drop the cb subscript for brevity. Due to momentum conservation, the kernels F (n) (q 1 , . . . , q n ) scale as k 2 when the sum k = i q i goes to zero. Therefore, the dominant contribution to the hard limit above comes from the diagram containing F (3) , and the coefficient c is given by We expect the k 2 -scaling guaranteed by momentum conservation also in the 2F scheme, but the c-coefficient can change and acquire a time-dependence inherited from the kernels. To analyse the UV and determine c(η) in the 2F model, we compute the hard limit of the one-loop integral numerically. This can be done as follows: Define the one-loop integrand p 1L by P 1L (k, η) = e 4∆η Λ dq q 2 P 0 (q) dΩ q p 1L (k, q; η) . (4.10) It contains the different diagrams contributing to the one-loop correction with the corresponding kernels and an additional linear power spectrum. In the hard limit, the leading contribution to p 1L should scale as F (3) (k, q, −q) ∼ k 2 /q 2 , therefore we can factorize the integration: To evaluate p h 1L , we fix the loop momentum to a large value q Λ and integrate over the angle numerically. Comparing to Eq. (4.7), we find . (4.12) In Fig. 6 we show the results for the c-coefficient at η = 0 for the different neutrino masses, using Λ = 1h Mpc −1 and fixing q = 10h Mpc −1 . The constant c = 61/210 in the 1F scheme is shown for comparison. Furthermore, we display dashed lines where the linear kernel is divided out in order to isolate the scaling of the F (3) kernel in the hard limit, i.e. c/F  the limit k k FS , neutrinos behave as dark matter, and the system can be mapped onto a set of standard dark matter fluid equations for which the hard limit yields c constant. We do not recover the 1F limit of 61/210, because the EdS approximation is relaxed and because the kernel F (3) (k, q, −q; η) is sensitive to the dynamics at q k FS where the growth rate is dampened due to neutrino freestreaming. For the largest neutrino mass with the largest free-streaming wavenumber, this limit is visible around k = 10 −3 h Mpc −1 in Fig. 6.
On scales much smaller than the free-streaming scale, neutrinos do not cluster and the CDM+baryon and neutrino fluid system effectively decouple. The CDM+baryon fluid is then equivalent to the standard case (1F), with the EdS approximation relaxed and with reduced growth due to the effect of the neutrinos on the background evolution. Hence c approaches a constant for k k FS . Our aim is to obtain a renormalized power spectrum in the mildly non-linear regime, k 0.1h Mpc −1 , therefore we assume a large scale separation k FS k q 5 . In this limit c is constant, which implies that the counterterm (4.6) can absorb the cutoff dependence of the one-loop correction, as we show explicitly below. For the lowest neutrino masses, we see indeed a plateau in the value of c around k = 0.1h Mpc −1 in Fig. 6. The measured value is read off at the dashed line and quoted in Tab. 2. For the largest neutrino mass, the assumption of large scale separation k FS k starts to break down (from Eq. (2.15) we have e.g. k FS (k = 0.1h Mpc −1 , η = 0) = 0.07h Mpc −1 ), and the plateau emerges only on even smaller scales. Nevertheless, we read off the value of c and include the model in our analysis for illustration. In practice, the limit k FS k q is well satisfied for relevant values of the sum of neutrino masses.
The counterterm (4.6) gives the additional contribution to the power spectrum, which is exactly of the form to absorb the UV sensitivity (4.7). The renormalized one-loop CDM+baryon power spectrum therefore reads P ren 1L (k, η, Λ) = P 1L (k, η; Λ) + P ctr 1L (k, η; Λ) . (4.14) Demanding that the renormalized result is cutoff-independent yields the renormalization group equa- at leading power in gradients k 2 /k 2 NL , which has the solution where γ 1 is an initial condition.

Renormalization of the two-loop correction
Next we discuss the UV-sensitivity of the two-loop correction and how we can treat it in the EFT. The second order contributions to the effective stress term are [87,88] where the number superscript indicates the perturbative order and the tidal tensor is defined as Three additional EFT parameters e 1 , e 2 and e 3 appear at this order. In writing down Eq. (4.17) we have used again that in massless neutrino cosmologies (1F) θ (1) = −Hf δ (1) . In addition, the operator ∆θ (2) can be written as a linear combination of the existing operators, hence it is also redundant. These simplifications do not hold in general for massive neutrino cosmologies (2F). While the linear CDM+baryon perturbations could be related far below the free-streaming scale, ∆θ (2) is not redundant in general. Furthermore, there are additional operators including the neutrino perturbations at second order as well as mixed operators with linear CDM+baryon and neutrino perturbations. Given the smallness of the neutrino perturbations around the weakly non-linear scales those should be negligible however. Even in the massless case, solving the perturbative equations of motion in the presence of the second order counterterms (4.17) is quite complex. Moreover, in practice, there are significant degeneracies between the EFT parameters when calibrating to simulations. Therefore, making an ansatz for a linear relation between the parameters and fitting the overall amplitude is typically sufficient [71]. One also easily overfits the data when introducing too many free parameters (we return to this issue in Sec. 5). All in all, we follow the prescription of [71] for renormalizing the two-loop correction and extend it to the 2F scheme. This works well for the power spectrum and an analogous method has also been applied successfully to the bispectrum at two-loop [39]. The upshot is that we do not need to know the specific form of the second order effective stress term, also in the massive neutrino case.
Note that for renormalizing the two-loop, one in principle also needs counterterms from the third order effective stress term. However, for the power spectrum, they have been argued to yield contributions that are degenerate with those arising at lower order [89]. The explicit third order correction to the stress tensor is not needed within the approach followed here.
To study the UV limit of the two-loop correction, it is useful to distinguish between the doublehard (hh) limit where both loop momenta are hard compared to the external momentum k, and the single-hard (h) limit where one of the momenta are hard compared to the external and the other loop momentum. We discuss first the double-hard limit.

Double-hard limit
In the 1F case, the leading contribution to the double-hard region of the two-loop comes from the propagator correction diagram containing F (5) [71]. This follows from the k 2 -scaling of the kernels in the limit k/q 1 ∝ k/q 2 1, and also holds in the 2F case guaranteed by momentum conservation. The other diagrams enter in this limit with additional powers of k 2 and can be neglected at leading order in gradients. The F (5) -diagram contribution scales as k 2 P 0 (k) and is therefore degenerate with the UV-limit at one-loop and can be absorbed by the counterterm (4.13). We can split the corresponding EFT parameter at NNLO into a one-and two-loop part, Following [71], we use a renormalization scheme where we choose γ 2L 1 so that it exactly cancels the double-hard region of the two-loop, while γ 1L 1 is calibrated to simulations and contains the finite part of the counterterm. We drop the superscript and use γ 1 ≡ γ 1L 1 in the numerical analysis below. This means in practice that we subtract the double-hard contribution from the two-loop correction: The double-hard contribution P hh 2L can be computed analytically with EdS-kernels in the 1F scheme (see e.g. [71]), and numerically in the massive neutrino models by fixing the loop momenta to large values. Alternatively, one can fit the low-k region by computing b hh = lim k→0 P 2L /k 2 /P 0 and subtracting P hh 2L = b hh k 2 P 0 . The double-hard limit is realized in this region because the loop integral has sole support from wavenumbers q 1 , q 2 k. We use the latter method. In the left panel of Fig. 7 we show the ratioP 2L /k 2 /P 0 for the 1F scheme and for the different neutrino masses using the 2F scheme (more precisely, we show the corresponding IR resummed quantities defined in Eqs. (4.47) and (4.48) below). The results are obtained using a cutoff Λ = 1h Mpc −1 . In the 1F scheme, we see that the two-loop attain the k 2 P 0 (k)-scaling below k 10 −2 h Mpc −1 . Due to the presence of the free-streaming scale in the 2F neutrino models, they only recover the asymptotic k 2 P (k)-scaling for even smaller k. We read off the limit b hh at k = 10 −3 h Mpc −1 (dashed vertical line in Fig. 7). Note that b hh defines the renormalization point, and changing its value only leads to a shift in the EFT parameter γ 1 . We check that adjusting the wavenumber at which we read off the limit b hh has insignificant impact on our conclusions, up to the shift. Therefore we find it convenient to evaluate the limit for k k FS , while we will obtain the corresponding limit b h for the single-hard two-loop contribution at k k FS as we discuss below. 6

Single-hard limit
Next, we consider the limit in which one of the loop momenta becomes large: q 1 k, q 2 ∼ k (or equivalently with q 1 , q 2 switched). The counterterms that renormalize this limit are one-loop diagrams with insertions of EFT operators. Computing those diagrams is rather complicated even in the 1F case, and would in the 2F scheme potentially involve further operators from the generalization of (4.17) including neutrino perturbations. Therefore we opt for the approach put forward in [71], using an ansatz for the single-hard counterterm with one extra EFT parameter.
Analogously to the one-loop case, we start by defining the two-loop integrand p 2L : In the single-hard limit q 1 → ∞, the leading contribution comes from kernels F (n) (. . . , q 1 , −q 1 , . . . ) ∝ 1/q 2 1 . The two-loop integral can thus be factorized in this limit giving the single-hard result The factor 2 in (4.22) accounts for the equivalent contribution where q 2 is hard. We compute the limit (4.23) both in the 1F and 2F schemes numerically by fixing the hard loop momentum to a large value q 1 Λ and evaluating the resulting integral using Monte Carlo integration.  Following the prescription of [71], we assume that we can correct for the spurious UV contribution in the single-hard limit by a shift in the value of the displacement dispersion σ 2 d , i.e.
where the γ 2 -term corresponds to the EFT correction and N is a factor introduced for convenience that we define below. The replacement yields an additional counterterm (4.25) to the power spectrum. This term could in principle immediately be added to the renormalized power spectrum, however we notice that part of it is still degenerate with the double-hard limit. The definition (4.23) integrates over a hard region q 2 k for external wavenumbers far below the cutoff, which yields a contribution proportional to k 2 P 0 (k). Such a contribution is absorbed by the γ 1counterterm defined above, however as discussed for the double-hard limit, we adopt a renormalization scheme where γ 2L 1 exactly cancels the double-hard region of the two-loop, now including also the double-hard contribution from the counterterm (4.25). In practice this means that we subtract the degenerate double-hard part from the single-hard contribution in analogy to the subtracted two-loop correction: To obtain the double-hard limit p hh 2L , we use the same method as for the full two-loop correction: we evaluate numerically the limit b h = lim k→0 p h 2L /k 2 /P 0 (k). In this limit the integral (4.23) only has support where q 2 k and since q 1 is hard we recovered the double-hard limit and p hh 2L = b h k 2 P 0 (k). The subtracted single-hard limitp h 2L divided by the double-hard scaling k 2 P 0 (k) is shown in the right panel of Fig. 7 (we plot the corresponding IR resummed quantities defined in the next section). We set the cutoff to Λ = 1h Mpc −1 and the hard loop momentum to q 1 = 10h Mpc −1 .
In the 1F scheme, the low-k limit has indeed the k 2 P 0 scaling, and we subtract the value at k = 10 −2 h Mpc −1 so that the single-hard limitp h 2L does not give a contribution degenerate with the one-loop counterterm. For the 2F scheme we also find the k 2 P 0 scaling is approached, however for very small wavenumbers certain deviations occur. They may in part be attributed to large numerical cancellations due to our choice q 1 = 10h Mpc −1 , and we checked that the k 2 P 0 scaling extends to lower k for q 1 = 5h Mpc −1 , without changing the result within the relevant range k 0.05h Mpc −1 wherep h 2L gives a significant contribution to the renormalized two-loop power spectrum. We therefore use a subtraction point k 0.015, 0.027, 0.04h Mpc −1 for M ν = 0.1, 0.2, 0.4 eV, to read off b h and subtract the double-hard contribution (indicated by the vertical dashed lines in Fig. 7). As for the full two-loop correction, the definition of the double-hard subtraction point and corresponding value b h determines the renormalization point, hence a change of b h only results in a shift of γ 1 . We check in particular that choosing 0.027h Mpc −1 as subtraction points also for M ν = 0.1 and 0.4 eV has negligible impacts on the results.
All in all, we have the two-loop counterterm where P ren NNLO (k, η; Λ) = P tree (k, η) + P ren 1L (k, η; Λ) + P ren 2L (k, η; Λ) . By removing the double-hard contribution to the two-loop correction (i.e. by using the subtracted quantitiesP 2L andp h 2L ), there are no degenerate UV contributions between the one-and two-loop corrections, and we can solve the RGE (4.29) independently for each term. Furthermore, the one-loop RGE (4.15) and its solution stay the same.
The remaining cutoff-dependence comes from the single-hard region of the two-loop correction, which should be corrected for by the two-loop counterterm, valid at leading power in gradients. This result reflects the assumption that the single-hard region is regulated by a shift in the displacement dispersion σ 2 d . The solution is where γ 2 is an initial condition. 7 We choose N = 1/c, with c defined in Eq. (4.9) (and evaluated in the 2F scheme using Eq. (4.12)), in order to treat γ 1 and γ 2 on equal footing (cf. RGE for γ 1 (4.15)).
To check that we indeed obtain cutoff-independent loop corrections, we show the various loop and counterterm results in Fig. 8, for M ν = 0.1 eV using two cutoffs Λ = 0.8 (dashes lines) and 1h Mpc −1 (solid lines). The top panels show the absolute contributions, while the bottom panels show the fractional difference between the renormalized power spectra using the two cutoffs. We see that the cutoff-dependence of the renormalized one-loop correction is on the permille level up to k 0.4h Mpc −1 ; for increasing k the neglected next-to-leading gradient corrections O(k 4 /k 4 NL ) become more and more important. The double-hard limit gives the largest contribution to the bare two-loop correction, as indicated by the bare (blue) and subtracted (yellow) graphs in the right panel  of Fig. 8. The single-hard limit is regulated by the counterterm (green) yielding the renormalized two-loop correction (red). We display its fractional difference using the two cutoffs in the bottom plot, where the large relative deviations at small k come from the renormalized correction being very close to zero. The difference is less than 1% up to about k = 0.25h Mpc −1 .

IR resummation
A well-known issue with the Eulerian perturbative description is that the effect of large bulk flows on the BAO peak is inadequately modelled [25,90]. Nevertheless, this effect can relatively straightforwardly be taken into account non-perturbatively by resumming the contributions from long-wavelength displacements, known as IR resummation [27][28][29][30]. We use IR resummation for both schemes 1F and 2F in order to model the BAO scales as accurately as possible.
The large bulk flows affect only the BAO wiggles and not the broadband part of the power spectrum, therefore it is useful to split the linear power spectrum into a wiggly and non-wiggly part P 0 (k) = P nw (k) + P w (k) . (4.33) We perform this split using the method described in [91,92], i.e. we obtain P nw by Fourier transforming the power spectrum to real space, removing the BAO peak, and smoothly interpolating the power spectrum without the peak. Following [30], the next step is to compute the damping factor where j n are spherical Bessel functions, k osc = h/(110 Mpc) is the wavenumber that corresponds to the BAO period and k s is the separation scale of long and short modes, introduced to have separate perturbative expansions in the two regimes. We use k s = 0.2h Mpc −1 , as recommended by [30]. 8 The IR resummed power spectrum at N -th order is where the bracket notation P L-loop [X] indicates that X is the linear input power spectrum that enters in the L-loop integrals. The exponential damping factor resums the infrared contributions, while the sum over n corrects for overcounting (the IR contributions are also included perturbatively order by order in the loop evaluations) [30]. 9 In the last equality we defined a shorthand notation for the input power spectrum in the IR resummed result P IR,N −L 0 , i.e. the quantity in the bracket, which depends on the difference N − L. Specifically, at the orders in perturbation theory in which we work one has P IR LO = P tree P nw (k) + e −k 2 Σ 2 P w (k) , (4.36) Note that the tree-level power spectrum is non-trivial in the 2F scheme and given by In deriving Eq. (4.35) in the Eulerian theory, the following property of the kernels in the soft limit is key [30]: · · · q · q n q 2 n (4. 40) in the limit where q 1 , . . . , q m are much smaller than the remaining arguments q m+1 , . . . , q n and where q = q i . This factorization follows from Galilean invariance [93], which is certainly also a symmetry in the 2F scheme. Therefore the IR resummed equation (4.35) is valid also in the 2F case.
In practice, we use the following simplified expression to evaluate IR resummed loops which is valid up to corrections of O(P 2 w ) as well as diagrams involving P w inside a hard loop. Neglecting these corrections is a good approximation because P w is small compared to P nw , and it oscillates around zero so that its integral vanishes. 8 We checked that changing the separation scale to ks = 0.1h Mpc −1 has negligible impacts on our results and conclusions. 9 We have lim N →∞ P N −L 0 = Pnw + Pw = P 0 such that the exact spectra with and without IR resummation agree; the difference arises only when working at finite order in perturbation theory.

All together: summary
We end this section by writing down the final expressions for the power spectrum up to NNLO including both IR resummation and EFT corrections. We have where the counterterms are given by with c defined in Eq. (4.9) and the subtracted quantities defined as The single-hard contribution p h 2L was defined in Eq. (4.23) and the RGEs for the EFT parameters were given in Eqs. For completeness, we note again that in the 1F scheme the time-dependence is factored out, so to obtain the power spectrum at z = 0 we can use η ini = 0 and ∆η = 0. In the 2F scheme on the other hand, we chose the matching between the Boltzmann hierarchy and the two-fluid description at η ini corresponding to z = 25. Therefore, to obtain the quantities P IR,N −L 0 , we take the linear power spectrum from CLASS at z = 0 or z = 25 depending on the scheme, split it into its wiggly and non-wiggly components, compute the damping factor Σ 2 and use the definition in Eq. (4.35).

Comparison to N-body simulations
To calibrate the EFT parameters and perform a comparison of the 1F and 2F schemes in the effective theory, we utilise the Quijote N-body simulations [72]. In particular, we use the set of massive neutrinos simulations with 512 3 CDM particles and 512 3 neutrino particles in a cubic box of size (1 Gpc/h) 3 . Moreover, we use the simulations with pair fixed initial conditions that significantly reduce the cosmic variance for the power spectrum. The cosmological parameters are Ω m = 0.3175, Ω b = 0.049, h = 0.6711, n s = 0.9624, σ 8 = 0.834 and the neutrinos mass is M ν = 0.1, 0.2 and 0.4 eV in the different models respectively. We estimate the power spectrum and uncertainty at z = 0 using the average and standard deviation of the 500 realizations.
On the largest scales, the binning method used to obtain the power spectrum from the simulations become significant and has to be taken into account when comparing simulation and theory. We use the Pylians library 10 to perform the same binning on the linear power spectrum. Instead of binning also the loop corrections, we can equivalently compare the theoretical results to the "unbinned" Nbody power spectrum The EFT parameters are calibrated by minimizing the following χ 2 at NLO and NNLO: where ∆P data (k) is the estimated uncertainty of the N-body power spectrum. We have k min = 0.0089h Mpc −1 and at k max = 0.1 and 0. , is not fitted sinceγ 1 is fixed from the two-parameter NNLO case. We explain the motivation for this case below. The third case, NNLOγ 2 =γ 1 , is motivated by the assumption that in the EFT, the displacement dispersion is corrected in an universal manner as σ 2 d → σ 2 d + γ 1 /c both at one-and two-loop.

Results
We begin by comparing the full 2F solution for the various NLO/NNLO calibration cases above with the most phenomenologically relevant neutrino mass M ν = 0.1 eV. In Fig. 9 we show the power spectrum using three pivot scales k max = 0.103, 0.148 and 0.204h Mpc −1 , normalized to the N-body results. The gray band indicates estimated uncertainty from the N-body simulation and the red shading corresponds to an ansatz for the expected theoretical uncertainty at NLO (light shading) and NNLO (darker shading) [94]. In all panels, there is a "bump" feature in the relative difference around k = 0.01 -0.05h Mpc −1 of a few permille (i.e. well within the uncertainty of the N-body data) that arises due to finite N-body box size and binning uncertainty not captured by the correction (5.1). Nevertheless, all cases can be fitted to lie within the N-body uncertainty up to even the largest pivot scale. Depending on the loop order, we are subject to overfitting above a certain k max . In particular, given the expected theoretical uncertainties, we see that the NLO result differs much less from N-body beyond k 0.10h Mpc −1 than expected from the missing two-loop correction, and similarly the NNLO result remains in good agreement above k 0.10h Mpc −1 , indicating overfitting. Ideally, we would be able to calibrate the EFT parameters precisely below k = 0.05h Mpc −1 where the perturbative uncertainty is small, but in practice the moderate simulation uncertainty (even after cosmic variance is reduced), a small number of grid points and finite box size effects complicates the picture.
To gauge the extent of overfitting at NLO, we introduce a 0-parameter NLO case where the γ 1 parameter is fixed to that in the NNLO case with the same pivot scale. At NNLO, there is a larger window in which the theoretical uncertainty is still small and hence more data points to measure γ 1 precisely. Moreover, in our renormalization scheme γ 1 remains the same after adding the two-loop correction, up to numerical uncertainty and overfitting. Therefore, the degree of overfitting is reduced in the NLOγ 1 =γ  The χ 2 per degree of freedom (d.o.f.) is shown for the different fits in Fig. 10. Since we are ultimately interested in assessing the effect of non-linear neutrino perturbations on the power spectrum, i.e. comparing the performance of the 1F and 2F schemes, we chose a simple estimator for the N-body uncertainty, for example not taking into account correlations between different bins. This is reflected by the small χ 2 /d.o.f.
1. We find that theγ 1 =γ 2 ansatz at NNLO proposed by [71] does not work as well when including the exact time-dependence and effect of neutrinos, and even performs worse than the NLO case at k 0.12h Mpc −1 .
The best-fit EFT parametersγ 1 andγ 2 as a function of k max are displayed in Fig. 11. We show the results from both schemes 1F and 2F, indicating 1σ uncertainty in the shaded bands. At NLO, we find that theγ 1 measurement is consistent with a constant up to about k max 0.11h Mpc −1 where the parameter starts exhibiting a running. At NNLO (with both parameters varied), this plateau for γ 1 extends to k max 0.14h Mpc −1 , and the measurement within 1σ is consistent with a constant up to around 0.2h Mpc −1 . Theγ 2 parameter (rightmost panel) can only be constrained beyond k max 0.15h Mpc −1 , because the two-loop counterterm (i.e. the subtracted single-hard limit (4.27)) is exceedingly small for lower wavenumbers. We find that when both parameters are allowed to vary  Figure 11. Measured parameters γ1, γ2 and γ = γ1 = γ2 as a function of kmax. We show the results from the 1F scheme (blue) and 2F scheme (green) for a neutrino mass sum Mν = 0.1 eV. The shaded regions indicate 1σ uncertainty.  at NNLO, they assume different values, giving another confirmation that the NNLO 1-parameter γ 1 =γ 2 model cannot match the data with the same accuracy. In the second leftmost plot we show the best-fit parameter in this model. It features a considerable running compared to the other cases even at low k max 0.1h Mpc −1 , suggesting a degree of overfitting already from large scales. Before commenting on the difference in the calibrated EFT parameters in the 1F and 2F schemes as is seen in Fig. 11, we discuss the power spectrum and χ 2 in the two schemes.
Next, we compare the 1F and 2F schemes for M ν = 0.1 eV. We find similar conclusions for the other neutrino mass cosmologies. The power spectra for the two schemes are shown in Fig. 12, using k max = 0.15h Mpc −1 , and the χ 2 /d.o.f. as a function of the pivot scale is shown in the left panel of Fig. 13. We find that both schemes can broadly model the data equally accurately, with some minimal differences. This suggests that the difference between the schemes for the bare density power spectrum (arising mostly due to relaxing the EdS approximation, see Fig. 3 and discussion in Sec. 3) can to large extent be absorbed into the counterterms. In particular, the difference can be captured by the k 2 P 0 -counterterm leading to a shift in the valueγ 1 . This is exactly what we see in Fig. 11. Both in the NLO and NNLO cases, we have an approximate shift ∆γ 1 ≈ −0.2 Mpc 2 /h 2 between the 1F and 2F schemes. The shift is consistent with the findings of [39] for the one-loop bispectrum with the EdS approximation relaxed. Furthermore, we find that the two-loop EFT parameterγ 2 is unchanged in the 2F scheme as compared to 1F, as seen for k max > 0.15h Mpc −1 in the rightmost panel of Fig. 11.
Finally, we consider the three cosmologies with M ν = 0.  (2-parameter) in the full 2F solution. The reduced χ 2 for the different models is displayed in the right panel of Fig. 13. By accident, the "bump" feature from finite box effects is slightly smaller for M ν = 0.2 eV than the other neutrino masses, leading to a smaller χ 2 for low k max . Moreover, for M ν = 0.4 eV, the linear 2F solution is off by almost a percent for k < 0.01h Mpc −1 (see Fig. 2), yielding a slightly larger χ 2 around k min . Apart from these insignificant discrepancies, we find that all neutrino models can match the N-body results accurately using the 2F scheme. This is also seen from the power spectra plots in Fig. 1.
The best-fit EFT parameters for the different neutrino mass models are shown in Fig. 14. In the left plot we use the 1-parameter NLO approach while in the middle and right plots we show the NNLO 2-parameter model. It is not surprising that the parameters acquire different values for the different cosmologies: the part of the counterterms correcting for the inaccuracy of SPT is altered due to different values of the displacement dispersion σ 2 d and the c-coefficient, and the part accounting for the impact of actual short-scale physics changes due to the impact of massive neutrinos on strongly non-linear scales. At NNLO, we measure a shift ∆γ 1 −0.14 Mpc 2 /h 2 at k max = 0.14h Mpc −1 and ∆γ 2 −0.2 Mpc 2 /h 2 at k max = 0.17h Mpc −1 when doubling the sum of neutrino masses. We quote the value of the c-coefficient as well as the measured EFT parameters γ 1 (Λ) and γ 2 (Λ) for Λ = 1h Mpc −1 and Λ → ∞ for the different neutrino cosmologies in Tab. 2. In summary, we can accurately match the Quijote N-body data both at NLO and NNLO using the full 2F solution, however being subject to a certain degree of overfitting. The window where the EFT parameters can be precisely measured is relatively small and involve uncertainties from the N-body simulation. We compared the 1F and 2F schemes and found that the difference of the unrenormalized power spectrum can largely be absorbed in the one-loop counterterm, leading to a shift in the EFT parameterγ 1 .

Conclusions
In this work we have performed a precision calculation of the CDM+baryon density and velocity power spectrum at next-to-next-to-leading (two-loop) order taking the full time-and scale-dependence of neutrino perturbations beyond the linear level as well as exact ΛCDM (+M ν ) time-dependence into account. Our aim is to assess the accuracy of the commonly used approach in which the aforementioned effects are neglected.
To describe the impact of neutrinos in gravitational collapse, we use a hybrid two-component fluid scheme introduced in [61] and further developed in [57]. The idea is that neutrinos can only be described by a fluid well after they have become non-relativistic, because the lower and higher order multipoles effectively decouple by powers of T ν /m ν . Therefore, we can use the full, linear Boltzmann hierarchy until an intermediate redshift z = 25, where we map the equations onto a two-component fluid consisting of CDM+baryons (one joint component) and neutrinos. The large neutrino free-streaming is reflected in an effective neutrino sound velocity, which we compute from linear theory. The two-component fluid system can readily be realized in the framework for computing loop corrections in cosmologies with general time-and scale-dependence introduced in [57].
We compare the (unrenormalized) two-loop CDM+baryon density and velocity power spectra at z = 0 from the full solution including neutrino perturbations beyond the linear level and exact time-dependence (2F scheme) to the commonly used simplified treatment with only linear neutrino perturbations and EdS dynamics for CDM+baryons (1F scheme). For the density spectrum, the main difference arises due to the departure from EdS, leading to 0.7% deviation at k = 0.15h Mpc −1 (consistent with the results in [57]), approximately neutrino mass independent. On the other hand, for the velocity spectrum, the deviation in the 1F scheme is neutrino mass dependent; for the largest neutrino mass M ν = 0.4 eV the deviation is larger than a percent for k 0.08h Mpc −1 . For k = 0.15h Mpc −1 the deviation in the velocity or cross power spectra up to two-loop are 2.7%, 1.3%, 0.8% for M ν = 0.4, 0.2, 0.1 eV, respectively, and even somewhat larger, 2.7%, 2.0%, 1.6%, when going only up to one-loop, due to a partial cancellation between one-and two-loop terms. Therefore, the commonly used 1F approximation is questionable for the velocity spectra when aiming at percent accuracy, implying that the 2F scheme should be used to access the neutrino mass sensitivity in redshift space.
Furthermore, we check to what extent the two-loop power spectrum on mildly non-linear scales in the presence of free-streaming massive neutrinos can be reproduced by a massless model by adjusting the overall amplitude of fluctuations. We find that the massless model can reproduce the total matter density power spectrum of the smallest neutrino mass to within a percent, but for the larger neutrino masses the neutrino free-streaming scale k FS is larger and mimicking the power spectrum in the massless model is more challenging.
To assess the performance of the full 2F solution and the simpler 1F scheme as accurately as possible, and to compare to N-body data, we take into account EFT corrections as well as effects of large bulk flows from the IR. We find that the k 2 scaling of the kernels -expected from momentum conservation in the limit when the external wavenumber k goes to zero (or the other momenta go to infinity) -is slightly spoiled by the free-streaming scale, but recovered when k k FS or k k FS . Therefore, we assume a large scale separation k FS k q, where q is the scale of hard loop momenta, when analysing the UV limit of the 2F model. This assumption is appropriate for the smallest neutrino masses, but breaks down for the largest neutrino mass. At one-loop, we show that the usual k 2 P 0 (k) counterterm absorbs the cutoff-dependence. The double-hard region of the two-loop contribution can also be corrected by this counterterm. For the single-hard region, we employ a numerical treatment yielding a counterterm that exactly captures the scale-dependence of the UV-limit of SPT, with one extra EFT parameter.
We compare the perturbative results of both the 1F and 2F schemes to N-body data, at different orders in perturbation theory and with various number of free parameters. The power spectrum in the 2F scheme matches the N-body data within a percent up to k 0.12h Mpc −1 and k 0.16h Mpc −1 for NLO and NNLO, respectively, using a pivot scale k max = 0.148h Mpc −1 . We find that the 1F scheme can achieve similar accuracy, indicating that the deviation in the bare power spectrum can be absorbed by the counterterms. In particular, the main source of error in the 1F density power spectrum comes from the EdS approximation, and it is largely degenerate with the one-loop counterterm. We measure a shift in the EFT parameter ∆γ 1 = −0.2 Mpc 2 /h 2 between the 1F and 2F schemes that accounts for the discrepancy.
In total, we have scrutinized the effect of massive neutrino perturbations taking into account the full impact of time-and scale-dependent growth on non-linear kernels for CDM+baryon as well as neutrino density and velocity. We further capture the exact time-dependence of ΛCDM (+M ν ). The impact on the two-loop density power spectrum is to great extent degenerate with counterterms. The influence of time-and scale-dependent growth due to massive neutrinos is larger on the velocity spectrum, suggesting that the full 2F treatment for neutrinos is warranted to access the neutrino mass information encoded in redshift space distortions. This motivates an analysis of the power spectrum in redshift space within the 2F scheme, which we leave to future work.