Accurate fitting functions for peculiar velocity spectra in standard and massiveneutrino cosmologies
^{1}
AixMarseille Univ., Université de Toulon, CNRS, CPT, Marseille, France
email: jbel@cpt.univmrs.fr
^{2}
INAF – Osservatorio Astronomico di Brera, via E. Bianchi 46, 23807 Merate, Italy
^{3}
Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Magrans s/n, 08193 Barcelona, Spain
^{4}
Institut d’Estudis Espacials de Catalunya (IEEC), 08034 Barcelona, Spain
^{5}
Università degli Studi di Milano, via G. Celoria 16, 20133 Milano, Italy
^{6}
INAF – Osservatorio Astronomico di Trieste, Via Tiepolo 11, 34143 Trieste, Italy
^{7}
INFN – Sezione di Trieste, Via Valerio 2, 34127 Trieste, Italy
^{8}
INFN – Sezione di Milano, via G. Celoria 16, 20133 Milano, Italy
Received:
26
October
2018
Accepted:
29
November
2018
We estimate the velocity field in a large set of Nbody simulations including massive neutrino particles, and measure the autopower spectrum of the velocity divergence field as well as the crosspower spectrum between the cold dark matter density and the velocity divergence. We perform these measurements at four different redshifts and within four different cosmological scenarios, covering a wide range in neutrino masses. We find that the nonlinear correction to the velocity power spectra largely depends on the degree of nonlinear evolution with no specific dependence on the value of neutrino mass. We provide a fitting formula based on the value of the rms of the matter fluctuations in spheres of 8 h^{−1} Mpc, describing the nonlinear corrections with 3% accuracy on scales below k = 0.7 h Mpc^{−1}.
Key words: largescale structure of Universe / dark matter
© ESO 2019
1. Introduction
The analysis of the largescale structure of the universe provides crucial information on the evolution of the background matter density and its perturbations (see, e.g., Bernardeau et al. 2002; Bassett et al. 2010; Weinberg et al. 2013). In particular, analyses of cosmological probes sensitive to different cosmic epochs could lead to an explanation of the mysterious latetime acceleration of cosmic expansion (e.g., Amendola et al. 2018). At the same time, largescale structure is sensitive to the details of the standard model of particle physics, providing upper bounds on the neutrino mass scale (see, e.g., Lesgourgues & Pastor 2006).
Measured redshifts, which are used to estimate galaxy distances, are affected by galaxy peculiar velocities generated by the growth of cosmological matter perturbations. Their lineofsight component combines with the cosmological expansion, systematically modifying the derived galaxy distances and generating what are known as redshift space distortions (RSD). This effect turns the amplitude and isotropy of redshiftspace clustering statistics into a sensitive probe of the linear growth rate of structure f = dlnD/dlna (Kaiser 1987).
Combining measurements of the expansion history H(z) and the growth rate of structure f can evidence deviations from the standard theory of gravity, that is, General Relativity. This (Guzzo et al. 2008) has led to renewed interest in RSD over the past decade (see e.g., Sánchez et al. 2017; Pezzotta et al. 2017 for the most recent analyses and a summary of previous results).
Extracting the linear growth rate from RSD measurements of biased, nonlinear tracers as galaxies however requires an accurate modeling of redshiftspace galaxy clustering. Extensive work over the past decade has addressed this in the context of cosmological perturbation theory, both in the Eulerian and Lagrangian formulations, providing linear and nonlinear predictions for the redshiftspace galaxy power spectrum (see, e.g., Kaiser 1987; Scoccimarro 2004; Matsubara 2008a,b; Taruya et al. 2009, 2010, 2013; Reid & White 2011; Sato & Matsubara 2011; Seljak & McDonald 2011; Valageas 2011; Zhang et al. 2013; Zheng et al. 2013; Senatore & Zaldarriaga 2014; Uhlemann et al. 2015; Okumura et al. 2015; Perko et al. 2016; Bose & Koyama 2016; Bianchi et al. 2016; Vlah et al. 2016; Hand et al. 2017; Hashimoto et al. 2017; Fonseca de la Bella et al. 2017). Upadhye et al. (2016) in particular investigate RSD, taking into account massive neutrinos and dark energy.
In the formulation of Scoccimarro (2004), and several that followed, the largescale effect of RSD is described in terms of two main ingredients: the autopower spectrum of the velocity divergence field P_{θθ}, and the crossspectrum between the velocity divergence and the matter density contrast P_{δθ}. Although, perturbation theory appears to be a powerful tool to predict these quantities in the quasilinear regime, it presents severe limitations when extended to smaller scales. In fact, RSD are characterised by a peculiar coupling of large and smallscale clustering that represents a severe challenge to perturbative methods. Several assumptions usually made in the perturbative treatment, such as the irrotational nature of the velocity field, particularly relevant in the RSD modeling, are clearly not valid on small scales.
An insight into the nonlinear evolution of the velocity field is offered by numerical simulations, the standard tool for these kind of investigations. However, while extensive literature is dedicated to the description of the smallscale power spectrum and the characterisation of virialized structures as dark matter halos, the estimation of the velocity field in Nbody simulations presents peculiar challenges. This is due to the challenging estimation of the velocity field in lowdensity regions, resulting in a relatively limited number of studies on this specific topic (see, e.g., Bernardeau & van de Weygaert 1996; Bernardeau et al. 1997a; Pueblas & Scoccimarro 2009; Jennings et al. 2011, 2015; Koda et al. 2014; Zheng et al. 2015; Zhang et al. 2015; Yu et al. 2015; Hahn et al. 2015).
In particular, Jennings (2012), updating the previous results of Jennings et al. (2011), provides a fitting formula for the nonlinear auto and crosspower spectra of the velocity divergence, P_{θθ} and P_{δθ}, calibrated using cosmological Nbody simulations, in terms of the nonlinear matter power spectrum P_{δδ}. Through a fourparameter fit, they reach a 2% accuracy on P_{θθ} for z = 0 and on scales k < 0.65 h Mpc^{−1}; the fit is less accurate for P_{δθ}, and at higher redshifts. A similar but simpler fit for P_{δθ}, again as a function of P_{δδ} is proposed by Zheng et al. (2013). It is expected to be 2% accurate for all scales (and redshift) where the adimensional matter power spectrum Δ_{δδ}(k)≡4πk^{3}P_{δδ}(k) < 1. An alternative fit of the relation of P_{θθ} and P_{θδ} with P_{δδ} is represented by the simple exponential damping proposed instead by Hahn et al. (2015), which reaches a 5% accuracy for k < 1 h Mpc^{−1}.
Any accurate modeling of galaxy clustering aiming at percent precision on the recovered cosmological parameters should also account for a nonzero neutrino mass. There are two main reasons for this. On one side, cosmological observations currently provide the best upper limits on the sum of neutrino masses (Planck Collaboration XIII 2016; PalanqueDelabrouille et al. 2015) and proper modeling is therefore required to extract the most unbiased estimates. On the other hand, at these levels of precision, neglecting the presence of this subdominant component would potentially add a comparable systematic error on any constraint derived from galaxyclustering data (see e.g., Baldi et al. 2014). These are the motivations behind the significant effort spent over the past few years to model the effect of neutrino masses on largescale structure observables, in particular using numerical simulations (see, e.g., AliHaïmoud & Bird 2013; VillaescusaNavarro et al. 2014, 2018; Castorina et al. 2015; Carbone et al. 2016; Inman et al. 2015; Emberson et al. 2017; Zennaro et al. 2017; Liu et al. 2018).
In this paper we use the Dark Energy and Massive Neutrinos Universe (DEMNUni) set of Nbody simulations (Carbone et al. 2016) to model the velocity power spectra (P_{δθ} and P_{θθ}) both in a standard ΛCDM scenario and including a massive neutrino component. The DEMNUni runs represent some of best simulations in massive neutrino cosmologies, both in terms of volume and mass resolution (Castorina et al. 2015). We propose an accurate fitting formula involving a minimal number of free parameters, which we calibrate against simulations, showing that the main dependence on cosmology can be encapsulated in a dependence on the current rms clustering amplitude σ_{8} of the cold matter, that is, CDM and baryons. We focus on this component, as opposed to the total matter including neutrinos, as this appears to provide the simplest description of halo abundance and bias (Castorina et al. 2014).
In Sect. 2 we present our set of Nbody simulations while in Sect. 3 we describe the method used to estimate the cold dark matter velocity field. In Sect. 4 we discuss previous results obtained in past literature and compare them to our results, and in Sect. 5 we summarise our findings in.
2. Simulations
The DEMNUni simulations are a set of Nbody cold dark matter simulations produced with the aim of testing multiple cosmological probes in the presence of massive neutrinos and darkenergy scenarios beyond the standard ΛCDM. They represent a reliable tool for exploring the impact of neutrinos on a wide range of dynamical scales, and have been extended to scenarios including a dynamical darkenergy background, with different equations of state parameters (w_{0}, w_{a}), in order to study their degeneracy with the total neutrino mass at the nonlinear level. The technical implementation and detailed features of the simulations are presented in the description paper by Carbone et al. (2016) and the analysis of cold dark matter clustering by Castorina et al. (2015).
In this work we exploit the first set of simulations, DEMNUniI, describing several flat cosmological models characterised by various values of the total neutrino mass, ∑m_{ν} = 0, 0.17, 0.3, 0.53 eV, while keeping the total matter density parameter fixed at Ω_{m} = 0.32. This implies that the cold dark matter relative density Ω_{cdm} changes across the four simulations in order to keep the sum Ω_{cdm} + Ω_{ν} = Ω_{m} constant. The neutrino density is related to the total neutrino mass as Ω_{ν} = ∑m_{ν}/93.14/h^{2} eV (Lesgourgues & Pastor 2006). Further properties shared by the four cosmologies are the density parameter associated to the cosmological constant Ω_{Λ} = 0.68 and to the baryon density Ω_{b} = 0.05, the Hubble constant H_{0} = 67 km s^{−1}Mpc^{−1}, the primordial spectral index n_{s} = 0.96 and, most importantly, the scalar amplitude of the matter power spectrum A_{s} = 2.1265 × 10^{9}. As a consequence, while in the largescale limit the power spectra of cold dark matter tend to the same value in all cosmological models, the value of the rms of cold dark matter perturbations on spheres of radius 8 h^{−1} Mpc depends on ∑m_{ν}. The latter is denoted as σ_{8, c}, to distinguish it from the rms of total matter perturbations, σ_{8, m}.
The simulations were run on the FERMI supercomputer at CINECA^{1} (5 × 10^{6} CPU hours) using the tree particle hydrodynamical code GADGET3 modified to include massive neutrino particles by Viel et al. (2010). The latter regulates the assembly of N_{cdm} = 2048^{3} cold dark matter particles and N_{ν} = 2048^{3} neutrino particles (when present) within a cubic periodic universe of comoving size L = 2000 h^{−1}Mpc. The mass resolution of cold dark matter and neutrinos varies slightly over the four simulations (values are listed in Table 1), but in all cases it is large enough as to properly describe clustering in the nonlinear regime within the systematic error induced by neglecting baryonic effects. Initial conditions were set at redshift z_{in} = 99 using the Zel’dovich (1970) approximation and were evolved to z = 0, with a softening length ε = 20 h^{−1} kpc. During the runs, 62 snapshots were saved for each simulation, with equal logarithmic interval in the scale factor.
Parameters of the four νΛCDM simulations that change among the different realisations.
3. Measurements of the density and velocity spectra
The estimation of the velocity field in cosmological Nbody simulations has been investigated by Bernardeau & van de Weygaert (1996); Bernardeau et al. (1997b); van de Weygaert et al. (1998); Pueblas & Scoccimarro (2009), or more recently with the Delaunay Tessellation Field Estimator by RomanoDíaz & van de Weygaert (2007) and the Kriging method by Yu et al. (2015). Here we adopt a method close to the one proposed by Pueblas & Scoccimarro (2009), which implements the Delaunay tessellation to reconstruct the velocity field from test particles of cold dark matter. We apply the countincell technique to average both the velocity and density fields within spheres of radius R = 3.9 h^{−1} Mpc. Given the large number of particles in the simulation, we perform the Delaunay tessellation only locally around each cell, rather than applying it to the whole set of cold dark matter particles. We estimate the velocity field on a grid, dividing the simulation into 1024^{3} cubes of size 1.95 h^{−1} Mpc, which are used to index the particle positions.
In order to estimate the velocity within a given spherical cell, we start by considering the particles belonging to the eight subcubes forming a cubical cell that contains the spherical cell itself. We then count how many particles are contained within the sphere; if it is greater than 200 we choose to run the Delaunay tessellation over the cell particles inside the sphere, otherwise we run the Delaunay tesselation inside the cubical cell. In any case we ensure that the total volume covered by the Voronoi tetrahedra inside the spherical cell includes at least 93% of the cell volume, otherwise we automatically extend the radius in which we are keeping the particles to perform the Delaunay tesselation by a factor 3/2. We note that in practice the spherical cell size in unchanged. Only the effective volume of the Delaunay tesselation is varied in order to make sure that the volume fraction of tetrahedra within the spherical cell is representative at the 7% level of the volume of the cell.
The difficulty arising from the estimation of the velocity field convolved with our spherical (tophat) window function is related to the treatment of the tetrahedra which are lying on the boundary of each spherical cell. Our method consists in weighting the average velocity of a tetrahedron by its volume if it lies entirely within the spherical cell. Instead, for tetrahedra that extend outside of the cell boundary, we generate a set of 100 uniformly distributed random points inside the considered tetrahedron and assign them a velocity, linearly interpolated from those at each node of the tetrahedron. The same random points are also used to estimate the volume fraction of the tetrahedron lying inside the spherical cell. In this way, a velocity can be assigned to that specific sector of the spherical cell by averaging the velocities of the included random points. The result is then weighted by the corresponding volume fraction. This process is described in more detail in the following paragraphs and illustrated in Fig. 1.
Fig. 1. Tetrahedra in a spherical cell: the empty circles represent the 3D distribution of CDM particles in a cubic region of size 8 h^{−1} Mpc projected in 2D. They represent the input set of points on which we perform the Delaunay tessellation. The gray circle represents the spherical cell. For clarity, we only show, in dashed light blue lines, the tetrahedra which have their four vertices included in a slice of 1.5 h^{−1} Mpc (in the zdirection) around the center of the cell. As an example we highlight one tetrahedron overlapping the boundary of the spherical cell (dark blue in the bottom lefthand corner). The red and green filled circles represent the randomly distributed points over the tetrahedron, the green ones represent those points lying inside the spherical cell and that are then used to assign a velocity to the corresponding region of the tetrahedron. This is repeated for all tetrahedrons of the tessellation that overlap the boundary of the spherical cell. 

