Massive Neutrinos and Galaxy Clustering in f(R) Gravity Cosmologies

Abstract

Cosmic neutrinos have been playing a key role in cosmology since the discovery of their mass. They can affect cosmological observables and have several implications being the only hot dark matter candidates that we currently know to exist. The combination of massive neutrinos and an adequate theory of gravity provide a perfect scenario to address questions on the dark sector that have remained unanswered for years. In particular, in the era of precision cosmology, galaxy clustering and redshift-space distortions afford one of the most powerful tools to characterise the spatial distribution of cosmic tracers and to extract robust constraints on neutrino masses. In this chapter, we study how massive neutrinos affect the galaxy clustering and investigate whether the cosmological effects of massive neutrinos might be degenerate with f(R) gravity cosmologies, which would severely affect the constraints.

1. Introduction

From first principles, it is well known that a theory of gravity is needed to describe the spatial properties and dynamics of the large-scale structures (LSS) of the universe. The observational data collected for several decades provide strong support to the concordance model Lambda cold dark matter (ΛCDM), which yields a consistent description of the main properties of the LSS [1, 2, 3, 4]. However, since cosmological observations have entered in an unprecedented precision era, one of the current aims is to test some of the most fundamental assumptions of the concordance model of the universe. In this sense, the ΛCDM model assumes (i) the general theory of relativity (GR) as the theory describing gravitational interactions at large scales, (ii) the standard model of particles and (iii) the cosmological principle. Moreover, in this framework, the universe is currently dominated by dark energy (DE) in the form of a cosmological constant, responsible for the late-time cosmic acceleration [5, 6, 7] and by a cold dark matter (CDM) component that drives the formation and evolution of cosmic structures.

Recently, several shortcomings have been found in the ΛCDM scenario, like a possible tension in the parameter constraint of H ₀ and ${\sigma }_8$ when different probes are used [4, 8, 9, 10]. This has motivated the interest on theoretical models beyond GR. Among them, the models based on f(R) gravity are the favourite ones because of their generality and rich phenomenology [11, 12]. Moreover, modified gravity (MG) models represent one of the most viable alternatives to explain cosmic acceleration [13] that require satisfying simultaneously solar system constraints and to be consistent with the measured accelerated cosmic expansion and large-scale constraints [14, 15, 16, 17]. An extra motivation to study MG models is given by the fact that massive neutrinos, the only (hot) dark matter candidates we actually know to exist, can affect these observables and have several cosmological implications [18]. However, the degeneracy between some MG models and the total neutrino mass [19, 20, 21] give rise to a limitation of many standard cosmological statistics [22, 23]. In this context, the clustering analysis and redshift-space galaxy clustering have been proven to be a powerful cosmological probe to discriminate among MG scenarios with massive neutrinos, as will be discussed in the following sections.

Regarding the dynamics of background universe in the standard framework, it is well described by the Friedmann-Lemaître-Robertson-Walker (FLRW) metric, whose line element in natural units $c=1$ is given by

$${\mathit{ds}}^2=-{\mathit{dt}}^2+a^2\left(t\right)\left[\frac{{\mathit{dr}}^2}{1-{\mathit{kr}}^2}+r^{2}d{\theta }^2+r^2{sin}^2\mathit{\theta d}{\varphi }^2\right].$$ E1

where a(t) is the so-called scale factor, ( $r$ , $\theta$ , $\varphi$ ) are dimensionless spherical-polar coordinates and k defines the geometry of the universe under consideration to be flat ( $k=0$ ), open ( $k<0$ ) or closed ( $k>0$ ). The equations of motion that describe the time evolution of a(t) and the dynamic growth of the universe are called Friedmann equations and are given by

$${\left(\frac{\dot{a}}{a}\right)}^2=\frac{8\mathit{\pi G}}{3}\rho -\frac{k}{a^2}+\frac{\mathrm{\Lambda }}{3},$$ E2

$$2\frac{\ddot{a}}{a}+{\left(\frac{\dot{a}}{a}\right)}^2=-8\mathit{\pi Gp}-\frac{k}{a^2}+\mathrm{\Lambda },$$ E3

that can be re-expressed as

$$\dot{\rho }=-3\frac{\dot{a}}{a}\left(\rho +p\right),$$ E4

$$\frac{\ddot{a}}{a}=-\frac{4\mathit{\pi G}}{3}\left(\rho +3p\right)+\frac{\mathrm{\Lambda }}{3}$$ E5

