Investigating the degeneracy between modified gravity and massive neutrinos with redshift-space distortions
Bill S. Wright Kazuya Koyama Hans A. Winther Gong-Bo Zhao

Oscillations in the Sun with SONG: Setting the scale for asteroseismic investigationsBased on observations made at the Hertzsprung SONG telescope operated at the Spanish Observatorio del Teide on the island of Tenerife by the Aarhus and Copenhagen Universities and by the Instituto de Astrofísica de Canarias.
    M. Fredslund Andersen1
          P. Pallé2,3
          J. Jessen-Hansen1
		  K. Wang2,4
		  F. Grundahl1
		  T. R. Bedding5,1
		  T. Roca Cortes2,3
		  J. Yu5,1
		  S. Mathur2,3
		  R. A. Gacia6,7
		  T. Arentoft1
		  C. Régulo2,3
		  R. Tronsgaard8
  		  H. Kjeldsen1
          J. Christensen-Dalsgaard1
          
    Received January 31, 2019; accepted February 25, 2019
===============================================================================================================================================================================================================================================================================================================================





§ INTRODUCTION


Modifications to Einstein's theory of General Relativity (GR) can be considered when searching for an explanation of the late-time acceleration <cit.>. The simplest class of models, known as scalar-tensor theories <cit.>, 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 <cit.>, DESI <cit.>, WFIRST <cit.>, LSST <cit.>, and SKA <cit.>.

In this paper we will consider f(R) gravity <cit.>, and in particular the commonly studied Hu-Sawicki variant of f(R) <cit.>. 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 <cit.>, the presence of which demand that at least two of the neutrino states are massive <cit.>. 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_ν≳ 0.05 eV <cit.>. 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 λ_ fs <cit.>. 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 λ_ 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 <cit.>, even in models beyond ΛCDM that affect structure growth in a scale-independent way <cit.>.

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 <cit.>. 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 <cit.> and halo mass function <cit.> 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 <cit.>, 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 <cit.> and weak lensing tomographic information at multiple redshifts <cit.>; as well as techniques that are superior at distinguishing models such as machine learning <cit.>.

A different observable that has degeneracy breaking potential is that of redshift-space distortions (RSD) <cit.>. 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 <cit.>, whereas there is a strong stretching along the line-of-sight at non-linear scales commonly known as the Fingers-of-God (FoG) effect <cit.>. 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 <cit.>. 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  <cit.> 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 , the modified version of  developed in <cit.>, which is itself based on the approach presented in <cit.>.

We validate this implementation against simulations using the COmoving Lagrangian Acceleration (COLA) method <cit.>, 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 <ref> we explain our implementation of modified gravity and massive neutrinos in the Standard Perturbation Theory (SPT) formalism and  code. In Section <ref>, we show the results of tests validating our implementation against simulation results. In Section <ref> we use our new implementation to investigate the degeneracy and then conclude in Section <ref>.



§ 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  which computes large-scale structure observables using perturbation theory. For the redshift-space quantities,  depends on the TNS model of redshift-space distortions which is named after the authors of <cit.> (Taruya, Nishimichi, and Saito).



 §.§  and the TNS model
 
 <cit.> 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 <cit.>.