Open with DEXTER 
Each Voronoi tetrahedron can be described by a set of four points x_{1}, x_{2}, x_{3} and x_{4} where x_{i} = (x_{i}, y_{i}, z_{i}); they are therefore defined by the transformation matrix (see Pueblas & Scoccimarro 2009)
where Δx_{i} = x_{i} − x_{1}. As a result, the volume of the jth tetrahedron can be computed as w_{j} = Ψ/6. The matrix Ψ can also be seen as the matrix transforming the basis of the vectors composed of Δx_{2}, Δx_{3} and Δx_{4} into the cartesian basis. If the position s is taken as the tetrahedron basis, then the corresponding coordinates in the cartesian basis are obtained as x = x_{1} + Ψs. The same can be applied in order to interpolate linearly the velocity of a point located in s,
where
For tetrahedra entirely contained inside the sphere one can show that the volume average of the interpolated velocity field inside the tetrahedron is the arithmetic mean of the four velocities taken at each vertex; . As mentioned above, for tetrahedra crossing the cell boundary, we randomly populate their volume with a uniform distribution of N = 100 points (see Fig. 1) and we evaluate their corresponding velocities using Eq. (2). The volumeaveraged velocity assigned to the tetrahedron can be computed as
where N_{in} is the number of random points belonging to the spherical cell and the fraction of its volume inside the spherical cell can be estimated as . Finally, the volumeaveraged velocity assigned to the spherical cell is obtained with the sum
where N_{t} is the total number of tetrahedra with at least one vertex belonging to the spherical cell. We note that we checked explicitly that in the most clustered catalogue, that is, the one corresponding to the ΛCDM cosmology at z = 0, the number of objects in a spherical cell is always greater than 8, thus avoiding empty cells.
Once the velocity and density grids of 512^{3} regularly spaced sampling points have been built, we can Fourier transform them by means of a Fast Fourier Transform algorithm. Since a simple countincell density interpolation can be severely affected by aliasing when transforming to Fourier space, we employ an interlacing technique to reduce this spurious contribution (Hockney & Eastwood 1988; Sefusatti et al. 2016). Regarding the shot noise correction, we neglect it because the mean number of particles in each cell is which corresponds to a 3% contribution to the variance at z = 1.5 and for the M_{ν} = 0.53eV (snapshot having the lowest variance). We then compute the divergence of the velocity field θ_{k} by simply combining the three velocity grids as θ_{k} = i(k_{x}v_{x} + k_{y}v_{y} + k_{z}v_{z}). Then the density power spectrum P_{δδ}, the cross power spectrum P_{δθ}, and the divergence of the velocity power spectrum P_{θθ} are estimated by averaging over spherical shells in kspace. We note that we also average the modes, and assign the value of the angular average of the spectra to the kspace position of the corresponding mode average.
4. Results
Our goal is to provide accurate prescriptions to estimate the P_{θθ} and P_{δθ} auto and crossspectra in the regime where the perturbative approach fails in describing the velocity field and its power spectra (Crocce et al. 2012).
The fitting functions for the velocity power spectra adopted in Jennings et al. (2011) and Jennings (2012) describe the nonlinear velocity spectra P_{θθ} and P_{δθ} in terms of the nonlinear matter power spectrum P_{δδ} assuming a cosmologyindependent relation between these quantities at redshift zero and introducing a scaling relation to extend the results at higher redshift. However, the relations between P_{θθ} and P_{δδ} and between P_{θδ} and P_{δδ} are not universal and depend strongly, in the first place, on the amplitude of linear fluctuations, as measured for example by σ_{8}. This is particularly evident when comparing our set of massive neutrino cosmologies where the neutrino mass directly affects this quantity. In Fig. 2, the velocity spectra (auto and cross) are plotted as a function of the corresponding matter density power spectrum. These relations, far from being universal, clearly depend as much on redshift as they depend on the sum of neutrino masses, via the amplitude suppression induced by the latter. The cosmologyindependence of the fit proposed by Jennings et al. (2011) is perhaps justified by the fact that the different cosmological models considered in that paper share the same amplitude normalisation in terms of the σ_{8} parameter.
Fig. 2. Relation between the velocity cross (left panel) and auto (right panel) spectra with the density autospectrum for the different redshifts and cosmologies considered in this work. The shade of blue represents different values of the total neutrino mass, darker when the neutrino mass is increasing. The various redshifts (z = 0, 0.5, 1, 1.5) are respectively represented with dotdashed, long dashed, short dashed, and solid lines. The black dotted line represents for each redshift the linear mapping ( and taken in the ΛCDM limit, i.e., at small k). 

