COLA with massive neutrinos

TL;DR

This paper develops a computationally efficient method to incorporate massive neutrinos into COLA simulations, accurately modeling their scale-dependent effects on structure growth in cosmology.

Contribution

It introduces a way to compute scale-dependent growth factors for modified gravity and neutrinos in Lagrangian perturbation theory within COLA, matching full N-body results.

Findings

Achieves percent-level accuracy in matter power spectra up to k≈1 h/Mpc.

Successfully integrates massive neutrinos into COLA with minimal computational cost.

Demonstrates compatibility with existing scale-dependent growth models.

Abstract

The effect of massive neutrinos on the growth of cold dark matter perturbations acts as a scale-dependent Newton's constant and leads to scale-dependent growth factors just as we often find in models of gravity beyond General Relativity. We show how to compute growth factors for $\Lambda$CDM and general modified gravity cosmologies combined with massive neutrinos in Lagrangian perturbation theory for use in COLA and extensions thereof. We implement this together with the grid-based massive neutrino method of Brandbyge and Hannestad in $\texttt{MG-PICOLA}$ and compare COLA simulations to full $\it N$-body simulations of $\Lambda$CDM and $f(R)$ gravity with massive neutrinos. Our implementation is computationally cheap if the underlying cosmology already has scale-dependent growth factors and it is shown to be able to produce results that match $\it N$-body to percent level accuracy for both the total and CDM matter power-spectra up to $k\lesssim 1 h/$Mpc.

Tables

Table 1. Table 1 : The cosmological parameters for the simulations performed in this paper. Common to all simulations are Ω m = 0.3175 subscript Ω m 0.3175 \Omega_{\rm m}=0.3175 , Ω b = 0.049 subscript Ω 𝑏 0.049 \Omega_{b}=0.049 , n s = 0.966 subscript 𝑛 𝑠 0.966 n_{s}=0.966 , A s = 2.215 ⋅ 10 − 9 subscript 𝐴 𝑠 ⋅ 2.215 superscript 10 9 A_{s}=2.215\cdot 10^{-9} and h = 0.671 ℎ 0.671 h=0.671 .

$m_{\nu }$ (eV)	${\Omega }_{\mathrm{CDM}}$	${\Omega }_{\nu }$	${\sigma }_8$ ($\Lambda \mathrm{CDM}$)
0.0	0.2685	0.0	0.850
0.2	0.2637	0.0048	0.798
0.4	0.259	0.0095	0.752
0.6	0.2542	0.0143	0.712

Equations

$${{\Psi }_{\mathrm{s}}}_i=ϵ{{\Psi }_{\mathrm{s}}}_i^{(1)}+{ϵ}^2{{\Psi }_{\mathrm{s}}}_i^{(2)}+\dots ,$$

$${\delta }_{\mathrm{cb}}={\left|\frac{\partial (x,y,z)}{\partial (q_x,q_y,q_z)}\right|}^{-1}-1=ϵ{\delta }_{\mathrm{cb}}^{(1)}+{ϵ}^2{\delta }_{\mathrm{cb}}^{(2)}+\dots .$$

$${\delta }_{\mathrm{cb}}^{(1)}$$

$${\delta }_{\mathrm{cb}}^{(2)}$$

$$\frac{d^2}{d{\tau }^2}{\Psi }_{\mathrm{cb}i,i}^{(1)}=-{{\nabla }_x}^2{\Phi }_N.$$

$${{\nabla }_x}^2{\Phi }_N=\kappa {\delta }_{\mathrm{cb}}^{(1)},$$

$$\left(\frac{d^2}{d{\tau }^2}-\kappa \right){{\nabla }_q}^2{\varphi }_{\mathrm{cb}}^{(1)}(\vec{q},\tau )=0,$$

$$\frac{d^2{D}_{1,\mathrm{cb}}}{d{\tau }^2}-\kappa D_{1,\mathrm{cb}}=0.$$

$${\mathcal{F}}_q\left[{{\nabla }_x}^2{\Phi }_N(\vec{x},\tau )\right]=\kappa {\mu }_{\mathrm{MG}}(k,\tau ){\mathcal{F}}_q\left[{\delta }_{\mathrm{cb}}^{(1)}(\vec{x},\tau )\right],$$

$$\left(\frac{d^2}{d{\tau }^2}-\kappa {\mu }_{\mathrm{MG}}(k,\tau )\right){\varphi }_{\mathrm{cb}}^{(1)}(\vec{k},\tau )=0,$$

$$\frac{d^2{D}_{1,\mathrm{cb}}(k,\tau )}{d{\tau }^2}-\kappa {\mu }_{\mathrm{MG}}(k,\tau )D_{1,\mathrm{cb}}(k,\tau )=0,$$

$$D_{1,\mathrm{cb}}(k,\tau )={\left[1+{\left(\frac{D_{1,\mathrm{cb}}^{m_{\nu }=0}(\tau )}{1+y_{fs}(\chi ;f_{\nu })}\right)}^{0.7}\right]}^{p_{\mathrm{cb}}/0.7}{D}_{1,\mathrm{cb}}^{m_{\nu }=0}(\tau {)}^{1-p_{\mathrm{cb}}},$$

$$D_{1,\mathrm{cb}\nu }(k,\tau )={\left[f_{\mathrm{cb}}^{0.7/p_{\mathrm{cb}}}+{\left(\frac{D_{1,\mathrm{cb}}^{m_{\nu }=0}(\tau )}{1+y_{fs}(\chi ;f_{\nu })}\right)}^{0.7}\right]}^{p_{\mathrm{cb}}/0.7}{D}_{1,\mathrm{cb}}^{m_{\nu }=0}(\tau {)}^{1-p_{\mathrm{cb}}},$$

$$p_{\mathrm{cb}}\equiv \frac{1}{4}\left[5-\sqrt[]{1+24f_{\mathrm{cb}}}\right]\ge 0,$$

$$y_{fs}(\chi ;f_{\nu })=17.2f_{\nu }\left(1+0.488f_{\nu }^{-7/6}\right){\left(N_{\nu }\chi /f_{\nu }\right)}^2,$$

