Investigating the degeneracy between modified gravity and massive neutrinos with redshift-space distortions

TL;DR

This paper explores how redshift-space distortions can help distinguish between modified gravity and massive neutrinos, by modeling non-linear effects to break their degeneracy in the matter power spectrum.

Contribution

It introduces an extended perturbation theory model for non-linear RSD that accounts for both Hu-Sawicki f(R) gravity and massive neutrinos to improve discrimination.

Findings

Degeneracy can be broken if the linear growth rate is accurately recovered from RSD.

Non-linear modeling is essential for distinguishing modified gravity from neutrino effects.

The extended model improves the interpretation of RSD signals in cosmological tests.

Abstract

There is a well known degeneracy between the enhancement of the growth of large-scale structure produced by modified gravity models and the suppression due to the free-streaming of massive neutrinos at late times. This makes the matter power-spectrum alone a poor probe to distinguish between modified gravity and the concordance $\Lambda$CDM model when neutrino masses are not strongly constrained. In this work, we investigate the potential of using redshift-space distortions (RSD) to break this degeneracy when the modification to gravity is scale-dependent in the form of Hu-Sawicki $f(R)$. We find that if the linear growth rate can be recovered from the RSD signal, the degeneracy can be broken at the level of the dark matter field. However, this requires accurate modelling of the non-linearities in the RSD signal, and we here present an extension of the standard perturbation theory-based model for non-linear RSD that includes both Hu-Sawicki $f(R)$ modified gravity and massive neutrinos.

Tables

Table 1. Table 1 : Survey parameters for a DESI-like survey computed from the information for emission line galaxies (ELGs) in Table V of [ 50 ] . These parameters are used in the computation of the covariance matrices for fitting σ v subscript 𝜎 𝑣 \sigma_{v} in MG-Copter in the study of the degeneracy in Section  4 .

$z$	$V_s$ (Gpc3/$h^3$)	${\overline{n}}_g$ ($h^3$/Mpc3)
0.5	3.40	2.95$\times {10}^{-4}$
1.0	7.68	5.23$\times {10}^{-4}$
1.5	10.14	1.71$\times {10}^{-4}$

Equations

$$a\frac{\partial \delta (\vec{k})}{\partial a}+\theta (\vec{k})$$

$$a\frac{\partial \theta (\vec{k})}{\partial a}+\left(2+\frac{a{H}^{\prime }}{H}\right)\theta (\vec{k})$$

$$=-\frac{1}{2}\int \frac{d^3{\vec{k}}_1{d}^3{\vec{k}}_2}{(2\pi {)}^3}{\delta }_D\left(\vec{k}-{\vec{k}}_1-{\vec{k}}_2\right)\beta ({\vec{k}}_1,{\vec{k}}_2)\theta ({\vec{k}}_1)\theta ({\vec{k}}_2),$$

$$\alpha ({\vec{k}}_1,{\vec{k}}_2)$$

$$\beta ({\vec{k}}_1,{\vec{k}}_2)$$

$$-{\left(\frac{k}{aH}\right)}^2\Phi (\vec{k})=\frac{3{\Omega }_{\mathrm{m}}(a)}{2}\delta (\vec{k}),$$

$${\delta }_n(\vec{k},a)=\int d^3{\vec{k}}_1\dots d^3{\vec{k}}_n{\delta }_D(\vec{k}-{\vec{k}}_{1\dots n})F_n({\vec{k}}_1,\dots ,{\vec{k}}_n,a)\Delta ({\vec{k}}_1)\dots \Delta ({\vec{k}}_n),$$

$${\theta }_n(\vec{k},a)=\int d^3{\vec{k}}_1\dots d^3{\vec{k}}_n{\delta }_D(\vec{k}-{\vec{k}}_{1\dots n})G_n({\vec{k}}_1,\dots ,{\vec{k}}_n,a)\Delta ({\vec{k}}_1)\dots \Delta ({\vec{k}}_n),$$

$$\hat{\mathcal{L}}\left[\begin{array}{c}{F}_n({\vec{k}}_1,\dots ,{\vec{k}}_n)\\ G_n({\vec{k}}_1,\dots ,{\vec{k}}_n)\end{array}\right]=\sum _{j=1}^{n-1}\left[\begin{array}{c}-\alpha ({\vec{k}}_{1\dots j},{\vec{k}}_{j+1\dots n})G_j({\vec{k}}_1,\dots ,{\vec{k}}_j)F_{n-j}({\vec{k}}_{j+1},\dots ,{\vec{k}}_n)\\ -\frac{1}{2}\beta ({\vec{k}}_{1\dots j},{\vec{k}}_{j+1\dots n})G_j({\vec{k}}_1,\dots ,{\vec{k}}_j)G_{n-j}({\vec{k}}_{j+1},\dots ,{\vec{k}}_n)\end{array}\right],$$

$$\hat{\mathcal{L}}=\left[\begin{array}{cc}a\frac{d}{da}& 1\\ \frac{3{\Omega }_{\mathrm{m}}}{2}& a\frac{d}{da}+\left(2+\frac{a{H}^{\prime }}{H}\right)\end{array}\right].$$