Open with DEXTER 
We chose a different approach to fit the velocity power spectra measured in the DEMNUniI simulations. In our set of simulations, all cosmological models considered present the same amplitude for the totalmatter power spectrum in the largescale limit, matching current CMB constraints (Planck Collaboration XIII 2016). However, since the presence of massive neutrinos suppresses the growth of fluctuations below the freestreaming scale, we obtain quite different values for σ_{8}(z = 0) as a function of neutrino mass (see Table 1). We are thus able to test a wider range of amplitudes and shapes of the linear power spectrum, since neutrinos affect both of them.
Figure 3 shows the ratio of the measured power spectrum of the density, velocity, and cross to their respective linear predictions. Across all cosmological models and redshifts considered, we notice the usual increase in smallscale power for the density power spectrum P_{δδ} as well as the slower nonlinear growth for the velocity perturbations leading to a nonlinear suppression of P_{δθ} and particularly of P_{θθ} (see e.g., Bernardeau et al. 2002).
Fig. 3. Ratio between nonlinear (measured) P_{NL} and the linear (predicted) P_{Lin} power spectra in the ΛCDM case (∑m_{ν} = 0) and for the three neutrino masses (∑m_{ν} = {0.17, 0.3, 0.53} eV). The density–density, density–velocity divergence, and velocity divergence–velocity divergence spectra are represented by blue solid, red short dashed, and green long dashed lines, respectively. Each row shows the redshift evolution of these ratios (from top to bottom panels: z = 0, 0.5, 1. and 1.5.) 