$$\chi =\frac{k}{\mathrm{Mp}{\mathrm{c}}^{-1}}{\Theta }_{2.7}^2({\Omega }_0{h}^2{)}^{-1}=\frac{k}{19.0}({\Omega }_0{H}_0^2{)}^{-1/2}(1+z_{\mathrm{eq}}{)}^{-1/2},$$

$$P_{\mathrm{cb}\nu }(k,z)={\left[\frac{D_{1,\mathrm{cb}\nu }(k,z)}{D_{1,\mathrm{cb}\nu }(k,z=0)}\right]}^2{P}_{\mathrm{cb}\nu }(k,z=0).$$

$$\frac{d^2{\Psi }_{\mathrm{cb}i,i}^{(2)}}{d{\tau }^2}-{\Psi }_{\mathrm{cb}j,i}^{(1)}\frac{d^2{\Psi }_{\mathrm{cb}i,j}^{(1)}}{d{\tau }^2}=-{{\nabla }_x}^2{\Phi }_N,$$

$${{\nabla }_x}^2{\Phi }_N=\kappa {\delta }_{\mathrm{cb}}^{(2)}.$$

$$\left(\frac{d^2}{d{\tau }^2}-\kappa \right){{\nabla }_q}^2{\varphi }_{\mathrm{cb}}^{(2)}=-\frac{\kappa }{2}\left[({{\nabla }_q}^2{\varphi }_{\mathrm{cb}}^{(1)}{)}^2-({{\nabla }_q}_i{{\nabla }_q}_j{\varphi }_{\mathrm{cb}}^{(1)}{)}^2\right].$$

$$\frac{d^2{D}_{2,\mathrm{cb}}}{d{\tau }^2}-\kappa D_{2,\mathrm{cb}}=-\kappa D_{1,\mathrm{cb}}^2.$$

$${{\nabla }_q}^2{\varphi }_{\mathrm{cb}}^{(2)}(\vec{q},{\tau }_{\mathrm{ini}})=\frac{1}{2}\left[{\left({{\nabla }_q}^2{\varphi }_{\mathrm{cb}}^{(1)}(\vec{q},{\tau }_{\mathrm{ini}})\right)}^2-{\left({{\nabla }_q}_i{{\nabla }_q}_j{\varphi }_{\mathrm{cb}}^{(1)}(\vec{q},{\tau }_{\mathrm{ini}})\right)}^2\right].$$

$$\begin{gathered} {\mathcal{F}}_x\left[{{\nabla }_x}^2\Phi (\vec{x},\tau )\right]=\kappa {\mu }^{(1)}(k_{\mathrm{E}},\tau ){\delta }_{\mathrm{cb}}^{(1)}({\vec{k}}_{\mathrm{E}},\tau )+\kappa {\mu }^{(2)}(k_{\mathrm{E}},\tau ){\delta }_{\mathrm{cb}}^{(2)}({\vec{k}}_{\mathrm{E}},\tau ) \\ +a^4{H}^2\int \frac{d^3{\vec{k}}_1{d}^3{\vec{k}}_2}{(2\pi {)}^3}{\delta }_D(\vec{k}-\overrightarrow{k_{12}}){\gamma }_2^{\mathrm{E}}(\vec{k},\overrightarrow{k_1},\overrightarrow{k_2},\tau ){\delta }_{\mathrm{cb}}^{(1)}(\overrightarrow{k_1},\tau ){\delta }_{\mathrm{cb}}^{(1)}(\overrightarrow{k_2},\tau ), \end{gathered}$$

$$\begin{gathered} {\mathcal{F}}_q\left[{{\nabla }_x}^2\Phi (\vec{x},\tau )\right]=\kappa {\mu }_{\mathrm{MG}}(k,\tau ){\tilde{\delta }}_{\mathrm{cb}}^{(2)}(\vec{k},\tau ) \\ +a^4{H}^2\int \frac{d^3{\vec{k}}_1{d}^3{\vec{k}}_2}{(2\pi {)}^3}{\delta }_D(\vec{k}-\overrightarrow{k_{12}}){\gamma }_2(\vec{k},\overrightarrow{k_1},\overrightarrow{k_2},\tau ){\delta }_{\mathrm{cb}}^{(1)}(\overrightarrow{k_1},\tau ){\delta }_{\mathrm{cb}}^{(1)}(\overrightarrow{k_2},\tau ), \end{gathered}$$

$${\varphi }_{\mathrm{cb}}^{(2)}(\vec{k},\tau )=-\frac{1}{2k^2}\int \frac{{\mathrm{d}}^3{\vec{k}}_1{\mathrm{d}}^3{\vec{k}}_2}{(2\pi {)}^3}{\delta }_D(\vec{k}-{\vec{k}}_{12}){\delta }_{\mathrm{cb}}^{(1)}(\overrightarrow{k_1},{\tau }_{\mathrm{ini}}){\delta }_{\mathrm{cb}}^{(1)}(\overrightarrow{k_2},{\tau }_{\mathrm{ini}})D_{2,\mathrm{cb}}(\vec{k},\overrightarrow{k_1},\overrightarrow{k_2},\tau ),$$

$$\left[\frac{d^2}{d{\tau }^2}-\kappa {\mu }_{\mathrm{MG}}(k,\tau )\right]D_{2,\mathrm{cb}}(\vec{k},\overrightarrow{k_1},\overrightarrow{k_2},\tau )=-\kappa {\mu }_{\mathrm{MG}}(k,\tau )D_{1,\mathrm{cb}}(k_1,\tau )D_{1,\mathrm{cb}}(k_2,\tau )$$

$$\times \{1-\left[\frac{2{\mu }_{\mathrm{MG}}(k_1,\tau )-{\mu }_{\mathrm{MG}}(k,\tau )}{{\mu }_{\mathrm{MG}}(k,\tau )}\right]\frac{({\vec{k}}_1\cdot {\vec{k}}_2{)}^2}{k_1^2{k}_2^2}$$

$$+\frac{2a^4{H}^2{\gamma }_2(\vec{k},\overrightarrow{k_1},\overrightarrow{k_2},\tau )}{\kappa {\mu }_{\mathrm{MG}}(k,\tau )}\},$$