$$P_{ij}^{1-\mathrm{loop}}(k)=P_{ij}^{\mathrm{L}}(k)+P_{ij}^{13}(k)+P_{ij}^{22}(k),$$

$$⟨x_2(\vec{k})y_2({\vec{k}}^{\prime })⟩$$

$$⟨x_1(\vec{k})y_3({\vec{k}}^{\prime })+x_3(\vec{k})y_1({\vec{k}}^{\prime })⟩$$

$$P_{\delta \delta }^{22}(k)=$$

$$P_{\delta \theta }^{22}(k)=$$

$$P_{\theta \theta }^{22}(k)=$$

$$P_{\delta \delta }^{13}(k)=$$

$$P_{\delta \theta }^{13}(k)=$$

$$+3\frac{k^3}{(2\pi {)}^2}{G}_1(k)P_0(k){\int }_0^{\infty }{r}^2{P}_0(kr)F_3(k,r,x)dr,$$

$$P_{\theta \theta }^{13}(k)=$$

$$P^{(s)}(k,\mu )=D_{\mathrm{FoG}}\left[k\mu {\sigma }_v\right]\left\{P_{\delta \delta }(k)-2{\mu }^2{P}_{\delta \theta }(k)+{\mu }^4{P}_{\theta \theta }(k)+A(k,\mu )+B(k,\mu )\right\},$$

$$A(k,\mu )$$

$$B(k,\mu )$$

$$⟨\theta ({\vec{k}}_1)\left\{\delta ({\vec{k}}_2)-\frac{k_{2z}^2}{k_2^2}\theta ({\vec{k}}_2)\right\}\left\{\delta ({\vec{k}}_3)-\frac{k_{3z}^2}{k_3^2}\theta ({\vec{k}}_3)\right\}⟩$$

$$=(2\pi {)}^3{\delta }_D({\vec{k}}_1+{\vec{k}}_2+{\vec{k}}_3)B_{\sigma }({\vec{k}}_1,{\vec{k}}_2,{\vec{k}}_3),$$

$$F(\vec{p})=\frac{p_z}{p^2}\left\{P_{\delta \theta }(p)-\frac{p_z^2}{p^2}{P}_{\theta \theta }(p)\right\}.$$

$$D_{\mathrm{FoG}}\left[k\mu {\sigma }_v\right]=\exp \left(-k^2{\mu }^2{\sigma }_v^2\right).$$

$$-2\ln \mathcal{L}=\sum _n\sum _{l,l^{\prime }=0,2}\left(P_{l,\text{Copter}}^s(k_n)-P_{l,\text{COLA}}^s(k_n)\right){\mathrm{Cov}}_{l,l^{\prime }}^{-1}(k_n)\left(P_{l^{\prime },\text{Copter}}^s(k_n)-P_{l^{\prime },\text{COLA}}^s(k_n)\right)$$

$$-{\left(\frac{k}{aH}\right)}^2\Phi (\vec{k})=\frac{3{\Omega }_{\mathrm{m}}(a)}{2}\delta (\vec{k})\mu (k,a)+S(\vec{k}),$$

$$S(\vec{k})=$$

$$+\int \frac{d^3{\vec{k}}_1{d}^3{\vec{k}}_2{d}^3{\vec{k}}_3}{(2\pi {)}^3}{\delta }_D(\vec{k}-{\vec{k}}_{123}){\gamma }_3(\vec{k},{\vec{k}}_1,{\vec{k}}_2,{\vec{k}}_3,a)\Delta ({\vec{k}}_1)\Delta ({\vec{k}}_2)\Delta ({\vec{k}}_3).$$

$$\hat{\mathcal{L}}\left[\begin{array}{c}{F}_n({\vec{k}}_1,\dots ,{\vec{k}}_n)\\ G_n({\vec{k}}_1,\dots ,{\vec{k}}_n)\end{array}\right]=\sum _{j=1}^{n-1}\left[\begin{array}{c}-\alpha ({\vec{k}}_{1\dots j},{\vec{k}}_{j+1\dots n})G_j({\vec{k}}_1,\dots ,{\vec{k}}_j)F_{n-j}({\vec{k}}_{j+1},\dots ,{\vec{k}}_n)\\ -\frac{1}{2}\beta ({\vec{k}}_{1\dots j},{\vec{k}}_{j+1\dots n})G_j({\vec{k}}_1,\dots ,{\vec{k}}_j)G_{n-j}({\vec{k}}_{j+1},\dots ,{\vec{k}}_n)-N_n(\vec{k},{\vec{k}}_1,\dots ,{\vec{k}}_n)\end{array}\right],$$

$$\hat{\mathcal{L}}=\left[\begin{array}{cc}a\frac{d}{da}& 1\\ \frac{3{\Omega }_{\mathrm{m}}}{2}\mu (k,a)& a\frac{d}{da}+\left(2+\frac{a{H}^{\prime }}{H}\right)\end{array}\right],$$

$$N_2=$$

$$N_3=$$

$$+{\gamma }_3(\vec{k},{\vec{k}}_1,{\vec{k}}_2,{\vec{k}}_3)F_1({\vec{k}}_1)F_1({\vec{k}}_2)F_1({\vec{k}}_3).$$