after eliminating the curvature term. These equations lead to the definition of the Hubble parameter as $H\equiv \frac{\dot{a}}{a}$ , which drives the expansion rate of the universe and usually is represented in terms of the dimensionless factor h, defined by the expression H₀ = 100 h km s⁻¹ Mp c⁻¹ at the present epoch. In order to have a full description of the background universe, it is necessary an equation of state (EoS) of the cosmic fluid, considering that it has three principal components: baryonic matter, dark matter and radiation (see Figure 1a). A first approximation consists in assuming a linear relationship between $\rho$ and $p$ ; thus, the equation of state can be written as $p=\mathit{w\rho }$ , where w is a parameter that in principle can be time-dependent $w=w\left(t\right)$ , but the simplest approach is to consider it as a constant. Under this assumption, Eq. (4) is easy to solve, resulting in $\rho \left(a\right)=a^{-3\left(1+w\right)}$ . The case with $w=-1$ corresponds to the so-called cosmological constant; it is obtained assuming a constant density energy so that the corresponding EoS is $p=-\rho$ . The case of a flat universe ( $k=0$ ) is interesting because of its agreement with many observational results [4]; it also implies a special value of the matter density in the universe that allows to introduce naturally a critical density ${\rho }_{\mathit{cr}}$ in terms of the Hubble parameter. The critical density is also useful to define the dimensionless density parameter ${\mathrm{\Omega }}_i$ for the various species i, related to radiation ${\mathrm{\Omega }}_r$ , matter ${\mathrm{\Omega }}_m$ , cosmological constant ${\mathrm{\Omega }}_{\mathrm{\Lambda }}$ and curvature ${\mathrm{\Omega }}_k$ . It is easy to verify that the sum of these parameters is equal to unity, as it can be expected from the Friedmann equations; in fact, ${\mathrm{\Omega }}_m+{\mathrm{\Omega }}_r+{\mathrm{\Omega }}_{\mathrm{\Lambda }}+{\mathrm{\Omega }}_k=1$ is known as the cosmic sum rule. Therefore, a Friedmann universe can be described by the cosmological parameters $\left(H_0,{\mathrm{\Omega }}_m,{\mathrm{\Omega }}_r,{\mathrm{\Omega }}_k,{\mathrm{\Omega }}_{\mathrm{\Lambda }}\right)$ such that the expansion rate as a function of the scale factor is given by

Figure 1.

(a) The current composition of the universe derived from Planck data [4], (b) evolution of the energy density parameter $\rho$ normalized to ${\rho }_{\mathit{cr},0}$ as a function of the scale factor a.

$$H^2\left(a\right)=H_0^2\left[{\mathrm{\Omega }}_r{a}^{-4}+{\mathrm{\Omega }}_m{a}^{-3}+{\mathrm{\Omega }}_k{a}^{-2}+{\mathrm{\Omega }}_{\mathrm{\Lambda }}\right].$$ E6

This equation is usually written as a dimensionless function defined by $E\left(a\right)=H\left(a\right)/H_0$ . The last constraints on cosmological parameters obtained by the Planck satellite show that $w=-1.03\pm 0.03$ , consistent with a cosmological constant, and ${\mathrm{\Omega }}_m=0.3153\pm 0.0073,{\mathrm{\Omega }}_{\mathrm{\Lambda }}=0.692\pm 0.012,{\mathrm{\Omega }}_k=0.000\pm 0.005$ , where ${\mathrm{\Omega }}_m$ contains the density of baryons $\left({\mathrm{\Omega }}_b\right)$ and cold dark matter $\left({\mathrm{\Omega }}_{\mathit{CDM}}\right)$ . Additionally, in the last few years, it became usual to include the energy density of neutrinos ${\mathrm{\Omega }}_{\nu }$ ; they contribute to the radiation density at early times but behave as matter after the non-relativistic transition at late times [24], so that for a flat universe, the total energy density is given by $\rho ={\rho }_{\gamma }+{\rho }_{\mathrm{CDM}}+{\rho }_{\mathrm{b}}+{\rho }_{\nu }+{\rho }_{\mathrm{\Lambda }}$ . Figure 1a shows the percentages, derived from [4] data, in which each species contributes to the total content of the universe. Figure 1b shows the evolution of the density parameter $\rho$ (in units of ${\rho }_{\mathit{cr},0}$ ) as a function of the scale factor a. For further details concerning the background universe, see Refs. [25, 26, 27].

2. Modified gravity and massive neutrinos

2.1 Modified gravity models

One of the most interesting modifications of GR is that which modifies the Einstein-Hilbert action by introducing a scalar function, f(R), as follows:

$$S=\int {\mathrm{d}}^{4}x\sqrt{-g}\left(\frac{R+f\left(R\right)}{16\mathit{\pi G}}+{\mathcal{L}}_m\right),$$ E7

where R is the Ricci scalar, G is the Newton’s gravitational constant, g is the determinant of the metric tensor $g_{\mathit{\mu \nu }}$ and ${\mathcal{L}}_m$ is the Lagrangian density of all matter fields. For a classification of f(R) theories of gravity and the assumptions needed to arrive to the various versions of f(R) gravity and GR, see, e.g., [28]. Thus, for a general f(R) model, one can consider a spatially flat FLRW universe with metric ${\mathit{ds}}^2=-{\mathit{dt}}^2+a^2\left(t\right)d{\mathbf{x}}^2$ , so that varying the Einstein-Hilbert action with respect to $g_{\mathit{\mu \nu }}$ , one can get a general form of the modified Einstein field equations. Consequently, the corresponding modified Friedmann equations are given by

$$3H^2\left(1+f_R\right)=8\mathit{\pi G}\left({\rho }_m+{\rho }_{\mathit{rad}}\right)+\frac{1}{2}\left(f_{R}R-f\right)-3H{\dot{f}}_R$$ E8

$$-2\left(1+f_R\right)\dot{H}=8\mathit{\pi G}\left({\rho }_m+\frac{4}{3}{\rho }_{\mathit{rad}}\right)+{\ddot{f}}_R-H{\dot{f}}_R,$$ E9

where $R=6\left(2H^2+\dot{H}\right)$ and the over-dot denotes a derivative with respect to the cosmic time t. In general, the background evolution of a viable f(R) is not simple as it has been shown by [29, 30, 31]. However it is possible to get an approximation in a way that is analogous to the DE models, by neglecting the higher derivative and the non-linear terms. By defining the growth rate as $f\equiv dln\delta /dlna$ , the equation that describes the growth of matter perturbations in terms of the density contrast $\delta$ in a f(R) model is approximated by [29, 32]