$$D_{2,\mathrm{cb}}(\vec{k},{\vec{k}}_1,{\vec{k}}_2,{\tau }_{\mathrm{ini}})$$

$$\frac{d{D}_{2,\mathrm{cb}}(\vec{k},{\vec{k}}_1,{\vec{k}}_2,{\tau }_{\mathrm{ini}})}{d\tau }$$

$${\varphi }_{\mathrm{cb}}^{(2)}(\vec{k},\tau )=$$

$$\times \int \frac{{\mathrm{d}}^3{\vec{k}}_1{\mathrm{d}}^3{\vec{k}}_2}{(2\pi {)}^3}{\delta }_D(\vec{k}-{\vec{k}}_{12}){\delta }_{\mathrm{cb}}^{(1)}(\overrightarrow{k_1},{\tau }_{\mathrm{ini}}){\delta }_{\mathrm{cb}}^{(1)}(\overrightarrow{k_2},{\tau }_{\mathrm{ini}})\left(1-\frac{(\overrightarrow{k_1}\cdot \overrightarrow{k_2}{)}^2}{k_1^2{k}_2^2}\right),$$

Full text

1 Introduction

Observations of neutrino flavour oscillations [1, 2] demand that at least two of the neutrino states are massive [3]. While the absolute mass of each of the three mass eigenstates is unknown, there are strong constraints on the difference in mass between the states which implies that at least one of them has a mass greater than $\mathtt{\sim}$0.06 eV.

Massive neutrinos are known to affect the formation of structure in the Universe by suppressing structure formation at small physical scales [4, 5, 6, 7, 8, 9]. This is in contrast to the effect of modified gravity models, which tend to enhance structure formation. With many current observations providing strong constraints on deviations from the $\Lambda$CDM ($\Lambda$-Cold-Dark-Matter) cosmology [10, 11], the opposing effects of massive neutrinos and modified gravity offer a potential alternative to $\Lambda$CDM when combined [12, 13, 14, 15, 16, 17, 18].

Multiple methods of including the effects of massive neutrinos in N-body simulations have been presented and used to make predictions for cosmological observables [19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]. High-resolution N-body simulations are still the only game in town if one wants to produce simulations that are accurate over a wide range of scales; however, they have the downside of being very computationally expensive to perform. For the next generation of galaxy surveys such as Euclid [30], eBOSS [31], DESI [32], WFIRST [33], LSST [34], and SKA [35], there is a need to develop large ensembles of galaxy mocks that include the effects of massive neutrinos in order to model the observables and their covariances. For such purposes faster methods are needed (see e.g. [36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50] and [51] for a recent comparison of some of these methods).

In previous work [52], we presented a code MG-PICOLA111The code can be found at https://github.com/HAWinther/MG-PICOLA-PUBLIC while the original L-PICOLA code on which it is based can be found at https://github.com/CullanHowlett/l-picola to perform fast, approximate numerical simulations of structure formation in models that have scale-dependent growth using the COmoving Lagrangian Acceleration (COLA) approach [43] (see also [44] for a similar approach). The COLA approach differs from traditional N-body simulations in that we solve for the perturbations about the particle trajectories predicted by Lagrangian perturbation theory. This allows us to take large N-body time-steps, thus saving a lot of computations (COLA is generally 2-3 orders of magnitude faster than N-body), while at the same time keeping accuracy on the largest scales, making the COLA approach well suited for producing large ensembles of galaxy mocks. In this paper we demonstrate how to extend this method to incorporate both modified gravity and massive neutrinos simultaneously in a quasi-nonlinear manner. The method is quasi-nonlinear because we only consider a linear neutrino density perturbation, such that the total density matter perturbation is defined $\delta_{\rm m}=(\overline{\rho}_{\rm cb}\delta_{\rm cb}+\overline{\rho}_{\nu}\delta^{(1)}_{\nu})/\overline{\rho}_{\rm m}$ where we track non-linearities only in $\delta_{\rm cb}$. Here the subscripts ‘m’, ‘cb’, and ‘$\nu$’ represent the total matter, cold dark matter + baryon, and massive neutrino components respectively. $f_{\rm cb}$ and $f_{\nu}$ are the fractional contributions of those components to the total matter density parameter $\Omega_{\rm m}$. The superscript numbers refer to the order in the perturbative expansion. We discuss the validity of only considering the linear neutrino density perturbation in Appendix B. We note that other fast, approximate methods can also handle massive neutrinos, namely PINOCCHIO (PINpointing Orbit-Crossing Collapsed HIerarchical Objects) [49, 50].

The paper is organised as follows. In Section 2 we recap how the first order growth factor is calculated in $\Lambda$CDM without massive neutrinos (Section 2.1), modified gravity without massive neutrinos (Section 2.2), and $\Lambda$CDM with massive neutrinos (Section 2.3), before explaining our extension for cosmologies with both modified gravity and massive gravity in Section 2.4 and presenting the results. In Section 3 we repeat this for the second order growth factors. In Section 4 we describe how these growth factors are implemented into the COLA approach. We present and discuss the outputs of our extended code in Section 5 and then conclude in Section 6.

2 First order growth factors

2.1 $\Lambda$CDM without massive neutrinos

In Lagrangian perturbation theory (LPT; see [53] for a review) the position of a particle of species ’s’ ${x_{\mathrm{s}}}_{i}$ is written in terms of its initial position ${q_{\mathrm{s}}}_{i}$ and a displacement field ${\Psi_{\mathrm{s}}}_{i}$ as ${x_{\mathrm{s}}}_{i}={q_{\mathrm{s}}}_{i}+{\Psi_{\mathrm{s}}}_{i}(\vec{q},\tau)$. Here $\tau$ is defined such that $d\tau=dt/a^{2}$. We expand the displacement field in a perturbation series

[TABLE]