$$\mu (k,a)=$$

Full text

1 Introduction

Modifications to Einstein’s theory of General Relativity (GR) can be considered when searching for an explanation of the late-time acceleration [1, 2]. The simplest class of models, known as scalar-tensor theories [3, 4, 5], introduce a new scalar field that causes the required acceleration. This new scalar field also couples to matter, leading to a so-called fifth force. Such models typically invoke a screening mechanism to ensure that the fifth force would become negligible in high density environments, and thus not be observable in solar system tests. However, in other environments the fifth force is present and can enhance structure formation. This enhancement can be used to constrain such scalar-tensor theories of modified gravity (MG) with large-scale structure observations, and this is indeed a key goal of upcoming large-scale structure surveys such as Euclid [6], DESI [7], WFIRST [8], LSST [9], and SKA [10].

In this paper we will consider $f(R)$ gravity [11, 12], and in particular the commonly studied Hu-Sawicki variant of $f(R)$ [13]. This model has been well studied with N-body simulations as it is well understood theoretically and offers a phenomenology that is representative of a broad range of modified gravity models.

In GR with cold dark matter the growth rate of matter perturbations is scale-independent. A key signature of modified gravity is that the linear growth rate can be scale-dependent. However, a vital but often overlooked complication when searching for signatures of modified gravity in large-scale structure is the suppression of structure growth due to massive neutrinos. Neutrinos were first shown to have mass in observations of neutrino flavour oscillations [14, 15], the presence of which demand that at least two of the neutrino states are massive [16]. Even though particle physics experiments do not yet tell us the absolute mass of each of the three mass eigenstates, they do allow strong constraints to be placed on the difference in mass between the states, and these imply that at least one of the mass eigenstates has a mass $m_{\nu}\gtrsim 0.05~{}\mathrm{eV}$ [17]. As a consequence of having mass, the matter-radiation equality time will be delayed and the neutrinos will not cluster at scales below their free-streaming length $\lambda_{\rm fs}$ [18]. The delay to the time of matter-radiation equality lowers the amplitude of the matter perturbations at the start of matter-domination, and the free-streaming of massive neutrinos causes the dark matter perturbations to feel a reduced gravitational potential below $\lambda_{\rm fs}$ and thus cluster less strongly than in a model with the same value of the matter density parameter but only massless neutrinos. The combination of these two effects leads to a scale-dependent suppression of structure growth; a signature which can be used to constrain the neutrino masses if it can be measured by the previously mentioned large-scale structure surveys [19, 20, 21, 22, 23, 24, 25, 26], even in models beyond $\Lambda$CDM that affect structure growth in a scale-independent way [27].

However, with the potential for scale-dependent enhancement of structure formation from modified gravity and the scale-dependent suppression due to massive neutrinos, there is a risk of degeneracy whereby large-scale structure in a universe with a strong modification to gravity and heavy neutrinos can be difficult to distinguish from that of a universe with GR and light neutrinos [28, 29, 30, 31, 32, 33, 34]. This degrades the ability of surveys to achieve their twin goals of testing gravity and constraining the neutrino masses in any theories of gravity beyond GR. Indeed, it has been shown that the non-linear matter power spectrum [30] and halo mass function [35] in $f(R)$ models are difficult to distinguish from their equivalents in GR when the neutrino masses are allowed to vary. The DES Collaboration considers neutrino mass and extensions to GR in the same analysis [36], although they only state the resulting constraints on the MG parameters and not the neutrino masses. There are some promising signs that certain observables may be better at reducing or even breaking this degeneracy, such as higher-order weak lensing statistics [37] and weak lensing tomographic information at multiple redshifts [38]; as well as techniques that are superior at distinguishing models such as machine learning [39, 40].

A different observable that has degeneracy breaking potential is that of redshift-space distortions (RSD) [41]. RSD occur when the distances to tracers are computed using their observed redshifts without accounting for the effect of the tracers’ peculiar velocities on the redshifts which adds to the contribution from the Hubble flow. On linear scales this results in a slight squashing along the line-of-sight [42], whereas there is a strong stretching along the line-of-sight at non-linear scales commonly known as the Fingers-of-God (FoG) effect [43]. For combinations of MG parameters and neutrino masses whose enhancement and suppression of structure growth produce matter power spectra that are difficult to differentiate between, the structure growth rate can still be different in each case and allow for models to be distinguished between. It has recently been shown that growth rate information imprinted in velocity statistics in real-space can be used to break the degeneracy [44]. However, real-space velocity statistics are not directly observable. Fortunately, because of the velocity information encoded in them, RSD observations can be used to extract the linear growth rate of structure $f$. However, in order to extract $f$ and break the degeneracy, it is necessary to accurately model the non-linearities of RSD with MG and massive neutrinos.

In this paper, we extend the cosmological perturbation theory code Copter [45] to include the effects of massive neutrinos in addition to those of modified gravity allowing us to accurately model non-linear RSD in scenarios with Hu-Sawicki $f(R)$ gravity and non-zero neutrino masses. We build on MG-Copter, the modified version of Copter developed in [46], which is itself based on the approach presented in [47].