$$-3{\mathrm{\Omega }}_m\left(a\right)\left(1-{\mathrm{\Omega }}_m\left(a\right)\right)\frac{\mathit{df}}{d{\mathrm{\Omega }}_m\left(a\right)}+f^2+\left(2-\frac{3}{2}{\mathrm{\Omega }}_m\left(a\right)\right)f=\frac{3}{2}\frac{G_{\mathrm{eff}}}{G}{\mathrm{\Omega }}_m\left(a\right).$$ E10

where G_eff is the effective gravitational constant that can be written as [33]. A plausible f(R) function able to satisfy the solar system constraints, to mimic the ΛCDM model at high-redshift regime where it is well tested by the CMB and, at the same time, to accelerate the expansion of the universe at low redshift but without a cosmological constant [34], suggests that

$$\underset{R\to \infty }{lim}f\left(R\right)=\text{const},\underset{R\to 0}{lim}f\left(R\right)=0,$$ E11

which can be satisfied by a broken power law function such that

$$f\left(R\right)=-m^2\frac{c_1{\left(\frac{R}{m^2}\right)}^n}{c_2{\left(\frac{R}{m^2}\right)}^n+1},$$ E12

where the mass scale m is defined as $m^2\equiv H_0^2{\mathrm{\Omega }}_{\mathrm{M}}$ and c ₁, c ₂ and n are non-negative free parameters of the model [34]. For this f(R) model, the background expansion history is consistent with the ΛCDM case by choosing $c_1/c_2=6{\mathrm{\Omega }}_{\mathrm{\Lambda }}/{\mathrm{\Omega }}_{\mathrm{M}}$ , where ${\mathrm{\Omega }}_{\mathrm{\Lambda }}$ and ${\mathrm{\Omega }}_{\mathrm{M}}$ are the dimensionless density parameters for vacuum and matter, respectively.

Nowadays it is generally accepted that some MG theories such as the [34] f(R) are strongly degenerated in a wide range of their observables with the effects of massive neutrinos; see, e.g., [19, 20, 21, 35]. This represents a serious challenge constraining cosmological models from current and future galaxy surveys requiring robust and reliable methods to disentangle both phenomena. Furthermore, for some specific combinations of the f(R) function and of the total neutrino mass $m_{\nu }\equiv \mathrm{\Sigma }{m}_{\nu ,i}$ , standard statistics would not distinguish them from the standard ΛCDM expectations (see [21, 22, 36]). In addition, since the degeneracy is mostly driven by the non-linear behaviour of both the MG and the massive neutrinos effects on the LSS, the linear tools are not suitable to properly disentangle the combined parameter space [36].

2.2 Massive neutrinos and the large-scale structure

Motivated by the apparent violation of energy, momentum and spin in $\beta$ -decay processes, Pauli proposed the existence of neutrinos in 1930 to keep the conservation laws safe. Eventually, 26 years after they have been theoretically postulated, the neutrinos were detected for the first time by Cowan [37]. Neutrinos are classified in three ‘flavours’ in the standard model of particles; they were considered to be massless for some time until the discovery of the neutrino oscillation phenomena, i.e. related to the change of flavour [38]. Since then, it is known that at least two of the three neutrino families are massive, in contrast to the particle standard model assumption; however, measuring the absolute masses of the neutrinos is not easy, which makes this a very active field of research today, both for cosmology and particle physics. In a cosmological context, the neutrinos leave detectable imprints on observations that can then be used to constrain their properties; in particular, the presence of massive neutrinos impacts the background evolution of the universe and the growth of structures [39]. In the early universe, massive neutrinos are relativistic and indistinguishable from the massless ones, behaving like photons, meaning that their energy density drops like $\propto a^{-4}$ . In this stage, neutrinos are in thermal equilibrium, and their momentum follows the standard Fermi-Dirac distribution

$$n\left(p\right)\mathrm{d}p=\frac{4\pi g_{\nu }}{{\left(2\pi \mathrm{\hslash }c\right)}^3}\frac{p^2\mathrm{d}p}{e^{\frac{p}{k_B\mathit{T\nu }}}+1},$$ E13

where $n\left(p\right)$ is the number of cosmic neutrinos with momentum between p and $p+\mathit{dp}$ , $g_{\nu }$ is the number of neutrino spin states, $T_{\nu }\left(z\right)$ is the neutrino temperature at redshift z and $k_B$ is the Boltzmann constant. In principle, in the momentum distribution function, the chemical potential should be also included; however, it has been shown to be negligible for cosmological neutrinos [40]. The temperature of the cosmic neutrino background and the one from CMB are related by $T_{\nu }\left(z=0\right)={\left(\frac{4}{11}\right)}^{1/3}{T}_{\gamma }\left(z=0\right)$ [41], such that the temperature of the neutrino background at certain redshift z is given by $T_{\nu }\left(z\right)\cong 1.95\left(1+z\right)\mathrm{K}$ . Then, when the average momentum of neutrinos drops below a certain mass, they become non-relativistic, and their energy density drops like $\propto a^{-3}$ , behaving like baryons and cold dark matter. Figure 2 shows the evolution of the massive neutrino density, normalised to the today’s critical density, as a function of the scale factor from its early stage to the late universe.

Figure 2.