and by assuming ${\Psi_{\mathrm{s}}}_{i}$ is curl-free we can write ${\Psi_{\mathrm{s}}}_{i}^{(n)}={\nabla_{q}}_{i}\phi^{(n)}_{\rm s}$ where $\phi^{(n)}_{\rm s}$ is a scalar field. We expand the cold dark matter (CDM) + baryon density contrast in a perturbation series, and use the Jacobian of the transformation between ${x_{\mathrm{cb}}}_{i}$ and ${q_{\mathrm{cb}}}_{i}$ to write

[TABLE]

Thus the CDM+baryon density contrast can be written in terms of the displacement-field order by order as

[TABLE]

[TABLE]

where $\Psi_{\mathrm{s}\ i,j}=\partial{\Psi_{\mathrm{s}}}_{i}/\partial{q}_{j}$. To first order, the Lagrangian equation of motion for the CDM+baryon component is

[TABLE]

In $\Lambda$CDM the only matter in the late-time Universe sourcing the Newtonian potential $\Phi_{N}$ is CDM+baryons, such that we can write $\delta_{\rm{m}}=\delta_{\rm{cb}}$. To first order the Poisson equation is thus

[TABLE]

where $\operatorname{\kappa}=4\pi G\overline{\rho}_{\rm{cb}}a^{4}=\frac{3}{2}\Omega_{\rm{cb}}H_{0}^{2}a$. The equations above combine to yield

[TABLE]

and ${\nabla_{q}}^{2}\phi^{(1)}_{\rm{cb}}(\vec{q},\tau_{\rm ini})=-\delta^{(1)}_{\rm{cb}}(\vec{q},\tau_{\rm ini})$ follows from assessing Eq. (2.3) at $\tau_{\rm ini}$. We can separate out the time dependence with $\phi^{(1)}_{\rm{cb}}(\vec{q},\tau)=D_{1,\rm{cb}}(\tau)\phi^{(1)}_{\rm{cb}}(\vec{q},\tau_{\rm ini})$ where the growth factor $D_{1,\rm{cb}}$ only depends on time and satisfies the simple ODE

[TABLE]

The initial conditions are set such that $D_{1,\rm{cb}}(\tau_{\rm ini})=1$ and $\frac{dD_{1,\rm{cb}}(\tau_{\rm ini})}{d\tau}=\left.\left(\frac{1}{a}\frac{da}{d\tau}\right)\right|_{\tau=\tau_{\rm ini}}$ corresponding to the growing mode in a matter dominated universe (Einstein-de Sitter).

2.2 Scale-dependent modified gravity without massive neutrinos

For modified gravity theories where the growth factor is scale-dependent the situation becomes a little more complicated than in $\Lambda$CDM. We will here consider a general, first order parametrisation of the gravitational potential in Fourier space

[TABLE]

where $\vec{k}$ is the Fourier coordinate with respect to coordinate $\vec{q}$, and $\mu_{\mathrm{MG}}$ is the effective Newton’s constant for a modified gravity theory. The addition of the scale-dependent $\mu_{\mathrm{MG}}$ in the Poisson equation means that, unlike in $\Lambda$CDM, we can no longer separate time and space. However, we can instead separate time for each Fourier mode. Inserting Eq. (2.9) into the Fourier transform of Eq. (2.5) with respect to $q$ and using the Fourier transform of Eq. (2.3) with respect to $q$ yields

[TABLE]

which allows us to make the split $\phi^{(1)}_{\rm{cb}}(\vec{k},\tau)=D_{1,\rm{cb}}(k,\tau)\phi^{(1)}_{\rm{cb}}(\vec{k},\tau_{\rm ini})$ where the growth factor satisfies

[TABLE]

with the initial field given by $k^{2}\phi^{(1)}_{\rm{cb}}(\vec{k},\tau_{\rm ini})=\delta^{(1)}_{\rm{cb}}(\vec{k},\tau_{\rm ini})$. We assume that at early times the modified gravity effects are negligible such that the initial conditions are still those for an Einstein-de Sitter universe: $D_{1,\rm{cb}}(k,\tau_{\rm ini})=1$, $\frac{dD_{1,\rm{cb}}(k,\tau_{\rm ini})}{d\tau}=\left.\left(\frac{1}{a}\frac{da}{d\tau}\right)\right|_{\tau=\tau_{\rm ini}}$. This second order differential equation can be solved numerically at each set of $(k,\tau)$ values for a given model with specified $\mu_{\mathrm{MG}}(k,\tau)$. We list $\mu_{\mathrm{MG}}(k,\tau)$ formulae for different modified gravity models in Appendix A.

2.3 $\Lambda$CDM with massive neutrinos

In order to account for the additional effect of massive neutrinos on the first order growth factor $D_{1,\rm{cb}}^{m_{\nu}=0}$ for matter perturbations in $\Lambda$CDM without massive neutrinos, we follow the practice of [54] and use the following fitting formulae to calculate the first order growth factor of CDM+baryon perturbations $D_{1,\rm{cb}}$ and total matter perturbations $D_{1,\rm{cb}\nu}$ in cosmologies with massive neutrinos:

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

and $f_{i}=\Omega_{i}/\Omega_{\mathrm{m}}$ is the density ratio for species $i$, $N_{\nu}$ is the number of massive neutrino species, and $\Theta_{2.7}$ is a measure of the CMB temperature at $z=0$ using $T_{\mathrm{CMB}}=2.7\Theta_{2.7}\mathrm{K}$.

This method allows us to insert $\mathrm{\Lambda}$CDM $D_{1,\rm{cb}}^{m_{\nu}=0}$ values that we calculate by solving Eq. (2.8) into Eqs. (2.12, 2.13) to calculate $D_{1,\rm{cb}}^{m_{\nu}}$ and $D_{1,\rm{cb}\nu}^{m_{\nu}}$ which then include the effect of massive neutrinos at linear order.

An alternative method for computing $D_{1,\rm{cb}}$ values would be to use the output of a Boltzmann solver code, such as CAMB [55]. In our implementation, we allow for both methods to be used as per the user’s preference.