We validate this implementation against simulations using the COmoving Lagrangian Acceleration (COLA) method [48], which is a fast approximate simulation method, and then investigate whether the degeneracy between the two effects is broken by RSD at the level of the dark matter field.

The paper is organised as follows. In Section 2 we explain our implementation of modified gravity and massive neutrinos in the Standard Perturbation Theory (SPT) formalism and MG-Copter code. In Section 3, we show the results of tests validating our implementation against simulation results. In Section 4 we use our new implementation to investigate the degeneracy and then conclude in Section 5.

2 Implementation

In order to model the combined effect of modified gravity and massive neutrinos on real- and redshift-space power spectra with low computational expense, it is necessary to include both effects in a semi-analytical code such as Copter which computes large-scale structure observables using perturbation theory. For the redshift-space quantities, Copter depends on the TNS model of redshift-space distortions which is named after the authors of [49] (Taruya, Nishimichi, and Saito).

2.1 MG-Copter and the TNS model

MG-Copter [46] solves the equations of Standard Perturbation Theory (SPT) to acquire the real-space power spectra up to 1-loop order based on the approach developed by [47].

Starting from the continuity and Euler equations, assuming fluid quantities to be irrotational such that velocity field $\vec{v}$ can be expressed in terms of the velocity divergence $\theta=\left(\vec{\nabla}\cdot\vec{v}\right)/aH$, and transforming to Fourier space yields

[TABLE]

where $a=1/(1+z)$ is the scale factor, $y^{\prime}=\partial y/\partial a$ and the kernels $\alpha$ and $\beta$ are given by

[TABLE]

The Poisson equation completes the above modified continuity and Euler equations

[TABLE]

where $\Omega_{\rm m}(a)=8\pi G\rho_{\rm m}/3H^{2}$. We want the $n^{\rm th}$ order solutions of Eqs. 2.1 and 2.2 to be of the form

[TABLE]

[TABLE]

where $\vec{k}_{1\ldots n}=\vec{k}_{1}+\ldots+\vec{k}_{n}$. Inserting these forms of the solutions into Eqs. 2.1 and 2.2 yields a generalised system of equations for the $n^{\rm th}$ order kernels [47]

[TABLE]

where

[TABLE]

MG-Copter solves this system of equations to compute the kernels $F_{i}$ and $G_{i}$. The power spectra up to 1-loop are given as

[TABLE]

where the 1-loop corrections are defined by

[TABLE]

where $x$ and $y$ can be $\delta$ or $\theta$. Working these through, the final expressions for the 1-loop corrections in terms of the $z=0$ linear power spectrum $P_{0}(k)=P^{\rm L}(k,z=0)$ are, for the 22 correction,

[TABLE]

while for the 13 correction we have

[TABLE]

With the SPT real-space power spectra computed up to 1-loop order, MG-Copter can then input these to the TNS model to calculate the redshift-space power spectrum $P^{(s)}(k)$.

The TNS model for the redshift-space power spectrum $P^{(s)}$ as a function of scale $k$ and line-of-sight (LoS) angle parameter $\mu=\cos(\theta)$ is given by Eq. (18) of [49], which we reproduce here with subtle changes due to the different definition of $\theta$:

[TABLE]

where $D_{\rm FoG}$ is the Fingers-of-God damping function which we will discuss later. It is generally a function of $k$, $\mu$, and the velocity dispersion $\sigma_{v}$. The power spectra $P_{\delta\delta}(k)$, $P_{\delta\theta}(k)$, and $P_{\theta\theta}(k)$ correspond to the density auto-correlation, density-velocity divergence cross correlation, and the velocity divergence auto-correlation respectively. $A(k,\mu)$ and $B(k,\mu)$ are correction terms given by

[TABLE]

where $B_{\sigma}$ is the cross bispectrum defined by

[TABLE]

and $F(\vec{p})$ is defined as

[TABLE]

Throughout we use an exponential form for the Fingers-of-God damping factor:

[TABLE]

The velocity dispersion $\sigma_{v}$ is a free parameter and needs to be fitted to some other $P^{(s)}$ data, for example from simulations as we do here. To do this, we minimise the likelihood function

[TABLE]

for the first two multipoles. Expressions for the covariance matrix between the different multipoles $\mathrm{Cov}_{l,l^{\prime}}$ are given in Appendix C of [49]. We do not consider non-Gaussianity in this covariance but we do include the effect of shot-noise. For the validation of our implementation of massive neutrinos in MG-Copter presented in Section 3 we assume an ideal survey with survey volume $V_{\rm s}=10~{}{\rm Gpc}^{3}/h^{3}$ and galaxy number density $\bar{n}_{\rm g}=4\times 10^{-3}h^{3}/{\rm Mpc}$. For the study of the degeneracy in Section 4, we want to model a slightly more realistic scenario, so we assume a DESI-like survey with $V_{\rm s}$ and $\bar{n}_{\rm g}$ as given in Table 1 and redshift bin width $\Delta z=0.2$. These values are computed using the information for emission line galaxies (ELGs) in Table V of [50].

Thus the TNS model can be used to compute $P^{(s)}(k,\mu)$ with the input of $P_{\delta\delta}$, $P_{\delta\theta}$, $P_{\theta\theta}$ at 1-loop order from MG-Copter.