Evolution of the massive neutrino density, normalised with respect today’s critical density, as a function of the scale factor a, for a neutrino mass m_𝜈  = 0.1 eV. The pink solid line shows the correct solution for ${\rho }_{\nu }\left(a\right)$ , while the green and blue dashed lines show the limit at early and late time, respectively. The early time limit is proportional to a ⁴, underlining that at early times neutrinos are relativistic and behave like radiation. The late time limit, instead, is proportional to a ³, following the same evolution of matter. The short vertical dashed line (black) marks the non-relativistic transition for neutrinos of this given mass. The figure was computed by using the CAMB code.

After neutrinos with mass $m_{\nu }$ decouple from the rest of the plasma at redshift $z_{\mathit{nr}}$ , as shown by Eq. (14)

$$1+z_{\mathit{nr}}\left(m_{\nu }\right)\simeq 1890\left(\frac{m_{\nu }}{1\mathrm{eV}}\right),$$ E14

the number density per flavour is fixed by the temperature, so that the universe is currently filled by a relic neutrino background, uniformly distributed, with a density of 113 part/cm³ per species and an average temperature of 1.95 K. As neutrinos are non-relativistic particles at late times, they contribute to the total matter density of the universe ${\mathrm{\Omega }}_M$ , so that ${\mathrm{\Omega }}_M={\mathrm{\Omega }}_{\mathrm{C}\mathit{DM}}+{\mathrm{\Omega }}_b+{\mathrm{\Omega }}_{\nu }$ , where ${\mathrm{\Omega }}_{\mathrm{C}\mathit{DM}}$ , ${\mathrm{\Omega }}_b$ and ${\mathrm{\Omega }}_{\nu }$ are the dimensionless density parameters for CDM, baryons and neutrinos, respectively. The density background is affected by massive neutrinos such that a perturbation in the density field is well described by [39, 42] as follows:

$${\delta }_m=\left(1-f_{\nu }\right){\delta }_{\mathit{CDM}}+f_{\nu }{\delta }_{\nu },\text{where}{f}_{\nu }\equiv \frac{{\mathrm{\Omega }}_{\nu }}{{\mathrm{\Omega }}_{\mathit{CDM}}+{\mathrm{\Omega }}_{\nu }},$$ E15

${\delta }_{\nu }$ being the neutrino perturbations and ${\mathrm{\Omega }}_{\nu }$ the density contribution related to massive neutrinos that can be expressed in terms of the total neutrino mass, $m_{\nu }\equiv {\sum }_i{m}_{{\nu }_i}$ , as follows:

$${\mathrm{\Omega }}_{\nu }=\frac{{\mathrm{\Sigma }}_i{m}_{{\nu }_i}}{93.14h^2\mathit{eV}}.$$ E16

Currently, several observations provide limits on the total neutrino mass under the assumption of standard GR [43]. Depending of their mass, neutrinos can affect different quantities such as the matter-radiation equality and at the same time imprint features in cosmological observables like the clustering, the matter power spectrum, the halo mass function and the redshift-space distortions [18, 44, 45, 46, 47].

3. Halo mass function and clustering analysis

Considering the impact of MG and massive neutrinos in the clustering, a powerful cosmological test to discriminate among these scenarios is provided by the redshift-space distortions (RSD), that is, the shift in the position of the tracers due to their peculiar motions. For this purpose, cosmological simulations have become a powerful tool for testing theoretical predictions and to lead observational projects. In this context, the formation and evolution of cosmic structures can be understood as a dynamical system of many particles, which trace the underlying mass distribution in a certain cosmological model. The N-body simulations, methods and algorithms have progressed continuously, achieving a high resolution to resolve finer structures with millions of particles, reducing the gap between theory and observations. For a detailed description on fundamentals of cosmological simulations, see, e.g., [48, 49, 50, 51].

Since the formation and evolution of cosmic structures is based on the growth of small fluctuations in the density field, it is expected that the amplitude of these initial perturbations have the correct value at late times to match the observed clustering today. An analytical development based on perturbation theory makes possible to follow the growth of structures to a certain extent using the linear approximation, being valid as long as $\mathit{\delta \rho }\ll \rho$ . Nevertheless, these calculations are limited and cannot be extrapolated to explain completely the observational data; they break down on a scale where the density contrast $\delta \sim 1$ . Moreover, beyond the linear regime, the observed structures have a density contrast in a wide range from cosmic voids with $\delta \sim -1$ to $\delta \sim {10}^6$ and larger. It makes necessary a more elaborated description of the perturbations in the non-linear regime, which can be achieved using higher-order perturbation theory or numerical simulations [48].

In this section we show a complementary analysis to the one performed in [52] in order to investigate the clustering in the context of modified gravity with massive neutrinos. We used a subset of the DUSTGRAIN-pathfinder runs [36], which implement the Hu-Sawicki f(R) model including massive neutrinos and whose cosmological parameters are consistent with Planck 2015 constraints [53].

3.1 Halo mass function

As CDM haloes form from collapsing regions that detach from the background density field, their abundance can be related to the volume fraction of a Gaussian density smoothed on a radius R above a critical collapse threshold ${\delta }_c$ [54]. The comoving number density of the haloes is strictly related to underlining cosmological model, such that within a mass interval $\left[M,M+\mathit{dM}\right]$ , the halo mass function is given by

$$\frac{\mathrm{d}n\left(M,z\right)}{\mathrm{d}M}=f\left(\sigma \left(M,z\right)\right)\frac{\overline{\rho }}{M}\frac{\mathrm{d}ln\sigma {\left(M,z\right)}^{-1}}{\mathrm{d}M},$$ E17