Figure 1 displays a comparison between the outputs of Eq. (2.12) and CAMB. Specifically, the figure plots the ratio of the first order CDM+baryon growth factor at $z=0$ and $z=1$, $D_{1,\rm{cb}}(k,z=0)/D_{1,\rm{cb}}(k,z=1)$. The comparison is done for three GR+massive neutrino cosmologies with $m_{\nu}=[0.2,0.4,0.6]$ eV. The ratio has also been normalised to the $\Lambda$CDM case without massive neutrinos. Figure 1 shows that the output of Eq. (2.12) matches ${\tt{CAMB}}$ to an accuracy of $<1\%$ for neutrino masses $m_{\nu}\lesssim 0.6$ eV up to $k=1.0\ h/\mathrm{Mpc}$.

2.4 Scale-dependent modified gravity with massive neutrinos

Eisenstein and Hu originally only considered $\mathrm{\Lambda}$CDM $D_{1,\rm{cb}}^{m_{\nu}=0}$ values as input to Eqs. (2.12, 2.13). In this work we extend their idea by using $D_{1,\rm{cb}}^{m_{\nu}=0}$ values for modified gravity cosmologies without massive neutrinos as input instead, these having been calculated numerically by solving Eq. (2.11). This adds the effect of massive neutrinos to the modified gravity models, enabling us to compute growth factor values for modified gravity + massive neutrino (MG+$m_{\nu}$) cosmologies at linear order.

To test the fitting formula method for growth factors in MG+$m_{\nu}$ cosmologies, we first verified that our extension to the fitting formula method gives accurate values of $D_{1,\rm{cb}}$ and $D_{1,\rm{cb}\nu}$ by comparison with the output from MGCAMB, which is an extension of CAMB for modified gravity models [56, 57]. These comparisons can be seen in Figures 2-4 for $D_{1,\rm{cb}}$ in the Hu-Sawicki $f(R)$, symmetron, and dilaton modified gravity models. As in Figure 1, the plots show the ratio of the first order CDM+baryon growth factor at $z=0$ and $z=1$, $D_{1,\rm{cb}}(k,z=0)/D_{1,\rm{cb}}(k,z=1)$, the ratios have been normalised to the $\Lambda$CDM case without massive neutrinos, and the comparison is made for three different neutrino masses $m_{\nu}=[0.2,0.4,0.6]$ eV. Figure 2 shows that the ability of the fitting formula to recover MGCAMB first order growth values for $f(R)+m_{\nu}$ cosmologies decreases as $m_{\nu}$ increases. However, for the F4 model of $f(R)$ gravity (where $|f_{R0}|=10^{-4}$) that we consider in Figure 2 the output of Eq. (2.12) matches that of MGCAMB to an accuracy of $<1\%$ up to $k=1.0\ h/\mathrm{Mpc}$ even for $m_{\nu}=0.6$ eV. Similarly, for the values of parameters we have considered here, Eq. (2.12) matches MGCAMB to an accuracy of $<1\%$ up to $k=1.0\ h/\mathrm{Mpc}$ for both the symmetron and dilaton cosmologies with neutrino masses $m_{\nu}\lesssim 0.6$ eV.

We also wanted to test whether the growth factors calculated using this method could be used to accurately ‘backscale’ the linear total matter (CDM+baryon+massive neutrino) power-spectra at $z=0$ so that they closely matched the linear total matter power-spectra output at earlier $z$ by MGCAMB directly. This was done using the relationship

[TABLE]

In Figures 5-8, we display the resulting backscaled total matter linear power-spectra for the GR, $f(R)$, symmetron, and dilaton gravity models with $m_{\nu}$=[0, 0.2, 0.4, 0.6]eV at $z=1.3$. As for the first order growth factors, we find that, for the values of the model parameters considered, using the fitting formula method to backscale the $z=0$ linear total matter power-spectrum to $z=1.3$ recovers the same result as is output by MGCAMB directly at $z=1.3$ to an accuracy of $<1\%$ up to $k=1.0\ h/\mathrm{Mpc}$ for each of the cosmologies with neutrino masses $m_{\nu}\lesssim 0.6$ eV.

3 Second order growth factors

3.1 $\Lambda$CDM without massive neutrinos

To second order, the Lagrangian equation of motion for the CDM+baryon component is

[TABLE]

and the second order Poisson equation is

[TABLE]

Using $\Psi^{(n)}_{i}={\nabla_{q}}_{i}\phi^{(n)}$ and Eq. (2.4), the above two equations combine to yield

[TABLE]

Again we can separate out the time dependency completely with $\phi^{(2)}_{\rm{cb}}(\vec{q},\tau)=D_{2,\rm{cb}}(\tau)\phi^{(2)}_{\rm{cb}}(\vec{q},\tau_{\rm ini})$ and write

[TABLE]

For an Einstein-de Sitter Universe the physically relevant solution has $D_{2,\rm{cb}}=-\frac{3}{7}D_{1,\rm{cb}}^{2}$ so the initial conditions are taken to be $D_{2,\rm{cb}}(\tau_{\rm ini})=-\frac{3}{7}$ and $\frac{dD_{2,\rm{cb}}(\tau_{\rm ini})}{d\tau}=-\frac{6}{7}\left.\left(\frac{1}{a}\frac{da}{d\tau}\right)\right|_{\tau=\tau_{\rm ini}}$. The initial field $\phi^{(2)}_{\rm{cb}}(\vec{q},\tau_{\rm ini})$ satisfies

[TABLE]

3.2 Scale-dependent modified gravity without massive neutrinos

In a theory where the growth factor is scale-dependent, there are several additional complications at second order. When taking the Fourier transform of Eq. (3.1) with respect to $q$ we will need an expression for $\mathcal{F}_{q}\left[{\nabla_{x}}^{2}{\Phi}_{N}(\vec{x},\tau)\right]$ which will contain frame-lagging terms above linear order [58]. For the general case where the scale-dependence at nth order is encapsulated in the effective Newton’s constant $\mu^{(n)}(k,\tau)$ the Fourier transform with respect to $x$ up to second order is

[TABLE]