2.2 Adding modified gravity

Modified gravity models, like the $f(R)$ gravity model we consider here, have been previously added to Copter in [46], resulting in MG-Copter. The 1-loop real-space power spectra are affected by the inclusion of modified gravity in SPT, but the TNS model of Eq. (2.18) is still applicable without changes. We shall reproduce here the essentials of the implementation of modified gravity in the SPT part of MG-Copter.

The modifications to gravity can be included in the Poisson equation, which up to $3^{\rm rd}$ order becomes

[TABLE]

where $\mu(k,a)=G_{\rm eff}(k,a)/G$ is an effective Newton’s constant111Not to be confused with the line-of-sight angle parameter $\mu$, which will always be presented without arguments. and the non-linear source term $S(\vec{k})$ up to $3^{\rm rd}$ order is

[TABLE]

While the effective Newton’s constant $\mu(k,a)$ is generally responsible for the (scale-dependent) growth of linear perturbations, at the fully non-linear level modified gravity models typically include a screening mechanism that will affect the growth of non-linearities, and the $\gamma_{2}$ and $\gamma_{3}$ terms provide the leading order description of this screening in perturbation theory.

Using the same form for the $n^{\rm th}$ order solutions as in Eqs. 2.6 and 2.7, the new system of equations for the $n^{\rm th}$ order kernels is

[TABLE]

where

[TABLE]

and

[TABLE]

In this work we investigate Hu-Sawicki $f(R)$ gravity, which has a single free parameter $|f_{R0}|$. Hu-Sawicki $f(R)$ gravity produces an enhanced, scale-dependent growth of density perturbations relative to GR, but the built-in chameleon screening mechanism ensures that the modifications become negligible in high density environments (which typically coincide with small scales). Hu-Sawicki $f(R)$ only becomes active at late times, and thus the modifications to GR are negligible in the early Universe. For this theory, the extra terms in Eq. (2.25) are given as

[TABLE]

where

[TABLE]

2.3 Adding massive neutrinos

We have added support for massive neutrinos to the code MG-Copter developed in [46]. Note that massive neutrinos were also added to the original Copter code in [51] using a similar approach. In our implementation, we follow the method of [52, 53] and include massive neutrinos at the level of the linear real-space power spectra $P^{\rm L}$, $P_{\delta\theta,\mathrm{L}}=f(k)P^{\rm L}$, and $P_{\theta\theta,\mathrm{L}}=f^{2}(k)P^{\rm L}$ without modifying the higher order SPT kernels. This allows us to take $P^{\rm L}(k)$ and $f(k)$ from CAMB [54] (or MGCAMB [55, 56] for MG+$m_{\nu}$) as input to MG-Copter; note that a small modification to CAMB/MGCAMB is necessary to get scale-dependent growth rate $f(k)$ as output. This method for including massive neutrinos is general enough to handle the various hierarchies of neutrino mass eigenstates [57], but for simplicity in the results that follow we have modelled the massive neutrinos as a single massive eigenstate with mass $m_{\nu}$ and two massless eigenstates.

The free-streaming of massive neutrinos causes suppression of $P^{\rm L}(k)$ relative to the case with massless neutrinos for scales smaller than the neutrino free-streaming scale after the time at which massive neutrinos become non-relativistic – see Fig. 1 of [58] for an example. A linear approximation gives the amplitude of suppression to be $-8f_{\nu}$ where $f_{\nu}=\Omega_{\nu}/\Omega_{\rm m}$ is the fraction of total matter in massive neutrinos [18]. Scale-dependent suppression also affects $f(k)$, although the amplitude of this effect is much smaller, as can be seen in Fig. 5 of [58].

The expressions for the 1-loop power spectra corrections in terms of the $z=0$ linear power spectrum $P_{0}(k)=P^{\rm L}(k,z=0)$ were given in Section 2.1. For our implementation, we want to take $P^{\rm L}(k,z)$ and $f(k,z)$ at the intended MG-Copter output redshift from CAMB/MGCAMB and use it as input to MG-Copter. Therefore we need to rewrite the expressions for the 1-loop power spectra in terms of $P^{\rm L}(k,z)$ instead of $P_{0}(k)$, using $F_{1}(k)=G_{1}(k)/f(k,z)$ and $P_{0}(k)=P^{\rm L}(k,z)/F_{1}^{2}(k)=f^{2}(k,z)P^{\rm L}(k,z)/G_{1}^{2}(k,z)$. The 22 correction terms are

[TABLE]

while the 13 correction terms are

[TABLE]

Note that in these expressions the only terms to contain massive neutrinos are $P^{\rm L}$ and $f$; all of the kernels $F_{i}$ and $G_{i}$ are unmodified. The $A$ and $B$ terms written in Eqs. 2.19 and 2.20 are also computed as convolutions of two linear power spectra with kernels, and thus are rewritten using the same method as for $P^{13}$ and $P^{22}$. We have implemented these equations in MG-Copter.

3 Validation