Open with DEXTER 
Following Hahn et al. (2015), we therefore choose to model the nonlinear corrections to the velocity spectra in terms of damping functions, in order to account for the suppression of power characterising the velocity divergence field. In addition, the level of such suppression will be, as observed, dependent on cosmology. In our approach, we do not assume a universal mapping independently of the considered cosmology. The main motivation supporting this is the empirical evidence that the velocity spectra are damped differently for different cosmological backgrounds.
In this section we explain our general fitting formulae and show how the chosen parameters depend on the value of the overall matter clustering. For simplicity, we first employ fitting functions featuring one single free parameter that accounts for the damping of the linear prediction in the nonlinear regime (see Fig. 3). As first approximation, one can model the velocity spectra using one single damping function such as
and
where the (only) two free parameters are the typical damping scales k_{δ} and k_{θ}. We note that refers to the nonlinear densitydensity CDM power spectrum computed from the Halofit calibration of Takahashi et al. (2012) while refers to the linear autospectrum of the velocity divergence which can be computed as .
The fit for the velocity power spectra is carried out using a leastsquares approach, that is, we compute the likelihood of the parameters given the measured spectra with the χ^{2} function defined as
where N is the total number of wavenumbers considered in the fit, is the measured auto or crossspectrum and is its variance at the i−th wave mode k_{i}. We limit the fitting range to k_{max} = 0.6 h Mpc^{−1}, and since we have only one realisation for each cosmology we neglect the shotnoise and the nonGaussian contributions to the covariance between wave modes. Under these assumptions the error can be approximated as
We note that our choice of using a fixed power α = 1 for the exponent, rather than having more degrees of freedom (i.e., exp(−(k/k^{*})^{α})), comes from a further test we carried out, which shows that the bestfit value of α is always close to 1 if we treat it as a free parameter (see also Hahn et al. 2015). We note that this is at odds with what one would expect when considering the propagator in renormalized perturbation theory (Bernardeau et al. 2008, 2012; Crocce et al. 2012), in which case α would be closer to 2.
The results of the fit are shown in the first and third columns of Fig. 4. One can see that the simple modeling provided by Eq. (6) is able to reproduce the cross power spectrum P_{δθ} in the nonlinear regime with an accuracy better than 5% up to k = 0.65 h Mpc^{−1}. The fit is not working as well for the auto spectrum P_{θθ} (third panel) especially at low redshift. The inaccuracy in this case can reach 7–8% in between k = 0.1 h Mpc^{−1} and k = 0.3 h Mpc^{−1}. Nonetheless, this approximation could be considered sufficient for analyses that do not require precision around the BAO scale of better than few percent. For more general applications, we improve the accuracy of the model by increasing the degrees of freedom of the fitting functions.
Fig. 4. Fractional deviation of the measured velocity spectra from our bestfitting models. Different rows mark different redshifts as specified on the left of the panels. The first and third columns show the fit of P_{δθ} using Eq. (6), and the fit of P_{θθ} using Eq. (7). The second and fourth columns show results corresponding to Eqs. (10) and (11). The black dashed line marks the 0, while the light/dark gray bands represent the 5%/3% deviations from the measurements. Finally, the red shaded color shows the results for various neutrino masses (darker when increasing the neutrino mass) using the fitted shape parameters while the blue shaded lines show the results when using the fitted dependance of the shape parameters with respect to σ_{8, m} (see Eq. (12)). 