where $\vec{k}_{\rm E}$ is the Fourier coordinate with respect to $\vec{x}$, $\delta_{\rm{cb}}(\vec{k}_{\rm E},\tau)=\mathcal{F}_{x}\left[\delta_{\rm{cb}}(\vec{x},\tau)\right](\vec{k}_{\rm E})$ and $\vec{k}_{12}=\vec{k}_{1}+\vec{k}_{2}$. We can show that the corresponding Fourier transform with respect to $q$ at second order for the specific case of scale-dependent modified gravity without massive neutrinos is [59, 58, 52]

[TABLE]

where we have defined $\gamma_{2}=\gamma_{2}^{\rm E}+\frac{3}{2}\Omega_{\rm m}(\tau)\left[\mu_{\rm MG}(k,\tau)-\mu_{\rm MG}(k_{1},\tau)\right]\frac{\vec{k}_{1}\cdot\vec{k}_{2}}{k_{2}^{2}}$ and $\tilde{\delta}^{(2)}_{\rm{cb}}=\mathcal{F}_{q}\left[\delta^{(2)}_{\rm{cb}}(\vec{q},\tau)\right](\vec{k})+\mathcal{O}(\epsilon^{3})$. There is also a frame-lagging effect on $\gamma_{2}^{\rm E}$, but this would be a third order term so it is not included here.

We write $\phi^{(2)}_{\rm{cb}}$ (in a way that will become clear later) as an expansion

[TABLE]

where $\delta^{(1)}_{\rm{cb}}(\vec{k},\tau_{\rm ini})$ corresponds to the initial density field. With $\phi_{\rm cb}^{(2)}$ written in this form, substituting Eq. (3.7) into the Fourier space version of Eq. (3.1) leads to the second order differential equation for $D_{2,\mathrm{cb}}$

[TABLE]

with initial conditions

[TABLE]

For most modified gravity models $\gamma_{2}$ is only a function of the wavenumber norms $k,k_{1},k_{2}$ and the dot-product $\vec{k_{1}}\cdot\vec{k_{2}}$. The $\delta_{D}$ function in the integral for $D_{2}$ enforces $\vec{k_{2}}=\vec{k}-\vec{k_{1}}$, so solving Eq. (3.2) becomes a three-dimensional problem such that we only need to solve it for the combinations of $k$, $k_{1}$ and $\cos\theta\equiv\frac{\vec{k_{1}}\cdot\vec{k_{2}}}{k_{1}k_{2}}$ that correspond to a valid triangle in Fourier space.

Evaluating the integral in Eq. (3.8) at each time-step, without being able to rely on fast Fourier transforms, would ruin the speed of the COLA approach. As in [52], we therefore settle on an approximation for this term. We replace Eq. (3.8) with

[TABLE]

and replace Eq. (3.2) with

[TABLE]

which has initial conditions $\hat{D}_{2,\rm{cb}}^{\rm ini}=-\frac{3}{7}$ and $\frac{d\hat{D}_{2,\rm{cb}}^{\rm ini}}{d\tau}=-\frac{6}{7}\left.\left(\frac{1}{a}\frac{da}{d\tau}\right)\right|_{\tau=\tau_{\rm ini}}$.

3.3 $\Lambda$CDM with massive neutrinos

When including massive neutrinos at second order instead of first order, we no longer have fitting functions available that can add the effect of massive neutrinos to a growth factor calculated for CDM+baryons in cosmologies without massive neutrinos. In the method that follows, we treat the massive neutrinos as entirely linear such that $\delta_{\nu}=\delta^{(1)}_{\nu}$. Thus the only non-linearity comes from the CDM+baryon component $\delta_{\rm{cb}}=\delta^{(1)}_{\rm{cb}}+\delta^{(2)}_{\rm{cb}}$. This idea has been implemented and tested in N-body simulations [19] and in Standard Perturbation Theory (SPT) [60, 61, 62]. Reference [63] raised the issue that the treatment of massive neutrinos as purely linear causes problems stemming from the violation of momentum conservation. We discuss the impact of this on our work in Appendix B.

Now that we are including massive neutrinos, the gravitational potential is sourced by both CDM+baryons and the neutrinos, and thus the Fourier space Poisson equation up to second order is

[TABLE]

We can expand the expression $\overline{\rho}_{\rm{m}}\left(\delta^{(1)}_{\rm{m}}+\delta^{(2)}_{\rm{m}}\right)$ as

[TABLE]

such that Eq. (3.14) becomes

[TABLE]

Therefore the Fourier transform with respect to $q$ instead of $x$ is

[TABLE]

Inserting this into the full equation of motion up to second order yields

[TABLE]

We can separate out the first and second order terms into two equations. The first order equation is

[TABLE]

We use the Fourier space versions of Eq. (2.3) to write

[TABLE]

The first order equation can be solved numerically at each $(k,\tau)$, where the initial conditions are given by inserting the EdS initial conditions into the fitting formula Eq. (2.12). The second order equation is

[TABLE]

We use the Fourier space versions of Eq. (2.4) to write

[TABLE]

If we define $\phi^{(2)}_{\rm cb}(\vec{k},\tau)$ as in Eq. (3.8), rewrite the final term in Eq. (3.20) using the first order solution, and write out the Fourier transforms explicitly, then the second order equation becomes

[TABLE]

As for the case of modified gravity without massive neutrinos, the integral in Eq. (3.8) would be very slow to evaluate at each time-step. To circumvent this problem, we make an approximation by redefining $\phi^{(2)}_{\rm cb}$ as

[TABLE]

and setting $\mu_{m_{\nu}}\rightarrow f_{\rm cb}$, which leads to the following equation for $\hat{D}_{2,\rm{cb}}(k,\tau)$:

[TABLE]