In order to validate our implementation of massive neutrinos in the MG-Copter code, we have tested its output against results from the fast, approximate N-body code MG-PICOLA, which is a modified version of L-PICOLA [59] that includes modified gravity [60] and massive neutrinos [61] and has been tested against full N-body simulations. In the legends of the figures that follow we shall refer to our modified MG-Copter code simply as Copter, and the MG-PICOLA code as COLA.

Throughout, we use paired-fixed MG-PICOLA simulations where we produce two simulations with fixed amplitudes, meaning the initial amplitudes of the Fourier modes of the density field are set to that of the ensemble average power spectrum, and paired, where the initial modes in the second simulation are mirrored compared to those of the first [62]. This procedure significantly reduces variance that arises from the sparse sampling of wavemodes without the need for averaging over a large number of density field realisations, and has been shown not to introduce a bias to the recovery of the mean properties of the Gaussian ensemble, despite the fixing introducing non-Gaussianity [63]. However, we also ran five additional MG-PICOLA simulations for each model with randomised realisations of the initial density field. The standard deviation in the power spectra of these additional five simulations is used for the error bars in the figures below unless explicitly stated otherwise. The modified gravity model considered here is the Hu-Sawicki $f(R)$ model, which has one free parameter $|f_{R0}|$ and we refer to $|f_{R0}|=10^{-5}$ and $|f_{R0}|=10^{-4}$ as F5 and F4 respectively. The velocity divergence field $\theta$ has been computed using the DTFE code [64]. The cosmological parameters used in this paper are the same as in [30]; $h=0.671$, $\Omega_{\rm m}=0.3175$, $\Omega_{\rm b}=0.049$, $A_{s}=2.215\times 10^{-9}$, and $n_{\rm s}=0.966$.

Note that a recent version update of MGCAMB improved the handling of massive neutrinos [65]. Although our results were produced using the previous version of MGCAMB, we have verified that for the parameters we use the difference in the linear power spectrum between the two versions is negligible.

We first study the comparison between MG-Copter and MG-PICOLA in the real-space power spectra, in Figs. 1, 2 and 3. Figure 1 shows the real-space non-linear power spectra at $z=1$ computed with both MG-PICOLA and MG-Copter. We display the density auto-correlation $P_{\delta\delta}$, the velocity divergence auto-correlation $P_{\theta\theta}$, and the density-velocity divergence cross-correlation $P_{\delta\theta}$, in the form $k^{3/2}P_{ij}$ for ease of viewing, for GR, F5, and F4 each with $0.0$eV, $0.06$eV, and $0.2$eV neutrinos. The error bars on the (paired-fixed) MG-PICOLA points are the standard deviation of the 5 additional (non-paired-fixed) MG-PICOLA simulations. In all cases, MG-Copter reproduces the results of the MG-PICOLA simulations very well up to the start of the quasi-non-linear scale around $k=0.1~{}h/{\rm Mpc}$ where perturbation theory begins to break down. The agreement between MG-Copter and MG-PICOLA persists to larger $k$ values for $P_{\theta\theta}$ and $P_{\delta\theta}$ than $P_{\delta\delta}$, which is consistent with the behaviour seen when MG-Copter was compared to full N-body simulations in Fig. 10 of [46].

Figure 2 displays the same data but presented as the ratio of the full non-linear power spectra to their linear components, which helps to show where the modelling of non-linearities with MG-Copter becomes inaccurate. Figure 3 again shows the same data but presented as the ratio of the power-spectra with and without massive neutrinos for both the $0.06$eV and $0.2$eV neutrinos. The scale up to which MG-Copter closely follows the results of the MG-PICOLA simulations is marginally improved due to taking the ratio between power spectra in two models.

Next, we look at the comparison between MG-PICOLA and MG-Copter with $\sigma_{v}$ fitted to the MG-PICOLA simulations in the non-linear redshift-space power spectra in Figs. 4, 5 and 6. Figure 4 shows the monopole $P_{0}$ and quadrupole $P_{2}$ of the redshift-space power spectra for GR, F5, and F4 gravity models each with $0.0$eV, $0.06$eV, and $0.2$eV neutrinos. We display the results computed from paired-fixed MG-PICOLA simulations and MG-Copter with the TNS velocity dispersion parameter $\sigma_{v}$ fitted to the MG-PICOLA simulations up to $k=0.15~{}h/{\rm Mpc}$ in the form $k^{3/2}P_{i}(k)$; the figure includes the best-fitting values of $\sigma_{v}$ (expressed in RSD displacement units ${\rm Mpc}/h$) and the reduced $\chi^{2}$ for each model. The error bars on the MG-PICOLA points are taken from the inverse covariance matrices used in the $\sigma_{v}$ fitting procedure, whose computation is described at the end of Section 2.1. The $\sigma_{v}$ fitting procedure prioritises recovering the monopole $P_{0}$, and thus the agreement between MG-Copter and MG-PICOLA is slightly worse for the quadrupole $P_{2}$. As expected, for each gravity model increasing the mass of the neutrinos leads to a decrease in the best-fitting value of $\sigma_{v}$ and the quality of the fit increases, while for a fixed neutrino mass increasing the strength of the modification of gravity from GR to F5 and then F4 leads to an increase in the best-fitting value of $\sigma_{v}$ and a slightly worse quality of fit. The reason for this behaviour is that enhancement to gravity leads to an increase in the velocities of galaxies around an overdensity, thus increasing the non-linear damping, while massive neutrinos have the opposite effect due to their suppression of structure formation. The quality of the fit is better when the non-linearity is smaller and vice versa. However, in all cases the quality of the fit of MG-Copter to MG-PICOLA is good up to quasi-non-linear scales.