Open with DEXTER 
In the case of the crosspower spectrum, it is sufficient to add only one parameter b, which we fit between k = 0.55 and k = 0.7 h Mpc^{−1},
For the autospectrum P_{θθ}, it turns out that two extra free parameters are required for a proper gain in accuracy. We therefore adopt a polynomial fitting for the damping function, as
involving three parameters, which we fit in the range 0.01 < k < 0.7 h Mpc^{−1}. The performances of the new fitting functions are shown in the second and fourth columns of Fig. 4. In this case, the measurements of P_{θθ} are reproduced with a maximum systematic error of 3% on all scales below k = 0.7 h Mpc^{−1}, and for all redshifts and neutrino masses considered. We therefore use a total of only five free parameters k_{δ}, b, a_{1}, a_{2} and a_{3} to mimic, with a 3% level accuracy at k < 0.7 h Mpc^{−1}, the nonlinear effects on both the auto and crossspectra P_{θθ} and P_{δθ}.
Let us now focus on the sensitivity of these parameters to cosmology and specifically to the overall amplitude of the matter power spectrum. From our set of four simulations at four different redshifts we are able to span a large range of possible values of σ_{8}, which we use as a proxy for the amount of nonlinearities. Regarding neutrino cosmologies it is necessary to choose whether we use the σ_{8} defined for cold dark matter only, σ_{8, c}, or the one defined for the total matter, σ_{8, m}. It has been shown (Castorina et al. 2015) that regarding the bias or the nonlinear effects on the density power spectrum P_{δδ}, what matters is the amplitude of the cold dark matter clustering and not the total one. We have analyzed the dependency of the fitted parameters with respect to both amplitudes σ_{8, c} and σ_{8, m} and we found that one should use the total matter σ_{8, m} parameter in order to assess the correct values of k_{δ}, b, a_{1}, a_{2} and a_{3}; at least, this choice is the one that lowers the residual dependency with respect to the neutrino mass. Figure 5 shows that the cosmological dependance of the fitting parameters can be mostly encapsulated into the σ_{8, m} parameter evaluated at the corresponding redshift. The most relevant example is the good match at redshift 1.5 of the ΛCDM simulation with the M_{ν} = 0.53eV simulation at redshift 1; these simulations share almost the same value of the total matter σ_{8, m}, while differing significantly in the power spectrum shape. For the two cases, we obtain the same values for the fitted parameters k_{δ} and k_{θ}.
Fig. 5. Top panel: dependence of the fitting parameters k_{δ}, b, k_{θ}, a_{1}, a_{2} and a_{3} on the total matter clustering amplitude σ_{8,m} at each redshift, in the corresponding simulation (symbols). Lines are showing the fit of the dependance. Bottom panel: residual between the fitted σ_{8,m} dependance and the measured one. 