where $f\left(\sigma \right)$ is the multiplicity function, $\sigma$ the RMS variance of the linear density field smoothed on scale R(M) and $\overline{\rho }$ the mean matter density. The product $f\left(\sigma \left(M,z\right)\right)\overline{\rho }$ quantifies the amount of mass contained in fluctuations of typical mass $M=\frac{4}{3}\pi R^3\overline{\rho }$ . The simplest argument to compute analytically the multiplicity function $f\left(\sigma \right)$ comes from the spherical collapse theory, following the [54] formalism, such that a perturbation is supposed to collapse when it reaches the threshold ${\delta }_c\simeq 1.68$ , by assuming that the probability distribution for a perturbation on a scale M is a Gaussian function with variance ${\sigma }_M^2$ , resulting in

$$f\left(\sigma \left(M\right)\right)=\sqrt{\frac{2}{\pi }}\frac{{\delta }_c}{\sigma \left(M\right)}exp\left(-\frac{{\delta }_c^2}{2{\sigma }^2\left(M\right)}\right).$$ E18

Another approach to determine $f\left(\sigma \right)$ is given by accurate fitting functions, like the proposed by [55], which extends phenomenologically the results of [54]. For [55] the function $f\left(\sigma \right)$ is expected to be universal to the changes in redshift and cosmology and is parameterized as follows:

$$f\left(\sigma \right)=A\left[{\left(\frac{\sigma }{b}\right)}^{-a}+{\sigma }^{-c}\right]exp\left(-d/{\sigma }^2\right),$$ E19

where A is an amplitude of the mass function and a, b, c and d are free parameters that depend on halo definition. The variance ${\sigma }^2$ is usually given by

$${\sigma }^2\left(R\left(M\right),z\right)=\int \frac{P\left(k,z\right)}{2{\pi }^2}{W}^2\left(\mathit{kR}\left(M\right)\right)k^2\mathrm{d}k,$$ E20

where $P\left(k\right)$ is the linear matter power spectrum as a function of the wave number k and W is the Fourier transform of the real-space top-hat window function of radius R. A fundamental feature of the mass function is that it decreases monotonically with increasing masses; furthermore, its dependency on cosmology is encoded in the variance ${\sigma }^2$ , as shown by the integrand of Eq. (20). From the point of view of N-body simulations, an approach to compute the mass function is given straightforward from Eq. (17), by counting the number of haloes $N_h$ above a certain mass threshold $M_{\mathit{min}}$ in a comoving volume V such as

$$N_h=A{\int }_{z_{\mathit{min}}}^{z_{\mathit{max}}}{\int }_{M_{\mathit{min}}}^{\infty }\frac{\mathrm{d}V}{\mathrm{d}z}n\left(M,z\right)\mathrm{d}M\mathrm{d}z,$$ E21

where A is the area, $z_{\mathit{min}}$ and $z_{\mathit{max}}$ are the redshift boundaries and $\mathit{dV}/\mathit{dz}$ is the comoving volume element.

Figure 3 shows the mass function of CDM haloes measured for all models of the DUSTGRAIN-pathfinder runs at six different redshifts $z=\mathrm{0,0.5,1,1.4,1.6,2}$ . Each panel contains the mass function, per each model as labelled, to track its evolution in redshift. As reference, the black dashed line represents the theoretical expectation by [55] for a flat ΛCDM model. As expected, massive haloes are less abundant with respect to smaller ones in a fixed comoving volume. The mass function decreases with redshift, since at earlier times the density field is smoother than at late times. The plot is logarithmic, meaning that the number density of large mass haloes falls off by several orders of magnitude over the range of redshifts shown. The f(R) models both with and without neutrinos reproduce in very well agreement this pattern, but only at really high masses, significant differences appear.

Figure 3.

The mass function of CDM haloes per each model of the DUSTGRAIN-pathfinder project at six different redshifts $z=\mathrm{0,0.5,1,1.4,1.6,2}$ as labelled. Each panel corresponds to one model as labelled, and the dashed line represents the theoretical expectation by [55], assuming a flat ΛCDM model.

Figure 4 compares the halo mass functions of the different DUSTGRAIN-pathfinder simulations, computed at $z=0$ (left column), $z=1$ (central column) and $z=1.6$ (right column). The lower panels show the percentage difference with respect to the ΛCDM model. It is possible to see that the effect of f(R) and massive neutrinos on the dynamical evolution of the matter density field results in different halo formation epochs and different number density of collapsed systems. In particular, the $\mathit{fR}4$ model (blue) is the most deviated model from the standard scenario, whereas the $\mathit{fR}6$ , $\mathit{fR}6\_0.06\mathit{eV}$ and $\mathit{fR}6\_0.1\mathit{eV}$ models mimic the ΛCDM behaviour over a wide range of masses.

Figure 4.

The mass function of CDM haloes for all the models of the DUSTGRAIN-pathfinder at three different redshifts: $z=0$ (left column), $z=1$ (central column) and $z=1.6$ (right column). The lower panels show the percentage difference with respect to the ΛCDM model. As in Figure 3 the dashed line represents the theoretical prediction by [55].

3.2 Clustering analysis

To quantify the halo clustering, we used the two-point correlation function (2PCF) that can be defined as the joint probability of finding a pair of objects at certain spatial separation given a tracer distribution. To measure the full 2PCF denoted as $\xi \left(r,\mu \right)$ , we used the Landy-Szalay estimator given by [56]:

$$\widehat{\xi }\left(r,\mu \right)=\frac{\mathit{DD}\left(r,\mu \right)-2\mathit{DR}\left(r,\mu \right)+\mathit{RR}\left(r,\mu \right)}{\mathit{RR}\left(r,\mu \right)},$$ E22

where $\mu$ is the cosine of the angle between the line of sight and the comoving halo pair separation, $r$ , and $\mathit{DD},\mathit{RR}$ and DR represent the normalised number of data-data, random-random and data-random pairs, respectively. This estimator is almost unbiased with minimum variance; this is the reason why it is preferred over the other estimators regarding the clustering measurements. Since the possible deviations from GR are more evident on small scales, we consider an intermediate non-linear range from 1  $h^{-1}\mathrm{Mpc}$ to 50  $h^{-1}\mathrm{Mpc}$ and random samples ten times larger than the halo ones. Then, in order to examine how significant is the RSD correction, it is convenient to expand the 2D 2PCF in the orthonormal basis of the Legendre polynomials $L_l\left(\mu \right)$ [57] such that

$$\xi \left(s,\mu \right)\equiv {\xi }_0\left(s\right)L_0\left(\mu \right)+{\xi }_2\left(s\right)L_2\left(\mu \right)+{\xi }_4\left(s\right)L_4\left(\mu \right),$$ E23

where each coefficient corresponds to the lth multipole moment:

$${\xi }_l\left(r\right)=\frac{2l+1}{2}{\int }_{-1}^{+1}\mathit{d\mu }\xi \left(r,\mu \right)L_l\left(\mu \right).$$ E24

Figure 5 shows the RSD effects on the iso-correlation curves of the 2D two-point correlation function (2PCF) in the plane ( $r_{\perp },r_{\parallel }$ ), where $r_{\perp }$ and $r_{\parallel }$ coordinates are, respectively, the perpendicular and parallel components along the line of sight of the observer. The 2PCF has been computed for the ΛCDM catalogue of the DUSTGRAIN-pathfinder simulations, in real space (left panel), and for the corresponding sample in redshift space (right panel). The contours are drawn at the iso-correlation levels $\xi \left(s_{\perp },s_{\parallel }\right)=\mathrm{0.3,0.5,1.0,1.4,2.2,3.6,7.2,21.6}$ . In real space the correlation function is undistorted, describing circular curves in this plane. In redshift space, the effect caused by the RSD is clearly visible on small scales, where the 2PCF is stretched in the direction of $r_{\parallel }$ (Fingers-of-God effect), and in the infall effect on large scales, the contours are squashed along the perpendicular direction (Kaiser effect).

Figure 5.

Contours of the 2D two-point correlation function (2PCF) in the plane ( $r_{\perp },r_{\parallel }$ ), the $r_{\perp }$ and $r_{\parallel }$ coordinates are, respectively, the perpendicular and parallel components along the line of sight of the observer. The 2PCF has been computed from an N-body simulation, in real space (left panel), and for the corresponding sample in redshift space (right panel). The contours are draw at the iso-correlation levels $\xi \left(s_{\perp },s_{\parallel }\right)=\mathrm{0.3,0.5,1.0,1.4,2.2,3.6,7.2,21.6}$ as indicated by the colour bar. Both the effect of redshift space distortions on small scales (Finger-of-God effect) and the infall effect at large scales (Kaiser effect) are clearly visible in the left panel. On small scales the 2PCF is stretched in the direction of $r_{\parallel }$ , and on large scales the contours are squashed along the perpendicular direction.

The symmetry of the 2PCF in real space means that the full clustering signal is encoded in the monopole moment ${\xi }_0\left(r\right)$ , while the rest of the multipole moments are statistically equal to zero. Figure 6 shows the monopole moment, ${\xi }_0\left(r\right)$ , of the CDM haloes for all models considered in the DUSTGRAIN-pathfinder project, at three different redshifts $z=\mathrm{0,1,2}$ . Subpanels show the percentage difference between MG models [f(R) with and without massive neutrinos] and the ΛCDM model. The monopole moments of the 2PCF of the $\mathit{fR}4$ and $\mathit{fR}4\_0.3\mathit{eV}$ models are the ones that deviate most from ΛCDM at low redshift. This behaviour is also present in the models $\mathit{fR}5$ , $\mathit{fR}5\_0.15\mathit{eV}$ and $\mathit{fR}5\_0.1\mathit{eV}$ , but it is less significant. In general, it is observed that massive neutrinos increase the clustering signal for all models, especially at high redshift, the $\mathit{fR}6$ , $\mathit{fR}6\_0.06\mathit{eV}$ and $\mathit{fR}6\_0.1\mathit{eV}$ models being degenerated with respect to ΛCDM. The quadrupole and hexadecapole moments are consistent with zero at 1 $\sigma$ error bar.

Figure 6.

The real-space 2PCF $r^2{\xi }_0$ of CDM haloes for all the models of the DUSTGRAIN-pathfinder project at three different redshifts: $z=0$ (left column), $z=1$ (central column) and $z=2$ (right column). From top to bottom, the panels show the $\mathit{fR}4$ , $\mathit{fR}5$ and $\mathit{fR}6$ models, respectively, compared to the results of the ΛCDM model. The error bars, shown only for the ΛCDM model for clarity reasons, are the diagonal values of the bootstrap covariance matrices used for the statistical analysis. Percentage differences between f(R), $f\left(R\right)+m_{\nu }$ and ΛCDM predictions are in the subpanels, while the shaded regions represent the deviation at 1 $\sigma$ confidence level.