For $f_{\rm cb}=1$ we recover the equation for $\Lambda$CDM without massive neutrinos. For a matter dominated Universe and for scales smaller than the neutrino free-streaming scale we have $D_{2,\rm cb}\simeq-\frac{3f_{\rm cb}}{3f_{\rm cb}+4(1-p_{\rm cb})^{2}}D_{1,\rm cb}^{2}$ which can be used to set the initial conditions when solving it numerically. The $\mu_{m_{\nu}}\rightarrow f_{\rm cb}$ approximation was previously made in [60, 61] for a $\Lambda$CDM+massive neutrino cosmology in Standard Perturbation Theory. They argued that the small value of $f_{\nu}$ suppresses the non-linear corrections to the above approximation of treating the massive neutrinos as an exclusively linear density perturbation. We have tested the effect of using $\mu_{m_{\nu}}(k)$ instead of $f_{cb}$ in Eq. (3.25) for the second order growth-factor in our COLA implementation (to be presented in the upcoming section). This change was found to have an negligible effect ($\lesssim 0.1-0.5$ % ) on the total matter power-spectrum for wavenumbers $k\lesssim 1\ h/$Mpc.

In Figure 9 we show the difference between $\mu_{m_{\nu}}$, calculated using Eq. (3.3), and $f_{\rm cb}$ in the range $0.01\leq k\leq 1.0\ h/\mathrm{Mpc}$ to highlight the consequences of making this approximation. Figure 9 shows that the approximation is less important as $k\rightarrow 1.0\ h/\mathrm{Mpc}$, but also shows that the approximation becomes less accurate as $m_{\nu}$ increases. Specifically, at $k=0.01\ h/\mathrm{Mpc}$, the percentage difference between $\mu_{m_{\nu}}$ and $f_{\rm cb}$ is 1.5/3.0/4.5$\%$ for $m_{\nu}=0.2/0.4/0.6$ eV respectively.

3.4 Scale-dependent modified gravity with massive neutrinos

In the case with both modified gravity and massive neutrinos, the Fourier space Poisson equation up to second order, with the potential $\Phi$ being sourced by both CDM+baryons and massive neutrinos, is

[TABLE]

Repeating the working of Sections 3.2 and 3.3, leads to

[TABLE]

where again we define $\gamma_{2}=\gamma_{2}^{\rm E}+\frac{3}{2}\Omega_{\rm m}(\tau)\left[\mu_{\mathrm{MG}}(k,\tau)\mu_{m_{\nu}}(k,\tau)-\mu_{\mathrm{MG}}(k_{1},\tau)\mu_{m_{\nu}}(k_{1},\tau)\right]\frac{\vec{k}_{1}\cdot\vec{k}_{2}}{k_{2}^{2}}$. Inserting this expression into the equation of motion yields the first order equation

[TABLE]

which can be solved numerically at each $(k,\tau)$, where the initial conditions are given by inserting the EdS initial conditions into the fitting formula Eq. (2.12). The second order equation is

[TABLE]

To speed up the calculation, we again make an approximation by redefining $\phi^{(2)}_{\rm cb}$ as

[TABLE]

and setting $\mu_{m_{\nu}}\rightarrow f_{\rm cb}$, which leads to the following equation for $\hat{D}_{2,\mathrm{cb}}$:

[TABLE]

For modified gravity models which only have late-time effects we can use the same initial conditions as for $\Lambda$CDM discussed above. Figure 9 shows that there is a negligible difference between the F4 model of $f(R)$ gravity and GR in the comparison between the values of $\mu_{m_{\nu}}$ and $f_{\rm cb}$ in the range $0.01\leq k\leq 1.0\ h/\mathrm{Mpc}$.

4 COLA Implementation

The implementation of massive neutrinos in the particle mesh part of the COLA algorithm is the grid-based method suggested in [19]. This method has been demonstrated to produce a matter power-spectrum that is accurate to $<1\%$ for neutrino masses $\sum m_{\nu}\lesssim 0.6$ eV. When we create the initial conditions for the CDM (CDM+baryon) particles we use the same initial seed to create a realisation of massive neutrinos using

[TABLE]

where $T_{\nu}$ and $T_{\rm cb}$ are the transfer functions of massive neutrinos and CDM+baryons respectively.

The massive neutrinos are kept in Fourier space for the duration of the simulation and are added to the source of the Poisson equation

[TABLE]

where the neutrino density at a given time $\tau$ is computed as

[TABLE]

In Appendix B we show a comparison of this scheme to an alternative scheme of modeling the non-linear neutrino density.

As for the COLA specific part we use the (scale-dependent) growth factors discussed in the previous sections. We have compared computing the growth factors using the fitting functions to directly solving the growth-ODEs with $\mu_{\nu}(k,\tau)=f_{cb}+f_{\nu}\frac{T_{\nu}(k,\tau)}{T_{cb}(k,\tau)}$ computed using transfer functions from CAMB or its alternatives. The difference between these two approaches was found to be negligible.

As long as the cosmological model we simulate already has scale-dependent growth then the additional computational cost of adding in neutrinos this way is almost negligible, but it does require some extra memory as we need to store the initial neutrino density field.

5 Simulation Results

We ran 5 COLA N-body simulations in a box of $B=512$ Mpc$/h$ with $N=512^{3}$ particles. The simulations were performed using the MG-PICOLA code. A smaller box of $B=256$ Mpc$/h$ with the same number of particles was used to check the convergence of the results, and Figure 10 shows, through comparison to the full N-body simulations of [14], that the CDM matter power-spectrum in our simulations can be trusted to percent level up to $k\sim 0.5-0.7~{}h/$Mpc. As found in our previous paper, the relative enhancement of the power-spectrum (i.e. when considering ratios of power-spectra as shown in the figures below) is accurate to larger $k$ values. The cosmological parameters for the simulations can be found in Table 1 and these are the same parameters as used by [14] where they performed combined massive neutrino and modified gravity simulations using a modified version of the simulation code Gadget [64, 65]. We will use these simulations to compare our results below. Ideally we would like to have run our simulations using exactly the same initial seed as the N-body simulations, however this was not available at the time we wrote this paper and we leave such a detailed comparison to future work.

In Figure 11 we show the suppression of the power-spectrum in $\Lambda$CDM as a function of neutrino mass in our simulations compared to the N-body results of [14]. We show power-spectrum results for both CDM+baryons $P_{\rm cb}(k,z)=\left<|\delta_{\rm{cb}}(\vec{k},z)|^{2}\right>$ and the total matter content (CDM+baryons+massive neutrinos) $P_{\rm m}(k,z)=\left<|f_{\rm{cb}}\delta_{\rm{cb}}(\vec{k},z)+f_{\nu}\delta_{\nu}(\vec{k},z)|^{2}\right>$.