Open with DEXTER 
On the contrary, when using the cold dark matter σ_{8,c} parameter (see Table 1) the residuals are increasing because of spurious neutrino mass dependence when going from one cosmology to another. This effect can be explained by the fact that if massive neutrinos have a weak effect on the dark matter clustering in the nonlinear regime, they add, instead, a relevant contribution to the velocity field which is felt by dark matter particles. As a result, it seems that the preferred dependence on the total matter σ_{8,m} comes from the fact that nonlinearities in the velocity field are generated by the whole matter distribution (cold dark matter plus neutrinos). This confirms the need for running cosmological simulations including massive neutrino particles in order to generate a velocity field correctly treated in the nonlinear regime, especially for what concerns RSD analyses.
The final step of our fitting process is to fit the dependence of the shape parameters k_{δ}, b, a_{1}, a_{2}, a_{3} and k_{θ} with respect to the total matter σ_{8,m}. To this purpose we limit ourselves to linear or quadratic fitting depending on the parameters, finding
where σ_{8,m} refers to the linear rms of total matter fluctuations computed at the required redshift. Therefore, the cross and autospectra P_{δθ} and P_{θθ} can be computed as follows: compute the linear and nonlinear cold dark matter power spectrum at the required redshift, evaluate the linear σ_{8,m} as
where R = 8h^{−1}Mpc, and P_{m} is the linear total matter power spectrum. Finally, compute the k_{δ}, b, k_{θ}, a_{1}, a_{2} and a_{3} parameters from Eq. (12) and the velocity spectra from Eqs. (10) and (11) (or (6) and (7) depending on the required accuracy). In order to summarize the overall accuracy of our fitting formula, we show in Fig. 4 the comparison between the intrinsic accuracy of the fitting formula (for the shape) in red shade and the final accuracy obtained assuming the additional fitted dependency on σ_{8,m} in blue shade. One can see that the accuracy below k = 0.8 h Mpc^{−1} is about 3% for the crosspower spectrum while the autopower spectrum reaches a similar accuracy below k = 0.7 h Mpc^{−1}.
5. Summary
We set up an original algorithm in order to estimate the velocity field in cosmological Nbody simulations. From those measurements performed on sixteen particle distributions spanning four different cosmological models and four redshifts, that is, a series of 16 snapshots, we have shown that the mapping from the nonlinear CDM density power spectrum P_{δδ} to the nonlinear velocity spectra P_{δθ} and P_{θθ} cannot be considered as universal but shows a clear dependence on the amplitude of dark matter clustering (see Fig. 2).
Adopting a very simple modeling involving only two free parameters, we managed to reproduce our measurements with a precision of 3% below k = 0.6 h Mpc^{−1} for the crosspower spectrum P_{δθ}. However, we reached only a 5% accuracy at redshifts higher than or equal to z = 0.5 and for k lower than 0.7 h Mpc^{−1}. We then present an improved version of the fitting function which involves a total of five free parameters (two for the crossspectrum and three for the autospectrum) and allows us to obtain an overall accuracy of 3% for wave modes lower than 0.7 h Mpc^{−1}.
Finally, we we showed evidence of the dependence of the shape parameters of the proposed fitting functions on the total matter rms fluctuations in spheres of radius 8 h^{−1} Mpc and proposed simple fitting forms for the shape parameters with respect to the σ_{8,m} parameter evaluated at the considered redshift. This preferred dependence with respect to the total matter amplitude rather than the cold dark matter one confirms the relevance of the neutrino perturbations on the cold dark matter velocity field in the nonlinear regime.
In a future paper we shall generalize these results using the second set of the DEMNUni cosmological simulations that include dynamical dark energy, extending our study of the dependence of the shape parameters on the rms amplitude of clustering σ_{8,m}.
Acknowledgments
This work were developed in the framework of the ERC Darklight project, supported by an Advanced Research Grant of the European Research Council (n. 291521). CC, LG and ES also acknowledge a contribution from ASI through grant I/023/12/0. ES is also supported by INFN INDARK PD51 grant. The DEMNUniI simulations were carried out in the framework of the “The Dark Energy and MassiveNeutrino Universe” using the Tier0 IBM BG/Q Fermi machine of the Centro Interuniversitario del NordEst per il Calcolo Elettronico (CINECA). We acknowledge a generous CPU allocation by the Italian SuperComputing Resource Allocation (ISCRA).
References
 AliHaïmoud, Y., & Bird, S. 2013, MNRAS, 428, 3375 [NASA ADS] [CrossRef] [Google Scholar]
 Amendola, L., Appleby, S., Avgoustidis, A., et al. 2018, Liv. Rev. Rel., 21, 2 [Google Scholar]
 Baldi, M., VillaescusaNavarro, F., Viel, M., et al. 2014, MNRAS, 440, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Bassett, B., & Hlozek, R. 2010, in Baryon acoustic oscillations, ed. P. RuizLapuente, 246 [Google Scholar]
 Bernardeau, F., & van de Weygaert, R. 1996, MNRAS, 279, 693 [NASA ADS] [CrossRef] [Google Scholar]
 Bernardeau, F., van de Weygaert, R., Hivon, E., & Bouchet, F. R. 1997a, MNRAS, 290, 566 [NASA ADS] [CrossRef] [Google Scholar]
 Bernardeau, F., van de Weygaert, R., Hivon, E., & Bouchet, F. R. 1997b, MNRAS, 290, 566 [NASA ADS] [CrossRef] [Google Scholar]
 Bernardeau, F., Colombi, S., Gaztañaga, E., & Scoccimarro, R. 2002, Phys. Rev., 367, 1 [Google Scholar]
 Bernardeau, F., Crocce, M., & Scoccimarro, R. 2008, Phys. Rev. D, 78, 103521 [NASA ADS] [CrossRef] [Google Scholar]
 Bernardeau, F., Crocce, M., & Scoccimarro, R. 2012, Phys. Rev. D, 85, 123519 [NASA ADS] [CrossRef] [Google Scholar]
 Bianchi, D., Percival, W. J., & Bel, J. 2016, MNRAS, 463, 3783 [NASA ADS] [CrossRef] [Google Scholar]
 Bose, B., & Koyama, K. 2016, J. Cosmol. Astropart. Phys., 8, 032 [NASA ADS] [CrossRef] [Google Scholar]
 Carbone, C., Petkova, M., & Dolag, K. 2016, J. Cosmol. Astropart. Phys., 7, 034 [NASA ADS] [CrossRef] [Google Scholar]
 Castorina, E., Sefusatti, E., Sheth, R. K., VillaescusaNavarro, F., & Viel, M. 2014, J. Cosmol. Astropart. Phys., 2, 49 [NASA ADS] [CrossRef] [Google Scholar]
 Castorina, E., Carbone, C., Bel, J., Sefusatti, E., & Dolag, K. 2015, J. Cosmol. Astropart. Phys., 7, 043 [NASA ADS] [CrossRef] [Google Scholar]
 Crocce, M., Scoccimarro, R., & Bernardeau, F. 2012, MNRAS, 427, 2537 [NASA ADS] [CrossRef] [Google Scholar]
 Emberson, J. D., Yu, H.R., Inman, D., et al. 2017, Res. Astron. Astrophys., 17, 085 [NASA ADS] [CrossRef] [Google Scholar]
 Fonseca de la Bella, L., Regan, D., Seery, D., & Hotchkiss, S. 2017, J. Cosmol. Astropart. Phys., 11, 039 [NASA ADS] [CrossRef] [Google Scholar]
 Guzzo, L., Pierleoni, M., Meneux, B., et al. 2008, Nature, 451, 541 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Hahn, O., Angulo, R. E., & Abel, T. 2015, MNRAS, 454, 3920 [NASA ADS] [CrossRef] [Google Scholar]
 Hand, N., Seljak, U., Beutler, F., & Vlah, Z. 2017, J. Cosmol. Astropart. Phys., 10, 009 [NASA ADS] [CrossRef] [Google Scholar]
 Hashimoto, I., Rasera, Y., & Taruya, A. 2017, Phys. Rev. D, 96, 043526 [NASA ADS] [CrossRef] [Google Scholar]
 Hockney, R. W., & Eastwood, J. W. 1988, Computer Simulation using Particles (New York: Taylor & Francis) [CrossRef] [Google Scholar]
 Inman, D., Emberson, J. D., Pen, U.L., et al. 2015, Phys. Rev. D, 92, 023502 [NASA ADS] [CrossRef] [Google Scholar]
 Jennings, E. 2012, MNRAS, 427, L25 [NASA ADS] [Google Scholar]
 Jennings, E., Baugh, C. M., & Pascoli, S. 2011, MNRAS, 410, 2081 [NASA ADS] [Google Scholar]
 Jennings, E., Baugh, C. M., & Hatt, D. 2015, MNRAS, 446, 793 [NASA ADS] [CrossRef] [Google Scholar]
 Kaiser, N. 1987, MNRAS, 227, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Koda, J., Blake, C., Davis, T., et al. 2014, MNRAS, 445, 4267 [NASA ADS] [CrossRef] [Google Scholar]
 Lesgourgues, J., & Pastor, S. 2006, Phys. Rev., 429, 307 [Google Scholar]
 Liu, J., Bird, S., Zorrilla Matilla, J. M., et al. 2018, J. Cosmol. Astropart. Phys., 3, 049 [NASA ADS] [CrossRef] [Google Scholar]
 Matsubara, T. 2008a, Phys. Rev. D, 78, 083519 [NASA ADS] [CrossRef] [Google Scholar]
 Matsubara, T. 2008b, Phys. Rev. D, 77, 063530 [NASA ADS] [CrossRef] [Google Scholar]
 Okumura, T., Hand, N., Seljak, U., Vlah, Z., & Desjacques, V. 2015, Phys. Rev. D, 92, 103516 [NASA ADS] [CrossRef] [Google Scholar]
 PalanqueDelabrouille, N., Yèche, C., Baur, J., et al. 2015, J. Cosmol. Astropart. Phys., 11, 011 [NASA ADS] [CrossRef] [Google Scholar]
 Perko, A., Senatore, L., Jennings, E., & Wechsler, R. H. 2016, ArXiv eprints [arXiv:1610.09321] [Google Scholar]
 Pezzotta, A., de la Torre, S., Bel, J., et al. 2017, A&A, 604, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pueblas, S., & Scoccimarro, R. 2009, Phys. Rev. D, 80, 043504 [NASA ADS] [CrossRef] [Google Scholar]
 Reid, B. A., & White, M. 2011, MNRAS, 417, 1913 [NASA ADS] [CrossRef] [Google Scholar]
 RomanoDíaz, E., & van de Weygaert, R. 2007, MNRAS, 382, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Sánchez, A. G., Scoccimarro, R., Crocce, M., et al. 2017, MNRAS, 464, 1640 [NASA ADS] [CrossRef] [Google Scholar]
 Sato, M., & Matsubara, T. 2011, Phys. Rev. D, 84, 043501 [NASA ADS] [CrossRef] [Google Scholar]
 Scoccimarro, R. 2004, Phys. Rev. D, 70, 083007 [NASA ADS] [CrossRef] [Google Scholar]
 Sefusatti, E., Crocce, M., Scoccimarro, R., & Couchman, H. M. P. 2016, MNRAS, 460, 3624 [NASA ADS] [CrossRef] [Google Scholar]
 Seljak, U., & McDonald, P. 2011, J. Cosmol. Astropart. Phys., 11, 039 [NASA ADS] [CrossRef] [Google Scholar]
 Senatore, L., & Zaldarriaga, M. 2014, ArXiv eprints [arXiv:1409.1225] [Google Scholar]
 Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152 [NASA ADS] [CrossRef] [Google Scholar]
 Taruya, A., Nishimichi, T., Saito, S., & Hiramatsu, T. 2009, Phys. Rev. D, 80, 123503 [NASA ADS] [CrossRef] [Google Scholar]
 Taruya, A., Nishimichi, T., & Saito, S. 2010, Phys. Rev. D, 82, 063522 [NASA ADS] [CrossRef] [Google Scholar]
 Taruya, A., Nishimichi, T., & Bernardeau, F. 2013, Phys. Rev. D, 87, 083509 [NASA ADS] [CrossRef] [Google Scholar]
 Uhlemann, C., Kopp, M., & Haugg, T. 2015, Phys. Rev. D, 92, 063004 [NASA ADS] [CrossRef] [Google Scholar]
 Upadhye, A., Kwan, J., Pope, A., et al. 2016, Phys. Rev. D, 93, 063515 [NASA ADS] [CrossRef] [Google Scholar]
 Valageas, P. 2011, A&A, 526, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van de Weygaert, R., & Bernardeau, F. 1998, in Large Scale Structure: Tracks and Traces, eds. V. Mueller, S. Gottloeber, J. P. Muecket, & J. Wambsganss, 207 [Google Scholar]
 Viel, M., Haehnelt, M. G., & Springel, V. 2010, J. Cosmol. Astropart. Phys., 6, 015 [NASA ADS] [CrossRef] [Google Scholar]
 VillaescusaNavarro, F., Marulli, F., Viel, M., et al. 2014, J. Cosmol. Astropart. Phys., 3, 11 [NASA ADS] [CrossRef] [Google Scholar]
 VillaescusaNavarro, F., Banerjee, A., Dalal, N., et al. 2018, ApJ, 861, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Vlah, Z., Castorina, E., & White, M. 2016, J. Cosmol. Astropart. Phys., 12, 007 [NASA ADS] [CrossRef] [Google Scholar]
 Weinberg, D. H., Mortonson, M. J., Eisenstein, D. J., et al. 2013, Phys. Rev., 530, 87 [Google Scholar]
 Yu, Y., Zhang, J., Jing, Y., & Zhang, P. 2015, Phys. Rev. D, 92, 083527 [NASA ADS] [CrossRef] [Google Scholar]
 Zel’dovich, Y. B. 1970, A&A, 5, 84 [NASA ADS] [Google Scholar]
 Zennaro, M., Bel, J., VillaescusaNavarro, F., et al. 2017, MNRAS, 466, 3244 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, P., Pan, J., & Zheng, Y. 2013, Phys. Rev. D, 87, 063526 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, P., Zheng, Y., & Jing, Y. 2015, Phys. Rev. D, 91, 043522 [NASA ADS] [CrossRef] [Google Scholar]
 Zheng, Y., Zhang, P., Jing, Y., Lin, W., & Pan, J. 2013, Phys. Rev. D, 88, 103510 [NASA ADS] [CrossRef] [Google Scholar]
 Zheng, Y., Zhang, P., & Jing, Y. 2015, Phys. Rev. D, 91, 043523 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Parameters of the four νΛCDM simulations that change among the different realisations.