Figure 5 displays the same data as Fig. 4 but presented as the ratio of the full non-linear multipoles to their linear counterparts computed with the Kaiser RSD model [42], while Fig. 6 presents the data of Fig. 4 as the ratio of the non-linear power-spectra multipoles with and without massive neutrinos for both the $0.06$eV and $0.2$eV neutrinos. The error bars on the MG-PICOLA points in these two figures represent the standard deviation of the 5 additional MG-PICOLA simulations. As in real-space, the scale up to which MG-Copter closely follows the results of the MG-PICOLA simulations is slightly improved due to taking the ratio between power spectra in two models.

We also quantify the ability of MG-Copter to recover the redshift-space multipole results of MG-PICOLA through $\chi^{2}_{m_{\nu}}$; the difference between the redshift-space multipoles with and without neutrino mass. In Fig. 7 we display $\chi^{2}_{m_{\nu}}$ as a function of the maximum comparison scale $k_{\rm max}$ for GR, F5, and F4 each with $0.06$eV and $0.2$eV neutrinos at $z=1$. Here, MG-Copter is fitted to the MG-PICOLA simulations up to $k_{\rm max}$ with the covariance computed assuming an ideal survey as described at the end of Section 2.1. The agreement in $\chi^{2}_{m_{\nu}}$ between MG-PICOLA and MG-Copter fitted to MG-PICOLA is excellent in all cases. This implies that MG-Copter with $\sigma_{v}$ fitted to simulations is capable of capturing the effect of massive neutrinos accurately.

4 Degeneracy

With the inclusion of modified gravity and massive neutrinos in MG-Copter the degeneracy between the two effects can be investigated.

4.1 Real- and redshift-space

We start by studying the degeneracy between modified gravity and massive neutrinos in real space.

In Fig. 8 we display the ratio of real-space power spectra in F4 gravity with $0.06$eV neutrinos in the left panel and $0.2$eV neutrinos in the right panel to a fiducial model which we take to be GR with $0.06$eV neutrinos at $z=1$. We show results for the density auto-correlation $P_{\delta\delta}$, the velocity divergence auto-correlation $P_{\theta\theta}$, and the density-velocity divergence cross-correlation $P_{\delta\theta}$. The results of paired-fixed MG-PICOLA simulations and of MG-Copter are plotted. The error bars on the MG-PICOLA results are computed using the standard deviation over the five additional simulations. In all cases, the results of MG-Copter agree well with those of MG-PICOLA up to quasi-non-linear scales around $k=0.1~{}h/{\rm Mpc}$. The left panel, where the neutrino masses are the same in both GR and F4, shows the scale-dependent enhancement of the real-space power spectra provided by F4 gravity. However, when heavier neutrinos are added to the F4 case, as in the right panel, this enhancement is opposed by the suppression effect of the massive neutrinos. Indeed, the right panel shows that $P_{\delta\delta}$ is a poor probe to distinguish between GR with $0.06$eV neutrinos and F4 with $0.2$eV neutrinos in this particular case. However, the two models remain distinguishable in $P_{\delta\theta}$ and $P_{\theta\theta}$, showing that velocity information has the potential to break the degeneracy between modified gravity and massive neutrinos. This was recently shown using the results of full N-body simulations [44]. However, neither $P_{\delta\theta}$ nor $P_{\theta\theta}$ can be measured directly by observations. Instead, it is necessary to extract the velocity information that is encoded within redshift-space distortions, and it is to this we turn our attention. We shall refer to GR with $0.06$eV neutrinos and F4 with $0.2$eV neutrinos as our two degenerate models.

In Fig. 9 we plot the redshift-space monopole and quadrupoles in F4 gravity with $0.2$eV neutrinos normalised to GR with $0.06$eV neutrinos computed with both MG-Copter and MG-PICOLA. For each model the MG-Copter result has been produced by fitting $\sigma_{v}$ to the paired-fixed MG-PICOLA simulation up to $k=0.15~{}h/{\rm Mpc}$ with the covariance computed assuming a DESI-like survey as detailed at the end of Section 2.1. The error bars on the MG-PICOLA results are computed using the standard deviation over five simulations with a boxsize of $1~{}{\rm Gpc}/h$ for each model. Firstly, this plot shows that modelling the redshift-space monopole and quadrupole using MG-Copter with $\sigma_{v}$ fitted to MG-PICOLA simulations works well. Secondly, for our degenerate models, while the monopole is still a poor probe for distinguishing between the models, the quadrupole, by virtue of the encoding of velocity information, displays differences between the two models and thus has the potential to break the degeneracy.

4.2 Redshift evolution

Our method also allows us to investigate how the degeneracy evolves with redshift in both real- and redshift-space.