In Figure 12 we show the results from simulations where we have both massive neutrinos and modified gravity. For the particular $f(R)$ model we study here having a total neutrino mass of $m_{\nu}\sim 0.4$ eV is seen to lead to a power-spectrum very similar to that of a standard $\Lambda$CDM model with massless neutrinos. This illustrates the degeneracy of massive neutrinos (suppressing growth) and modified gravity (enhancing growth).

Our COLA implementation gives power-spectrum (both for CDM and for the total matter) results that agree to $\lesssim 1\%$ accuracy for $k\lesssim 1~{}h/$Mpc to full N-body simulations of [14] for both $\Lambda$CDM and $f(R)$.

In Figures 13 and 14 we show the results for the halo mass function computed using Rockstar [66]. We note that the results of [14] were computed using a different halo-finder (SUBFIND) so the results are not directly comparable; however, the enhancement of the halo mass-function generally shows a good agreement.

6 Conclusions

Massive neutrinos, just like many modified gravity theories, lead to scale-dependent growth of matter perturbations. The effect can be described by an effective gravitational constant $\mu_{m_{\nu}}(k,\tau)$. We have implemented massive neutrinos in the COLA N-body code MG-PICOLA which allows for scale-dependent growth factors. For the particle mesh part of the COLA algorithm we use the grid-based method of [19] where massive neutrinos are kept in Fourier space and evolved linearly according to the neutrino growth factors. This method has previously been shown to be a good approximation for $\sum m_{\nu}\lesssim 0.6$ eV. A comparison to full N-body simulations of massive neutrino cosmologies, both for $\Lambda$CDM and $f(R)$ gravity, shows that we can match N-body results to percent level accuracy in both the total and CDM matter power-spectra with this approach.

We have also shown that the Eisenstein-Hu fitting formulae for the growth factors in massive neutrino cosmologies are a good approximation for a wide range of modified gravity theories as long as we replace the $\Lambda$CDM growth factor by the modified gravity counterpart.

In order to be able to judge the accuracy of our scheme more directly we would need to do a comparison to full N-body simulations where we use exactly the same initial conditions in COLA. This was not available to us at the time this study was performed, but such a comparison is something we plan to do in the future. In future work we also intend to study the degeneracy of massive neutrinos and modified gravity in detail and find a way to distinguish this case from $\Lambda$CDM; for example through signatures in redshift-space distortions.

Acknowledgement

We would like to thank Marco Baldi for providing us with the data for the $f(R)$ + massive neutrino N-body simulations. 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.

Appendix A Modified gravity models

In this section we give a brief overview of the two modified gravity models we are using in this paper focusing on the equations that are needed for our COLA implementation. For a more thorough review of these models, and modified gravity in general, see [67, 68].

A.1 $f(R)$ gravity

For $f(R)$ gravity [69] the growth of linear perturbations is determined by

[TABLE]

where $m(a)$ depends on the model in question. For the Hu-Sawicki model [70], which is the $f(R)$ model we will consider in this paper, $m(a)$ is given by

[TABLE]

where

[TABLE]

The $\gamma_{2}^{E}$ factor222Note that $\gamma_{2}$ as used in Eq. (3.4) is defined via $\gamma_{2}(k,k_{1},k_{2},\tau)=\gamma_{2}^{\rm E}(k,k_{1},k_{2},\tau)+\frac{3}{2}\Omega_{\rm m}(\tau)[\mu_{\rm MG}(k,\tau)\mu_{\nu}(k,\tau)-\mu_{\rm MG}(k_{1},\tau)\mu_{\nu}(k_{1},\tau)]\frac{\vec{k}_{1}\cdot\vec{k}_{2}}{k_{2}^{2}}$. is likewise given by [59]

[TABLE]

where

[TABLE]

A.2 Symmetron

In the symmetron model [71]:

[TABLE]

[TABLE]

if $a>a_{\star}$ and [math] otherwise where

[TABLE]

[TABLE]

A.3 Dilaton

In the dilaton model [72]:

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

Appendix B Comparison to SPT and alternative schemes for modeling the non-linear neutrino density

In this appendix we show a comparison of our code with linear theory and standard perturbation theory (SPT). The SPT results were obtained by following the method of [61] using the Einstein-de Sitter approximation. We have also carried out a test using an alternative scheme to include the neutrino density in the Poisson equation Eq. (4.2). In this scheme we use

[TABLE]

where $\delta_{cb}$ is the non-linear CDM+baryon density contrast instead of $\delta_{\nu}=\delta_{\nu}^{\rm lin}$ in Eq. (4.2). This is what [63] calls the improved external source scheme.

Fig. (15) shows that our implementation gives a result for the total matter power-spectrum that lies between SPT and linear theory on quasi-linear scales. The differences we see with respect to SPT for the slightly larger wavenumbers is expected (see Fig. (10) in [63]) as SPT slightly underestimates the power on these scales. The alternative scheme we tested is seen to slightly overestimate the power on linear scales and generally performs a bit worse on linear scales and at low redshift than the scheme we are using in this paper, especially for larger values of the neutrino mass.

Additionally, [63] warns of inaccuracies at large scales when treating neutrinos in a purely linear way with $\delta_{\nu}=\delta_{\nu}^{\rm lin}$, where there is a large deviation from the full non-linear scheme at large scales as a consequence of violation of momentum conservation. Indeed such inaccuracies would appear from Eqs. (3.23) and (3.29), but our method does not suffer from them because the approximations we make to maintain the speed of our COLA approach in the final pair of equations in Sections 3.3 and 3.4 demand that $D_{2,\rm cb}(\vec{k},\vec{k}_{1},\vec{k}_{2},\tau)=\left(1-\frac{(\vec{k}_{1}\cdot\vec{k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}}\right)\hat{D}_{2,\rm cb}(k,\tau)\rightarrow 0$ as $k\rightarrow 0$.