All Figures
Fig. 1. Tetrahedra in a spherical cell: the empty circles represent the 3D distribution of CDM particles in a cubic region of size 8 h^{−1} Mpc projected in 2D. They represent the input set of points on which we perform the Delaunay tessellation. The gray circle represents the spherical cell. For clarity, we only show, in dashed light blue lines, the tetrahedra which have their four vertices included in a slice of 1.5 h^{−1} Mpc (in the zdirection) around the center of the cell. As an example we highlight one tetrahedron overlapping the boundary of the spherical cell (dark blue in the bottom lefthand corner). The red and green filled circles represent the randomly distributed points over the tetrahedron, the green ones represent those points lying inside the spherical cell and that are then used to assign a velocity to the corresponding region of the tetrahedron. This is repeated for all tetrahedrons of the tessellation that overlap the boundary of the spherical cell. 

Open with DEXTER  
In the text 
Fig. 2. Relation between the velocity cross (left panel) and auto (right panel) spectra with the density autospectrum for the different redshifts and cosmologies considered in this work. The shade of blue represents different values of the total neutrino mass, darker when the neutrino mass is increasing. The various redshifts (z = 0, 0.5, 1, 1.5) are respectively represented with dotdashed, long dashed, short dashed, and solid lines. The black dotted line represents for each redshift the linear mapping ( and taken in the ΛCDM limit, i.e., at small k). 