Another important feature to take into account in the clustering analysis is related to the bias that is introduced when the ΛCDM model is wrongly assumed to predict the DM clustering of a f(R) universe with massive neutrinos [58]. This effective halo bias, $⟨b⟩$ , allows to characterise the relation between the halo clustering and the underlying mass distribution. By using the theoretical effective bias proposed by [59] and its corresponding mass function, it is possible to disentangle the degeneracy with respect to ${\sigma }_8$ . The level of agreement with the measurements obtained from the f(R) and $f\left(R\right)+m_{\nu }$ scenarios give us a better understating of the discrepancy with the ΛCDM ones. To compute the theoretical mass function, we consider the linear CDM + baryon power spectrum as have been stated in the CDM prescription [60] and replacing ${\rho }_m$ with ${\rho }_{\mathrm{CDM}}$ [42]. These quantities can be obtained with CAMB or another Boltzmann code; the linear power-spectrum for CDM, $P_{\mathit{lin}}^{\mathrm{CDM}}\left(k\right)=T_{\mathrm{CDM}}^2/T_m^2{P}_{\mathit{lin}}^m\left(k\right)$ , is expressed in terms of $P_{\mathit{lin}}^{\mathrm{CDM}+b}\left(k\right)$ and $P_{\mathit{lin}}^m\left(k\right)$ , with $T_{\mathrm{CDM}}\left(k\right)$ and $T_b\left(k\right)$ being the transfer functions. It implies, as shown by [61], that the effect of neutrinos on the cluster abundance is well captured by rescaling the smoothed density field such that

$${\sigma }^2\to {\sigma }_{\mathrm{CDM}}^2\left(z\right)=\int \frac{P_{\mathrm{CDM}}\left(k,z\right)}{2{\pi }^2}{W}^2\left(\mathit{kR}\right)k^2\mathrm{d}k$$ E25

with the CDM power spectrum obtained by rescaling the total matter power spectrum with the corresponding transfer functions, $T_{\mathit{CDM}}$ and $T_b$ , weighted by the density of each species so that

$$P_{\mathrm{CDM}}\left(k,z\right)=P_{\mathrm{m}}\left(k,z\right){\left(\frac{{\mathrm{\Omega }}_{\mathrm{CDM}}{T}_{\mathrm{CDM}}\left(k,z\right)+{\mathrm{\Omega }}_{\mathrm{b}}{T}_{\mathrm{b}}\left(k,z\right)}{T_{\mathrm{m}}\left(k,z\right)\left({\mathrm{\Omega }}_{\mathrm{CDM}}+{\mathrm{\Omega }}_{\mathrm{b}}\right)}\right)}^2.$$ E26

The impact on the bias when the CDM prescription is not considered can be appreciated in Figure 7. In most cases this correction is small with the exception of the $\mathit{fR}4\_0.3\mathrm{eV}$ model. A detailed discussion on these results can be found in [52].

Figure 7.

The coloured solid lines represent the apparent effective halo bias, $⟨b⟩$ , as a function of redshift, averaged in the range 10 $h^{-1}\mathit{Mpc}$ $<r<$ 50 $h^{-1}\mathit{Mpc}$ . Black lines show the theoretical ΛCDM effective bias predicted by [59] (dashed), normalised to the ${\sigma }_8$ values of each DUSTGRAIN-pathfinder simulation, while the cyan-shaded areas show a 10% error. The upper set of panels shows the results, considering the total power spectrum, while the lower set of panels shows the results when the CDM+baryon power spectrum is used instead.

3.3 Modelling the redshift-space distortions

In a realistic case, spectroscopic surveys observe a combination of density and velocity fields in redshift space. The observed redshift is a combination of cosmological effects plus an additional term caused by the peculiar motions along the line of sight of the observer. This combination makes the redshift-space catalogues appear distorted with respect to the real-space ones, and they can be reproduced from N-body simulations since the positions and velocities are known. Currently, the modelling of the redshift-space distortions provide a powerful tool to test the gravity theory by exploring the spatial statistics encoded in the 2PCF, which is anisotropic due to the dynamic distortions.

We consider the first two even multipoles of the 2PCF, ${\xi }_0$ and ${\xi }_2$ , taking into account that odd multipole moments vanish by symmetry. The analysis consisted of modelling (i) the individual multipole moments and (ii) both multipole moments simultaneously. To derive cosmological constraints from the clustering signal and quantify the effects of f(R) gravity and massive neutrinos on RSD, we performed a Bayesian analysis to set constraints on the linear growth rate $f\left({\mathrm{\Omega }}_m\right)$ and the linear bias. All numerical tasks were performed with the CosmoBolognaLib¹ [62].

The Kaiser formula is a good description of the RSD only at very large scales, where non-linear effects can be neglected, but it does not describe accurately the non-linear regime. Thus, with the aim of extracting information from the RSD signal at non-linear regime and considering the increasing precision of recent and upcoming surveys, many more approaches have been proposed. There is a vast literature that shows the efforts to model the RSD beyond the linear Kaiser model [63, 64, 65, 66], some of them making use of a phenomenological description of the velocity field and others, instead, taking into account higher orders in perturbation theory since, in principle, there is no reason to stop at linear order. Other approaches do a combination of both frameworks. A simple alternative to model the redshift-space 2PCF at small scales consists of extending the Kaiser formula, by adding a phenomenological damping factor that plays the role of a pairwise velocity distribution. It can account for both linear and non-linear dynamics. Therefore, to construct the likelihood, we consider this model sometimes called dispersion model [67], which introduces a damping function to describe the distortions in the clustering at small scales (Fingers-of-God). This model is enough accurate to quantify the relative differences between f(R) models with massive neutrinos and ΛCDM. For the Bayesian analysis, the dispersion model is fully described by three parameters, $f{\sigma }_8$ , $b{\sigma }_8$ and ${\mathrm{\Sigma }}_S$ , that we constrain by minimising numerically the negative log-likelihood.