Starting from the continuity and Euler equations, assuming fluid quantities to be irrotational such that velocity field v can be expressed in terms of the velocity divergence θ=( ∇·v)/aH, and transforming to Fourier space yields

    a ∂δ(k)/∂ a + θ(k)
           = -∫d^3k_1 d^3k_2/(2π)^3δ_D( k - k_1 - k_2 ) α(k_1, k_2) θ(k_1)δ(k_2) , 
    
    
        a ∂θ(k)/∂ a + ( 2 + aH^'/H)θ(k)    - ( k/aH)^2 Φ(k)
        
       = -1/2∫d^3k_1 d^3k_2/(2π)^3δ_D( k - k_1 - k_2 ) β(k_1, k_2) θ(k_1)θ(k_2) ,

where a=1/(1+z) is the scale factor, y^'=∂ y/∂ a and the kernels α and β are given by

    α(k_1, k_2)    = 1 + k_1 ·k_2/|k_1|^2 ,
        
    
    β(k_1, k_2)    = (k_1 ·k_2)|k_1+k_2|^2/|k_1|^2|k_2|^2 .

The Poisson equation completes the above modified continuity and Euler equations

    -( k/aH)^2 Φ(k) = 3Ω_ m(a)/2δ(k) ,

where Ω_ m(a) = 8π Gρ_ m/3H^2. We want the n^ th order solutions of <ref>
to be of the form

    δ_n(k, a) = ∫ d^3k_1 … d^3k_n δ_D(k-k_1 … n)F_n(k_1,…,k_n, a) Δ(k_1)…Δ(k_n) ,


    θ_n(k, a) = ∫ d^3k_1 … d^3k_n δ_D(k-k_1 … n)G_n(k_1,…,k_n, a) Δ(k_1)…Δ(k_n) ,

where k_1 … n=k_1+…+k_n. Inserting these forms of the solutions into <ref> yields a generalised system of equations for the n^ th order kernels <cit.>

    ℒ̂[ F_n(k_1,…,k_n); G_n(k_1,…,k_n) ]
        = ∑_j=1^n-1[    -α(k_1 … j, k_j+1 … n) G_j(k_1,…,k_j) F_n-j(k_j+1,…,k_n); -1/2β(k_1 … j, k_j+1 … n) G_j(k_1,…,k_j) G_n-j(k_j+1,…,k_n) ] ,

where

    ℒ̂ = [                 ad/da                     1;               3Ω_ m/2 ad/da + ( 2 + aH^'/H) ] .

 solves this system of equations to compute the kernels F_i and G_i. The power spectra up to 1-loop are given as

    P^ 1-loop_ij(k) = P^ L_ij(k) + P_ij^13(k) + P_ij^22(k) ,

where the 1-loop corrections are defined by

    ⟨ x_2(k) y_2(k^') ⟩   = (2π)^3 δ_D(k+k^') P_xy^22(k) ,
        
    ⟨ x_1(k) y_3(k^') + x_3(k) y_1(k^') ⟩   = (2π)^3 δ_D(k+k^') P_xy^13(k) ,

where x and y can be δ or θ. Working these through, the final expressions for the 1-loop corrections in terms of the z=0 linear power spectrum P_0(k)=P^ L(k, z=0) are, for the 22 correction,

    P_δδ^22(k) =    2 k^3/(2π)^2∫_0^∞r^2 dr ∫_-1^1 P_0(kr)P_0(k√(1+r^2-2rx))F_2^2(k, r, x) dx ,
        
    
    
        P_δθ^22(k) =    2 k^3/(2π)^2∫_0^∞r^2 dr ∫_-1^1 P_0(kr)P_0(k√(1+r^2-2rx))F_2(k, r, x) G_2(k, r, x) dx ,
        
    
    
        P_θθ^22(k) =    2 k^3/(2π)^2∫_0^∞r^2 dr ∫_-1^1 P_0(kr)P_0(k√(1+r^2-2rx))G_2^2(k, r, x) dx ,

while for the 13 correction we have

    P_δδ^13(k) =    6 k^3/(2π)^2F_1(k)P_0(k) ∫_0^∞ r^2 P_0(kr)F_3(k, r, x) dr ,
        
    
    
        P_δθ^13(k) =    3 k^3/(2π)^2F_1(k)P_0(k) ∫_0^∞ r^2 P_0(kr)G_3(k, r, x) dr
        
       + 3 k^3/(2π)^2G_1(k)P_0(k) ∫_0^∞ r^2 P_0(kr)F_3(k, r, x) dr ,
        
    
    
        P_θθ^13(k) =    6 k^3/(2π)^2G_1(k)P_0(k) ∫_0^∞ r^2 P_0(kr)G_3(k, r, x) dr .

With the SPT real-space power spectra computed up to 1-loop order,  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 μ=cos(θ) is given by Eq. (18) of <cit.>, which we reproduce here with subtle changes due to the different definition of θ:

    P^(s)(k, μ) = D_ FoG[ kμσ_v ] { P_δδ(k) - 2μ^2 P_δθ(k) + μ^4 P_θθ(k) + A(k, μ) + B(k, μ) } ,

where D_ FoG is the Fingers-of-God damping function which we will discuss later. It is generally a function of k, μ, and the velocity dispersion σ_v. The power spectra P_δδ(k), P_δθ(k), and P_θθ(k) correspond to the density auto-correlation, density-velocity divergence cross correlation, and the velocity divergence auto-correlation respectively. A(k, μ) and B(k, μ) are correction terms given by

    A(k, μ)    = - kμ∫d^3 p/(2π)^3p_z/p^2{ B_σ(p, k - p, -k)-B_σ(p, k, -k-p) }, 
    
        B(k, μ)    = (kμ)^2 ∫d^3 p/(2π)^3 F(p)F(k-p) ,

where B_σ is the cross bispectrum defined by

    ⟨θ(k_1) {δ(k_2) - k^2_2z/k^2_2θ(k_2) }{δ(k_3) - k^2_3z/k^2_3θ(k_3) }⟩
     = (2π)^3 δ_D(k_1+k_2+k_3) B_σ(k_1,k_2,k_3)
         ,

and F(p) is defined as

    F(p)=p_z/p^2{ P_δθ(p) - p_z^2/p^2P_θθ(p) } .

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

    D_ FoG[ kμσ_v ] = exp( -k^2μ^2 σ_v^2 ) .

The velocity dispersion σ_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

    -2lnℒ = ∑_n ∑_l,l^'=0, 2( P^s_l, (k_n) - P^s_l, (k_n) ) Cov^-1_l,l^'(k_n) ( P^s_l^', (k_n) - P^s_l^', (k_n) )

for the first two multipoles. Expressions for the covariance matrix between the different multipoles Cov_l,l^' are given in Appendix C of <cit.>. 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  presented in Section <ref> we assume an ideal survey with survey volume V_ s=10  Gpc^3/h^3 and galaxy number density n̅_ g=4× 10^-3 h^3/ Mpc. For the study of the degeneracy in Section <ref>, we want to model a slightly more realistic scenario, so we assume a DESI-like survey with V_ s and n̅_ g as given in Table <ref> and redshift bin width Δ z=0.2. These values are computed using the information for emission line galaxies (ELGs) in Table V of <cit.>.



Thus the TNS model can be used to compute P^(s)(k, μ) with the input of P_δδ, P_δθ, P_θθ at 1-loop order from .



 §.§ Adding modified gravity
 

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

The modifications to gravity can be included in the Poisson equation, which up to 3^ rd order becomes

    -( k/aH)^2 Φ(k) = 3Ω_ m(a)/2δ(k)μ(k, a) + S(k)  ,

where μ(k, a)=G_ eff(k, a)/G is an effective Newton's constant[Not to be confused with the line-of-sight angle parameter μ, which will always be presented without arguments.] and the non-linear source term S(k) up to 3^ rd order is

    S(k) =    ∫d^3k_1 d^3k_2/(2π)^3δ_D(k-k_12)γ_2(k, k_1, k_2, a) Δ(k_1)Δ(k_2)
        
       + ∫d^3k_1 d^3k_2 d^3k_3/(2π)^3δ_D(k-k_123)γ_3(k, k_1, k_2, k_3, a) Δ(k_1)Δ(k_2)Δ(k_3) .

While the effective Newton's constant μ(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 γ_2 and γ_3 terms provide the leading order description of this screening in perturbation theory.

Using the same form for the n^ th order solutions as in <ref>, the new system of equations for the n^ th order kernels is

    ℒ̂[ F_n(k_1,…,k_n); G_n(k_1,…,k_n) ]
        = ∑_j=1^n-1[                        -α(k_1 … j, k_j+1 … n) G_j(k_1,…,k_j) F_n-j(k_j+1,…,k_n); -1/2β(k_1 … j, k_j+1 … n) G_j(k_1,…,k_j) G_n-j(k_j+1,…,k_n) - N_n(k, k_1,…,k_n) ] ,

where

    ℒ̂ = [                 ad/da                     1;        3Ω_ m/2μ(k, a) ad/da + ( 2 + aH^'/H) ] ,

and

    N_2 =   γ_2(k, k_1, k_2)F_1(k_1)F_1(k_2) ,
        
    
    
        N_3 =    γ_2(k, k_1, k_23)F_1(k_1)F_2(k_2, k_3)
        + γ_2(k, k_12, k_3)F_2(k_1, k_2)F_1(k_3)
        
       + γ_3(k, k_1, k_2, k_3)F_1(k_1)F_1(k_2)F_1(k_3) .


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. (<ref>) are given as

    μ(k, a) =    1 + ( k/a)^2 1/3Π(k, a) ,
        
    
    γ_2(k, k_1, k_2, a) =    - 9/48( kH_0/aH)^2 ( H_0^2Ω_m0/a^3)^2
        (Ω_m0-4a^3(Ω_m0-1))^5/a^15|f_R0|^2(3Ω_m0-4)^4
       ×1/Π(k, a)Π(k_1, a)Π(k_2, a) ,
        
    
    γ_3(k, k_1, k_2, k_3, a) =   ( kH_0/aH)^2 ( H_0^2Ω_m0/a^3)^3
        1/36Π(k, a)Π(k_1, a)Π(k_2, a)Π(k_3, a)Π(k_23, a)
       ×[ -45/8Π(k_23, a)/a^21|f_R0|^3((Ω_m0-4a^3(Ω_m0-1))^7/(3Ω_m0-4)^6) .
        
       . + H_0^2 ( 9/4a^15|f_R0|^2(Ω_m0-4a^3(Ω_m0-1))^5/(3Ω_m0-4)^4)^2 ]  ,

where

    Π(k, a) = ( k/a)^2 + H_0^2 (Ω_m0-4a^3(Ω_m0-1))^3/2|f_R0|a^9(3Ω_m0-4)^2 .




 §.§ Adding massive neutrinos
 

We have added support for massive neutrinos to the code  developed in <cit.>. Note that massive neutrinos were also added to the original  code in <cit.> using a similar approach. In our implementation, we follow the method of <cit.> and include massive neutrinos at the level of the linear real-space power spectra P^ L, P_δθ, L=f(k)P^ L, and P_θθ, L=f^2(k)P^ L without modifying the higher order SPT kernels. This allows us to take P^ L(k) and f(k) from  <cit.> (or  <cit.> for MG+m_ν) as input to ; note that a small modification to / 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 <cit.>, but for simplicity in the results that follow we have modelled the massive neutrinos as a single massive eigenstate with mass m_ν and two massless eigenstates.

The free-streaming of massive neutrinos causes suppression of P^ 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 <cit.> for an example. A linear approximation gives the amplitude of suppression to be -8f_ν where f_ν=Ω_ν/Ω_ m is the fraction of total matter in massive neutrinos <cit.>. Scale-dependent suppression also affects f(k), although the amplitude of this effect is much smaller, as can be seen in Fig. 5 of <cit.>.

The expressions for the 1-loop power spectra corrections in terms of the z=0 linear power spectrum P_0(k)=P^ L(k, z=0) were given in Section <ref>. For our implementation, we want to take P^ L(k, z) and f(k, z) at the intended  output redshift from / and use it as input to . Therefore we need to rewrite the expressions for the 1-loop power spectra in terms of P^ L(k, z) instead of P_0(k), using F_1(k)=G_1(k)/f(k, z) and P_0(k)=P^ L(k, z)/F_1^2(k)=f^2(k, z)P^ L(k, z)/G_1^2(k,z). The 22 correction terms are



    P_δδ^22(k) =    2 k^3/(2π)^2∫_0^∞r^2 dr ∫_-1^1 P^ L(kr, z)P^ L(k√(1+r^2-2rx), z)
        
       ×F_2^2(k, r, x)/F_1^2(kr)F_1^2(k√(1+r^2-2rx))dx ,
        
    
    
        P_δθ^22(k) =    2 k^3/(2π)^2∫_0^∞r^2 dr ∫_-1^1 P^ L(kr, z)P^ L(k√(1+r^2-2rx), z)
        
       × f(kr, z) f(k√(1+r^2-2rx), z) G_2(k, r, x)/G_1(kr) G_1(k√(1+r^2-2rx))
       ×F_2(k, r, x)/F_1(kr)F_1(k√(1+r^2-2rx))dx ,
        
    
    
        P_θθ^22(k) =    2 k^3/(2π)^2∫_0^∞r^2 dr ∫_-1^1 P^ L(kr, z)P^ L(k√(1+r^2-2rx), z)
        
       × f^2(kr, z)f^2(k√(1+r^2-2rx), z) G_2^2(k, r, x)/G_1^2(kr) G_1^2(k√(1+r^2-2rx))dx ,



while the 13 correction terms are



    P_δδ^13(k) =    6 k^3/(2π)^2 P^ L(k, z) ∫_0^∞ r^2 P^ L(kr, z) F_3(k, r, x)/F_1(k)F_1^2(kr)dr ,
        
    
    
        P_δθ^13(k) =    3 k^3/(2π)^2F_1(k) P^ L(k, z) ∫_0^∞ r^2 P^ L(kr, z) f(k, z)f^2(kr, z) G_3(k, r, x)/G_1(k)G_1^2(kr)dr
        
       + 3 k^3/(2π)^2f(k, z)P^ L(k, z) ∫_0^∞ r^2 P^ L(kr, z) F_3(k, r, x)/F_1(k)F_1^2(kr)dr ,
        
    
    
        P_θθ^13 =    6 k^3/(2π)^2P^ L(k, z) ∫_0^∞ r^2 P^ L(kr, z) f^2(k, z)f^2(kr, z) G_3(k, r, x)/G_1(k)G_1^2(kr)dr .



Note that in these expressions the only terms to contain massive neutrinos are P^ L and f; all of the kernels F_i and G_i are unmodified. The A and B terms written in <ref> 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 .



§ VALIDATION
 

In order to validate our implementation of massive neutrinos in the  code, we have tested its output against results from the fast, approximate N-body code , which is a modified version of  <cit.> that includes modified gravity <cit.> and massive neutrinos <cit.> and has been tested against full N-body simulations. In the legends of the figures that follow we shall refer to our modified  code simply as , and the  code as .

Throughout, we use paired-fixed  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 <cit.>. 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 <cit.>. However, we also ran five additional  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  and we refer to |f_R0|=10^-5 and |f_R0|=10^-4 as F5 and F4 respectively. The velocity divergence field θ has been computed using the  code <cit.>. The cosmological parameters used in this paper are the same as in <cit.>; h=0.671, Ω_ m=0.3175, Ω_ b=0.049, A_s=2.215× 10^-9, and n_ s=0.966.

Note that a recent version update of  improved the handling of massive neutrinos <cit.>. Although our results were produced using the previous version of , 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  and  in the real-space power spectra, in <ref>. Figure <ref> shows the real-space non-linear power spectra at z=1 computed with both  and . We display the density auto-correlation P_δδ, the velocity divergence auto-correlation P_θθ, and the density-velocity divergence cross-correlation P_δθ, in the form k^3/2P_ij for ease of viewing, for GR, F5, and F4 each with 0.0eV, 0.06eV, and 0.2eV neutrinos. The error bars on the (paired-fixed)  points are the standard deviation of the 5 additional (non-paired-fixed)  simulations. In all cases,  reproduces the results of the  simulations very well up to the start of the quasi-non-linear scale around k=0.1 h/ Mpc where perturbation theory begins to break down. The agreement between  and  persists to larger k values for P_θθ and P_δθ than P_δδ, which is consistent with the behaviour seen when  was compared to full N-body simulations in Fig. 10 of <cit.>.



Figure <ref> 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  becomes inaccurate. Figure <ref> again shows the same data but presented as the ratio of the power-spectra with and without massive neutrinos for both the 0.06eV and 0.2eV neutrinos. The scale up to which  closely follows the results of the  simulations is marginally improved due to taking the ratio between power spectra in two models.





Next, we look at the comparison between  and  with σ_v fitted to the  simulations in the non-linear redshift-space power spectra in <ref>. Figure <ref> 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.0eV, 0.06eV, and 0.2eV neutrinos. We display the results computed from paired-fixed  simulations and  with the TNS velocity dispersion parameter σ_v fitted to the  simulations up to k=0.15 h/ Mpc in the form k^3/2P_i(k); the figure includes the best-fitting values of σ_v (expressed in RSD displacement units Mpc/h) and the reduced χ^2 for each model. The error bars on the  points are taken from the inverse covariance matrices used in the σ_v fitting procedure, whose computation is described at the end of Section <ref>. The σ_v fitting procedure prioritises recovering the monopole P_0, and thus the agreement between  and  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 σ_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 σ_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  to  is good up to quasi-non-linear scales.



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





We also quantify the ability of  to recover the redshift-space multipole results of  through χ^2_m_ν; the difference between the redshift-space multipoles with and without neutrino mass. In Fig. <ref> we display χ^2_m_ν as a function of the maximum comparison scale k_ max for GR, F5, and F4 each with 0.06eV and 0.2eV neutrinos at z=1. Here,  is fitted to the  simulations up to k_ max with the covariance computed assuming an ideal survey as described at the end of Section <ref>. The agreement in χ^2_m_ν between  and  fitted to  is excellent in all cases. This implies that  with σ_v fitted to simulations is capable of capturing the effect of massive neutrinos accurately.





§ DEGENERACY
 

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



 §.§ Real- and redshift-space


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



In Fig. <ref> we display the ratio of real-space power spectra in F4 gravity with 0.06eV neutrinos in the left panel and 0.2eV neutrinos in the right panel to a fiducial model which we take to be GR with 0.06eV neutrinos at z=1. We show results for the density auto-correlation P_δδ, the velocity divergence auto-correlation P_θθ, and the density-velocity divergence cross-correlation P_δθ. The results of paired-fixed  simulations and of  are plotted. The error bars on the  results are computed using the standard deviation over the five additional simulations. In all cases, the results of  agree well with those of  up to quasi-non-linear scales around k=0.1 h/ 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_δδ is a poor probe to distinguish between GR with 0.06eV neutrinos and F4 with 0.2eV neutrinos in this particular case. However, the two models remain distinguishable in P_δθ and P_θθ, 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 <cit.>. However, neither P_δθ nor P_θθ 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.06eV neutrinos and F4 with 0.2eV neutrinos as our two degenerate models.



In Fig. <ref> we plot the redshift-space monopole and quadrupoles in F4 gravity with 0.2eV neutrinos normalised to GR with 0.06eV neutrinos computed with both  and . For each model the  result has been produced by fitting σ_v to the paired-fixed  simulation up to k=0.15 h/ Mpc with the covariance computed assuming a DESI-like survey as detailed at the end of Section <ref>. The error bars on the  results are computed using the standard deviation over five simulations with a boxsize of 1  Gpc/h for each model. Firstly, this plot shows that modelling the redshift-space monopole and quadrupole using  with σ_v fitted to  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.



 §.§ Redshift evolution


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





In Fig. <ref> we show the real-space power spectra in the ratio between the two degenerate models as in the right panel of Fig. <ref> but at z=0.5 (left panel) and z=1.5 (right panel). In Fig. <ref> we show the redshift-space power spectrum multipoles in the ratio between the two degenerate models as in Fig. <ref> 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 <ref> 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. <ref> we plot the difference between the redshift-space multipoles in the two degenerate models quantified through χ^2_MG+m_ν as a function of the maximum comparison scale k_ max. We show χ^2_MG+m_ν as computed by both  and  with σ_v fitted to  up to k_ max with the covariance computed assuming a DESI-like survey as detailed at the end of Section <ref>. The results from both methods agree with each other very well. We plot χ^2_MG+m_ν at three redshifts z=1.5, 1.0, 0.5 and it is clear from these results, along with those in <ref>, 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.



§ 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 . We find the implementation of modified gravity and massive neutrinos in  produces a good agreement for both real- and redshift-space power spectra with the simulation results from the code  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  to recover a fiducial combination of |f_R0| and m_ν from a simulation. An important open question for this endeavour is whether the process of fitting σ_v introduces a new degeneracy, where σ_v can dampen the redshift-space multipoles of a model with incorrect |f_R0| and m_ν 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 <cit.>, 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 <cit.>.



§ ACKNOWLEDGEMENT


We would like to thank Benjamin Bose for assistance with . 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.


JHEP