Open with DEXTER  
In the text 
Fig. 3. Ratio between nonlinear (measured) P_{NL} and the linear (predicted) P_{Lin} power spectra in the ΛCDM case (∑m_{ν} = 0) and for the three neutrino masses (∑m_{ν} = {0.17, 0.3, 0.53} eV). The density–density, density–velocity divergence, and velocity divergence–velocity divergence spectra are represented by blue solid, red short dashed, and green long dashed lines, respectively. Each row shows the redshift evolution of these ratios (from top to bottom panels: z = 0, 0.5, 1. and 1.5.) 

Open with DEXTER  
In the text 
Fig. 4. Fractional deviation of the measured velocity spectra from our bestfitting models. Different rows mark different redshifts as specified on the left of the panels. The first and third columns show the fit of P_{δθ} using Eq. (6), and the fit of P_{θθ} using Eq. (7). The second and fourth columns show results corresponding to Eqs. (10) and (11). The black dashed line marks the 0, while the light/dark gray bands represent the 5%/3% deviations from the measurements. Finally, the red shaded color shows the results for various neutrino masses (darker when increasing the neutrino mass) using the fitted shape parameters while the blue shaded lines show the results when using the fitted dependance of the shape parameters with respect to σ_{8, m} (see Eq. (12)). 

Open with DEXTER  
In the text 
Fig. 5. Top panel: dependence of the fitting parameters k_{δ}, b, k_{θ}, a_{1}, a_{2} and a_{3} on the total matter clustering amplitude σ_{8,m} at each redshift, in the corresponding simulation (symbols). Lines are showing the fit of the dependance. Bottom panel: residual between the fitted σ_{8,m} dependance and the measured one. 

Open with DEXTER  
In the text 