In Fig. 10 we show the real-space power spectra in the ratio between the two degenerate models as in the right panel of Fig. 8 but at $z=0.5$ (left panel) and $z=1.5$ (right panel). In Fig. 11 we show the redshift-space power spectrum multipoles in the ratio between the two degenerate models as in Fig. 9 but at $z=0.5$ (left panel) and $z=1.5$ (right panel). These figures demonstrate that the degeneracy evolves significantly with redshift, both in real- and redshift-space. Figure 10 shows that while our two degenerate models had similar matter power spectra at $z=1$ it is easier to distinguish between the two models with the matter power spectrum at other redshifts.

In Fig. 12 we plot the difference between the redshift-space multipoles in the two degenerate models quantified through $\chi^{2}_{\mathrm{MG}+m_{\nu}}$ as a function of the maximum comparison scale $k_{\rm max}$. We show $\chi^{2}_{\mathrm{MG}+m_{\nu}}$ as computed by both MG-PICOLA and MG-Copter with $\sigma_{v}$ fitted to MG-PICOLA up to $k_{\rm max}$ with the covariance computed assuming a DESI-like survey as detailed at the end of Section 2.1. The results from both methods agree with each other very well. We plot $\chi^{2}_{\mathrm{MG}+m_{\nu}}$ at three redshifts $z=1.5,\ 1.0,\ 0.5$ and it is clear from these results, along with those in Figs. 10 and 11, that the ability to distinguish between the redshift-space multipoles of these two models evolves with redshift. This emphasises the potential for data at multiple redshifts to break the degeneracy. The tomographic nature of weak lensing observations make them well suited to this task, and the combination of redshift-space distortion measurements with weak lensing observations could prove one of the best probes for breaking the modified gravity-massive neutrino degeneracy. However, it should be noted that systematics associated with weak lensing such as baryonic effects and intrinsic alignments may impact the effectiveness of such a probe.

5 Conclusions

In this paper, we have studied the potential for redshift-space distortions to break the degeneracy between the enhancement of structure growth provided by modifications to gravity and suppression of structure growth due to massive neutrinos, at the level of the dark matter field. For combinations of modified gravity parameters and neutrino masses that have similar matter power spectra at a given redshift, the growth rates are different and will remain distinguishable. This degeneracy-breaking growth rate information is encoded via velocities into redshift-space distortions. To carry out this work, we have modelled the effects of both modified gravity and massive neutrinos on real- and redshift-space power spectra with Standard Perturbation Theory through the code MG-Copter. We find the implementation of modified gravity and massive neutrinos in MG-Copter produces a good agreement for both real- and redshift-space power spectra with the simulation results from the code MG-PICOLA in the case of Hu-Sawicki $f(R)$ gravity.

We have then investigated the degeneracy and shown that the quadrupole of the redshift-space power spectrum retains enough of the velocity information to distinguish between GR with light neutrinos and Hu-Sawicki $f(R)$ with heavy neutrinos. The logical next step is to confirm that we can use the computationally inexpensive modelling of RSD in MG-Copter to recover a fiducial combination of $|f_{R0}|$ and $m_{\nu}$ from a simulation. An important open question for this endeavour is whether the process of fitting $\sigma_{v}$ introduces a new degeneracy, where $\sigma_{v}$ can dampen the redshift-space multipoles of a model with incorrect $|f_{R0}|$ and $m_{\nu}$ values in a way that makes them difficult to distinguish from those of the fiducial simulation. Future work will focus on extending the modelling of RSD with modified gravity and massive neutrinos to dark matter halos and galaxies in order to bring this method closer to being able to use RSD observations to jointly constrain modified gravity and massive neutrinos. It will be important to study how our conclusions change when we consider biased tracers instead of the underlying dark matter, as the bias parameters may also introduce additional degeneracies.

We have also briefly studied how the degeneracy evolves with redshift. There is a clear evolution of the degeneracy with redshift even for the matter power spectrum; for combinations of modified gravity and neutrino mass parameters that give comparable matter power spectra at one redshift, the matter power spectra at another redshift are in general likely to be distinguishable. The tomographic nature of weak lensing is particularly well suited to investigating this approach to breaking the degeneracy, although weak lensing systematics such as baryonic effects and intrinsic alignments could cause complications. Alternatively, if modified gravity is only a low redshift effect, a constraint on neutrino mass from clustering at higher redshift, for example from HI intensity mapping [66], would help break the degeneracy.

A work that appeared shortly after this paper studied RSD in halos from simulations with $f(R)$ gravity and massive neutrinos and reached conclusions broadly similar to our own [67].

Acknowledgement

We would like to thank Benjamin Bose for assistance with MG-Copter. BSW is supported by the U.K. Science and Technology Facilities Council (STFC) research studentship. KK and HAW are supported by the European Research Council through 646702 (CosTesGrav). KK is also supported by the UK Science and Technologies Facilities Council grants ST/N000668/1. GBZ is supported by NSFC Grants 1171001024 and 11673025, and the National Key Basic Research and Development Program of China (No. 2018YFA0404503). Numerical computations for this research were done on the Sciama High Performance Compute (HPC) cluster which is supported by the ICG, SEPNet, and the University of Portsmouth.