Figure 8 shows the normalised covariance matrices $\left(C_{i,j}/\sqrt{C_{i,i}{C}_{j,j}}\right)$ of the redshift-space monopole and quadrupole moments of CDM haloes with bootstrap errors resampling at three different redshifts $z=\mathrm{0.5,1.0}$ and 1.6. As it can be appreciated, the covariance matrices represent how the scatter propagates into the likelihood and on the final posterior probabilities of the parameters. Then, to assess the posterior distributions of the three model parameters, we perform a MCMC analysis. The fitting analysis is limited to the scale range $10\le r\left[\mathrm{Mpc}{\mathrm{h}}^{-1}\right]\le 50$ , assuming flat priors in the ranges $0\le f{\sigma }_8\le 2$ , $0\le b{\sigma }_8\le 3$ and $0\le {\mathrm{\Sigma }}_S\le 2$ . Figure 9 shows the $f{\sigma }_8$ - $b{\sigma }_8$ posterior constraints obtained from the MCMC analysis of ${\xi }_0$ , ${\xi }_2$ and ${\xi }_0+{\xi }_2$ for each mock catalogue of the DUSTGRAIN-pathfinder. The figure shows the constraints for all models considered in this work at $z=1.0$ using the monopole (orange), quadrupole (green) and monopole plus quadrupole (blue), while the intersection regions correspond to the joint analysis of the multipoles. Thus, the joint analysis of the redshift-space monopole and quadrupole is able to break the degeneracy between $f{\sigma }_8$ and $b{\sigma }_8$ even in the presence of massive neutrinos. In Figure 10, we show the theoretical behaviour of the linear distortion parameter and the growth factor as a function of the redshift for each family of models, assuming a flat ΛCDM model.

Figure 8.

Covariance matrices from the analysis of the redshift-space monopole and quadrupole moments of CDM haloes with Bootstrap errors. The models correspond to ΛCDM (upper panels), fR4 (central panels) and fR5 (bottom panels), at three different redshifts $z=\mathrm{0.5,1.0}$ and 1.6 from left to right.

Figure 9.

Contours at $1-2\sigma$ confidence level of the $f{\sigma }_8-b{\sigma }_8$ posterior distributions, obtained from the MCMC analysis in redshift space for 2PCF multipoles of CDM haloes. The contours correspond to $z=1.0$ with the monopole in orange and the quadrupole in green.

Figure 10.

Theoretical expectation of the linear distortion parameter $\beta$ (upper panel) and the growth rate $f{\sigma }_8$ (lower panel) as a function of the redshift for each family of mocks from the DUSTGRAIN-pathfinder runs, assuming a flat ΛCDM model.

From the $f{\sigma }_8$ - $b{\sigma }_8$ posterior contours, at $1-2\sigma$ confidence levels, the results suggest that the clustering information encoded in the two first non-null multipole moments of the 2PCF can discriminate the alternative MG models considered in this work at $z\gtrsim 1$ . At low redshifts the f(R) models studied are statistically indistinguishable from ΛCDM, and further studies are required to break this degeneracy.

4. Conclusions

In this chapter we have introduced the theoretical framework of modern cosmology in present massive neutrinos. We emphasize on the structure formation and on the statistical description of the density field as well as the measurements of galaxy clustering and discuss the redshift-space distortions and the differences between clustering in real and redshift space, considering that in the recent years, the spatial distribution of matter on cosmological scales has become one of the most efficient probes to investigate the properties of the universe, such as test gravity theories on large scales, to explore the dark sector and the origin of the accelerated expansion of the universe as well as a probe to constrain alternative cosmological models.

In the context of models based on modified gravity and massive neutrino cosmologies, we investigated the spatial properties of the large-scale structure by exploiting the DUSTGRAIN-pathfinder simulations that follow, simultaneously, the effects of f(R) gravity and massive neutrinos. These are two of the most interesting scenarios that have been recently explored to account for possible observational deviations from the standard ΛCDM model. In particular, we studied whether redshift-space distortions in the 2PCF multipole moments can be effective, breaking the cosmic degeneracy between these two effects. We analysed the redshift-space distortions in the clustering of dark matter haloes at different redshifts, focusing on the monopole and quadrupole moments of the two-point correlation function, both in real and redshift space. The deviations with respect to ΛCDM model have been quantified in terms of the linear growth rate parameter. We found that multipole moments of the 2PCF from redshift-space distortions provide a useful probe to discriminate between ΛCDM and modified gravity models, especially at high redshifts ( $z\gtrsim 1$ ), even in the presence of massive neutrinos. The linear growth rate constraints that we obtain from all the analysed $f\left(R\right)+m_{\nu }$ mock catalogues are statistically distinguishable from ΛCDM predictions at high redshifts.

Acknowledgments

We thank Lauro Moscardini, Federico Marulli and Alfonso Veropalumbo for the continuous development of the CosmoBolognaLib and their suggestions during this project. We also thank Carlo Giocoli and Marco Baldi for the crucial work performing the DUSTGRAIN-pathfinder runs.

Conflict of interest

The authors declare no conflict of interest.