An efficient implementation of massive neutrinos in non-linear structure formation simulations

Ali-Haïmoud, Yacine; Bird, Simeon

doi:10.1093/mnras/sts286

Abstract

Massive neutrinos make up a fraction of the dark matter, but due to their large thermal velocities, cluster significantly less than cold dark matter (CDM) on small scales. An accurate theoretical modelling of their effect during the non-linear regime of structure formation is required in order to properly analyse current and upcoming high-precision large-scale structure data, and constrain the neutrino mass. Taking advantage of the fact that massive neutrinos remain linearly clustered up to late times, this paper treats the linear growth of neutrino overdensities in a non-linear CDM background. The evolution of the CDM component is obtained via N-body computations. The smooth neutrino component is evaluated from that background by solving the Boltzmann equation linearized with respect to the neutrino overdensity. CDM and neutrinos are simultaneously evolved in time, consistently accounting for their mutual gravitational influence. This method avoids the issue of shot noise inherent to particle-based neutrino simulations, and, in contrast with standard Fourier-space methods, properly accounts for the non-linear potential wells in which the neutrinos evolve. Inside the most massive late-time clusters, where the escape velocity is larger than the neutrino thermal velocity, neutrinos can clump non-linearly, causing the method to formally break down. It is shown that this does not affect the total matter power spectrum, which can be very accurately computed on all relevant scales up to the present time.

neutrinos, cosmology: dark matter, cosmology: large-scale structure of Universe, cosmology: theory

INTRODUCTION

The determination of neutrino masses lies on the boundary between particle physics, astrophysics and cosmology. Atmospheric and solar neutrino oscillations have allowed the measurement of two mass-squared differences Δm²₁₂ ≡ m₂² − m²₁ = 7.54^{+ 0.46}_{− 0.39} × 10^{− 5} and |Δm²₁₃| ≡ |m₃² − (m²₁ + m₂²)/2| = 2.43^{+ 0.12}_{− 0.16} × 10^{− 3} eV² (central values and 2σ error bars from Fogli et al. 2012). These measurements do not allow us to pin down the individual neutrino masses, nor their hierarchy, since only the absolute value of Δm²₁₃ is measured – note, however, that next-generation oscillation experiments such as NOVA¹ aim to distinguish the neutrino hierarchy through second-order mixing effects. Current β-decay experiments imply an upper limit on the individual neutrino mass m_ν ≲ 2 eV (Kraus et al. 2005), and future experiments like KATRIN (KArlsruhe TRItium Neutrino experiment) should lower this limit by an order of magnitude (Eitel 2005).

Cosmological probes are sensitive to the absolute mass of neutrino species (see, for example, Lesgourgues & Pastor 2006 and Wong 2011 for recent reviews), and so far provide the tightest limits on the total neutrino mass M_ν ≡ ∑m_ν ≡ m₁ + m₂ + m₃, at the level of a few tenths of electron volts (eV). Increasingly sensitive future experiments aim to reach levels of precision at which a total neutrino mass ∑m_ν ∼ 0.1 eV or lower would be detectable (we refer the reader to Abazajian et al. 2011 for a recent compilation of current constraints and forecasts). Ultimately, the neutrino hierarchy may even be within reach (Takada, Komatsu & Futamase 2006; Jimenez et al. 2010).

Massive neutrinos affect the growth of perturbations through two effects. First, they change the background evolution. For a total neutrino mass less than a few eV, neutrinos are still relativistic at the time of matter-radiation decoupling but non-relativistic today. At fixed total density of non-relativistic matter today, the epoch of matter-radiation equality is moved to later times as the neutrino mass is increased. If one holds the density of baryons and CDM fixed instead, the angular-diameter distance and time of matter – dark energy equality are changed for different neutrino masses. In each case, the cosmic microwave background (CMB) anisotropy power spectrum is affected (Lesgourgues & Pastor 2006; Wong 2011), which has allowed CMB observations to set what is perhaps the most robust constraint on the neutrino mass, ∑m_ν < 1.3 eV (Komatsu et al. 2011). Upcoming data from the Planck satellite is expected to lower this bound by a factor of 2 (Planck Collaboration 2005; Abazajian et al. 2011).

The second effect of massive neutrinos is to slow down the growth of structure on scales smaller than their free-streaming length. This effect has long been understood in the linear regime (Bond, Efstathiou & Silk 1980; Bond & Szalay 1983; Ma & Bertschinger 1995; Lesgourgues & Pastor 2006 and references therein), and is implemented in all modern Boltzmann codes (see for example Lewis & Challinor 2002; Lesgourgues & Tram 2011). For a total neutrino mass of a few tenths of eV, this effect becomes manifest at scales ≲100h⁻¹ Mpc and gets stronger in the non-linear regime at scales of a few tenths to a few tens of Mpc (Lesgourgues & Pastor 2006).

Combined with CMB anisotropy measurements, various probes of the matter distribution have already constrained the total neutrino mass to ∑m_ν ≲ 0.2–0.3 eV. These include galaxy surveys such as SDSS (Sloan Digital Sky Survey; de Putter et al. 2012; Xia et al. 2012), cosmic shear surveys such as CFHT-LS (Canada-France-Hawaii Telescope Legacy Survey; Ichiki, Takada & Takahashi 2009), Lyman α forest measurements (Seljak, Slosar & McDonald 2006; Viel, Haehnelt & Springel 2010) and galaxy cluster surveys (Vikhlinin et al. 2009). As the survey samples grow larger, the sensitivity to the neutrino mass is expected to reach ∑m_ν ∼ 0.1 eV (Abazajian et al. 2011) or even lower (see Takada et al. 2006 for a detailed Fisher-matrix forecast of the constraining power of future high-redshift surveys). In addition, future measurements of CMB lensing (Hall & Challinor 2012) and 21 cm surveys (Mao et al. 2008) have the potential to reach yet lower bounds. Massive neutrinos also affect the growth of haloes (Brandbyge et al. 2010; Marulli et al. 2011).

Each one of the aforementioned methods has its own advantages and limitations. Systematic errors may arise from using biased tracers of the underlying density field, from difficulties in modelling complex baryonic physics or from foreground contamination. In addition, this level of precision requires an accurate modelling of corrections from non-linear structure growth. Here we focus on the latter problem. Cosmological neutrinos only interact gravitationally and, because of their large thermal velocities, effectively constitute a hot dark matter (HDM) component. While it poses significant practical difficulties, the non-linear clustering of massive collisionless particles is now a relatively well-understood problem. The growth of pure cold dark matter (CDM) structure can be computed through collisionless N-body simulations, whose accuracy is limited primarily by available computational resources, which have been steadily increasing over time (see e.g. fig. 1 of Dehnen & Read 2011). By contrast, the effect of neutrinos on the non-linear growth of matter perturbations has only recently been investigated in depth. This delicate problem has been approached from two angles. On the one hand, extensions to perturbation theory allow for semi-analytic calculations of the CDM power spectrum in the weakly non-linear regime. Combined with a linear treatment of neutrinos, such methods can reach a few per cent accuracy up to k ∼ 1 h Mpc⁻¹ at high redshifts (z ≳ 2) (Saito, Takada & Taruya 2008, 2009; Wong 2008; Lesgourgues et al. 2009). Shoji & Komatsu (2009) go beyond linear theory for neutrinos and solve for both the CDM and neutrino overdensities to third order in perturbation theory, making however the simplifying assumption that neutrinos can be described by an ideal fluid with a constant Jeans scale. Hannestad, Haugbølle & Schultz (2012) incorporated neutrinos into N-body simulations using a similar fluid approach.

Figure 1.

Open in new tab Download slide

Non-linear scale (solid) and free-streaming scales for various neutrino masses, as a function of redshift.

On the other hand, the ‘exact’ solution to the problem of CDM+neutrino clustering can in principle be obtained from N-body simulations where both CDM and neutrinos are simulated as particles. This approach was recently taken by Brandbyge et al. (2008), Viel et al. (2010) and Bird, Viel & Haehnelt (2012) (see also references therein for earlier works). The main limiting factor of such an approach is numerical: since neutrinos do not significantly cluster below their free-streaming scale, their power spectrum on small scales is dominated by shot noise for any reasonable number of simulated particles,² increasingly so for lower neutrino masses and at earlier times. A vast amount of computational power and memory is therefore wasted to extract a small amount of information on the actual neutrino clustering. In order to bypass the shot noise issue, Brandbyge & Hannestad (2009) included neutrino perturbations in Fourier space, while still treating the CDM component as particles in N-body simulations. The density field of their neutrino component was only computed using linear theory, however, and did not account consistently for non-linear CDM perturbations sourcing the gravitational potential in which neutrinos evolve (Bird et al. 2012). Recently, Brandbyge & Hannestad (2010) suggested a hybrid method, where neutrinos are initially treated using the Fourier method with linear theory, and converted to particles as time progresses. One of the remaining sources of errors of this implementation is that the Fourier-space part does not account for the non-linear growth of gravitational potentials in the neutrino evolution.

The aim of this work is to close this gap and provide an efficient method for accurately computing the effect of massive neutrinos on the non-linear matter power spectrum, with a minimal modification to well-tested pure CDM N-body simulation methods. To do so, we follow the CDM with the tree-PN N-body code gadget-3 (Springel 2005; Viel et al. 2010), and compute analytically the linearized neutrino overdensity sourced by the full gravitational potential. The CDM and neutrinos are evolved simultaneously and self-consistently, accounting for their mutual gravitational influence. Our method is therefore semilinear, in the sense that we account exactly for the non-linear growth of CDM overdensities and gravitational potentials, but effectively use a linear transfer function to obtain the neutrino overdensity from the gravitational potential. For neutrinos of the mass allowed by current data, our implementation is sufficient on its own for a complete description of the matter power spectrum up to the present time, and, for z ≥ 1, the neutrino power spectrum as well. It only fails to resolve neutrino clustering in the most massive galaxy clusters, in which neutrinos do cluster non-linearly (Ringwald & Wong 2004); however, it should be easy to perform zoomed simulations of these clusters using our method as a base. An ancillary advantage of our implementation is that it can easily be added into an N-body code as a small ‘patch’ and does not require running an independent Boltzmann code in parallel. Furthermore, our code can easily include the exact contribution of neutrinos to the background expansion, and the effect of the neutrino hierarchy, which are problematic or expensive to include in particle implementations. It can be applied to regimes where shot noise severely limits the reliability of the particle method, such as low neutrino masses, the 21cm forest or large-volume galaxy mock catalogue simulations. Perhaps the most important property of our code is that it allows massive neutrinos to be included into simulations with negligible cost in CPU and memory, and with no loss of accuracy. Thus, it can be useful for interpreting the results of essentially any probe of large-scale structure, including CMB lensing, the Lyman α forest, weak lensing experiments or galaxy surveys (Cooray 1999; Kaplinghat, Knox & Song 2003; Wang et al. 2005; Takada et al. 2006; Gratton, Lewis & Efstathiou 2008; Ichiki et al. 2009; Schlegel, White & Eisenstein 2009; Vallinotto et al. 2009; Beaulieu et al. 2010; Carbone et al. 2011).

This paper is organized as follows. In Section 2, we define our notation and discuss characteristic scales and the regime of validity of our main approximation. We derive our main equations in Section 3. The practical interfacing with our N-body code is detailed in Section 4. We describe our results and compare them against other methods in Section 5. We discuss potential applications of our method in Section 6 and conclude in Section 7.

NOTATION, CHARACTERISTIC TIME-SCALES AND LENGTH-SCALES

Notation

Throughout this paper, we use τ for the conformal time [dτ = dt/a, where t is the physical time and a = (1 + z)⁻¹ is the scale factor] and work in units where c = G = k_B = 1. Length-scales and wavenumbers are comoving unless otherwise stated. Overdots denote differentiation with respect to conformal time and gradients are with respect to comoving lengths. The Hubble rate today is H₀ = 100 h km s⁻¹ Mpc⁻¹. When referring to an individual neutrino mass, we use the lower case m_ν. When referring to the sum of neutrino masses, we use the upper case M_ν ≡ ∑m_ν. The temperature of the (unperturbed) neutrino background today is

\begin{equation} T_{\nu ,0} = 1.95 \textrm { K} = 1.68 \times 10^{-4} \textrm { eV}. \end{equation}

(1)

Neutrinos being non-relativistic today, they contribute a fraction f_ν of the total non-relativistic matter abundance, given by

\begin{equation} f_{\nu } = \frac{\Omega _{\nu }}{\Omega _{\rm M}} = \frac{1}{\Omega _{\rm M} \,h^2} \frac{\sum m_{\nu }}{93.14 \,{\rm eV} } \approx 0.022 \,\frac{0.147}{\Omega _{\rm M} \,h^2}\, \frac{\sum m_{\nu }}{0.3 \,{\rm eV}}. \end{equation}

(2)

Relativistic to non-relativistic transition

We denote by |$\boldsymbol p$| the physical momentum of a neutrino and |$f_0(\boldsymbol p)$| the unperturbed neutrino phase-space density, giving the number of neutrinos and antineutrinos per unit physical phase-space in the absence of gravitational clustering. Neutrinos decouple from the baryon plasma at a temperature of about 1 MeV (Lesgourgues & Pastor 2006). Given that their masses are known to be less than 1 eV, this means that they were ultrarelativistic at decoupling. This implies that their unperturbed phase-space density is the redshifted relativistic Fermi–Dirac distribution (since |$\boldsymbol p \propto a^{-1}$| in the absence of gravitational potentials),

\begin{equation} f_0(\boldsymbol p) = \frac{2}{h^3} \frac{1}{\exp (p/T_{\nu }) + 1}, \end{equation}

(3)

with T_ν(z) = (1 + z)T_{ν, 0}. A neutrino of mass m_ν becomes non-relativistic when its momentum falls below its mass, p ≲ m_ν. This occurs at a redshift z_nr such that

\begin{equation} 1 + z_{\rm nr} \approx \frac{m_{\nu }}{T_{\nu ,0}} \frac{T_{\nu }}{p} \approx 595 \frac{m_{\nu }}{0.1 \textrm { eV}} \frac{T_{\nu }}{p}. \end{equation}

(4)

The Fermi–Dirac distribution (3) is such that 50 per cent of neutrinos have a momentum p < 2.84 T_ν and 90 per cent have a momentum p < 5.47 T_ν. As a consequence, 50 per cent of neutrinos are non-relativistic for |$z \lesssim 200\,(m_{\nu }/0.1 \textrm { eV})$| and 90 per cent have become non-relativistic by |$z \lesssim 100\,(m_{\nu }/0.1 \textrm { eV})$|⁠. For the lowest masses considered (m_ν ≈ 0.05 eV), relativistic corrections to neutrino clustering may be of the order of unity at high redshifts (we start our simulations at z = 49); however, neutrinos are basically unclustered on all scales when they are quasi-relativistic, and the exact value of their very small inhomogeneities is then irrelevant for CDM clustering. We therefore treat neutrinos in the non-relativistic limit, i.e. assume p ≪ m_ν in all our derivations.³

Free-streaming scale

Neutrinos can free-stream across a comoving length-scale λ if the time it takes them to cross λ is much less than the Hubble time H⁻¹, i.e. if

\begin{equation} \frac{a \lambda }{v_{\nu }} \ll H^{-1}, \end{equation}

(5)

where v_ν(z) ∼ T_ν(z)/m_ν is the characteristic neutrino velocity. This equation defines a characteristic comoving free-streaming scalek_fs(z) ∼ aH/v_ν(z). A more detailed analysis in Section 3 will show that it is convenient to define the free-streaming scale as

\begin{equation} k_{\rm fs}(a) \equiv \left(\frac{3}{2} \Omega _{\rm M}(a) \overline{v^{-2}}\right)^{1/2} a H, \end{equation}

(6)

where |$\overline{v^{-2}}$| is the mean inverse velocity squared,

\begin{equation} \overline{v^{-2}} = \frac{2 \ln (2)}{3 \zeta (3)} \left(\frac{m_{\nu }a}{T_{\nu , 0}} \right)^2 \approx \left( 810\, {\rm km\, s^{-1}} \,(1+z)\frac{0.1\, \rm eV}{m_{\nu }}\right)^{-2}, \end{equation}

(7)

and Ω_M(a) is the relative contribution of the non-relativistic components to the total energy density at scale factor a. On scales much larger than the free-streaming scale, neutrinos behave like CDM. We shall show later on that for k ≫ k_fs, the (linearized) neutrino overdensity δ_ν relates to the total (possibly non-linear) matter overdensity δ_M as

\begin{equation} \delta _{\nu }(k \gg k_{\rm fs}) \approx \left(\frac{k_{\rm fs}}{k}\right)^2 \delta _{\rm M}. \end{equation}

(8)

Numerically, we obtain

\begin{equation} k_{\rm fs}(z) \approx \frac{0.08}{\sqrt{1+z}} \sqrt{\frac{\Omega _{\rm M}}{0.3}} \frac{m_{\nu }}{0.1 \, \textrm {eV}} \,h\, \textrm {Mpc}^{-1}. \end{equation}

(9)

We show the non-linear scale k_nl (defined such that the variance of CDM overdensity per log-k interval reaches unity on that scale) and free-streaming scale as a function of neutrino mass and redshift in Fig. 1. We see that for ∑m_ν < 0.6 eV, we always have k_fs < k_nl, with an increasing difference at large redshifts. Moreover, the non-linear matter power spectrum typically grows as k³P_M(k) ∝ k^α, with α ≲ 2 (Seljak 2000) and therefore we expect the power per log interval in the neutrino component to be decreasing for k ≳ k_nl > k_fs as k³P_ν(k) ∝ k^β, with β ≲ −2.

We therefore expect the neutrino component to be linearly clustered on all scales for masses below the current upper limit: for k < k_nl neutrinos cluster at most like the CDM, which is itself linearly clumped. For k > k_nl > k_fs, the neutrino power per logarithmic interval is a decreasing function of k.

This argument motivates our use of linear theory for the neutrino component. However, this does not give the full physical picture of neutrino clustering: the slowest neutrinos may in fact significantly cluster in massive haloes. Before moving to the core of our calculation, we first discuss in which conditions it may break down.

Neutrino capture in massive haloes

In this section, we qualitatively discuss under which conditions neutrinos may significantly cluster in a potential well of characteristic physical scale r₀ and characteristic depth ϕ₀ < 0. Let us consider, for simplicity, a spherical top hat potential, ϕ(r < r₀) = ϕ₀ and ϕ(r > r₀) = 0. We assume that r₀ is constant in time, but allow the depth ϕ₀ to vary slowly, on a Hubble time-scale, as is the case for the most massive haloes currently forming. We consider the fate of a particle on a purely radial trajectory. If the particle enters the potential at time t_i, with initial velocity v(t⁻_i) = v_i, its velocity upon entry becomes

\begin{equation} v\big(t_{\rm i}^+\big) = \sqrt{v_{\rm i}^2 + 2 |\phi _0(t_{\rm i})|}. \end{equation}

(10)

Since we have assumed the potential to be flat inside the halo, the particle's velocity is conserved until it reaches the other end, at time t_f. By then the gravitational potential has grown a little deeper, and the particle will escape only if its velocity is larger than the new escape velocity, i.e. if

\begin{equation} v_{\rm i}^2 + 2 |\phi _0(t_{\rm i})| > 2 |\phi _0(t_{\rm f})|. \end{equation}

(11)

Provided the crossing time is short compared to the evolution time-scale of the potential, we can Taylor-expand this equation and obtain the escape condition

\begin{equation} v_{\rm i}^2 > 2 (t_{\rm f} - t_{\rm i}) |\dot{\phi }_0|. \end{equation}

(12)

Inside the halo, provided 2|ϕ₀| ≫ v²_i, the velocity is approximately |$\sqrt{2 |\phi _0|}$| and the crossing time is therefore

\begin{equation} t_{\rm f} - t_{\rm i} \approx \frac{2 r_0}{\sqrt{2 |\phi _0|}}. \end{equation}

(13)

If we define the time-scale for variation of ϕ as

\begin{equation} \Delta t_{\phi } \equiv \frac{|\phi _0|}{|\dot{\phi }_0|}, \end{equation}

(14)

the escape condition for neutrinos becomes

\begin{equation} r_0 < \frac{1}{2} (H \Delta t_{\phi }) \frac{v}{H} \frac{v}{\sqrt{2|\phi _0|}}, \end{equation}

(15)

where we have purposefully inserted the Hubble parameter to make the order-unity parameter HΔt_ϕ appear. We see that for deep potential wells varying on the Hubble time-scale, the condition for neutrinos to truly free-stream and not be captured is more stringent than simply requiring their characteristic scale to be smaller than the (physical) free-streaming scale r_fs ∼ v/H. Equation (15) can also be turned into an escape condition for the neutrino momentum. For z ≈ 0, this is

\begin{eqnarray} &&\!\!\!\!\frac{p}{T_{\nu }} \gtrsim \frac{m_{\nu }}{T_{\nu , 0}} \frac{1}{\sqrt{H_0 \Delta t_{\phi }}} \left( 2 H_0 \,r_0 \sqrt{|2 \phi _0|}\right)^{1/2}\nonumber \\ &&\!\!\!\!\approx \left(H_0 \Delta t_{\phi }\right)^{-1/2} \frac{m_{\nu }}{0.1\ {\rm eV}} \left(\frac{r_0}{0.5 \ h^{-1}\ {\rm Mpc}}\right)^{1/2} \left(\frac{\sqrt{|\phi _0|}}{3000 \ {\rm km\,s^{-1}}}\right)^{1/2}. \nonumber \\ \end{eqnarray}

(16)

We have normalized the length-scale and depth of the gravitational well to values typical for the most massive haloes. The outcome of this analysis is that in massive haloes varying on a Hubble time-scale, neutrinos with momentum p ≲ T_ν are typically captured, while those with momentum p ≳ T_ν can escape. We emphasize that the cutoff value at p ≈ T_ν is a pure numerical coincidence, arising from the characteristic sizes and depths of massive haloes, and for neutrino masses of the order of 0.1 eV. Given that only about 6 per cent of neutrinos have a momentum p < T_ν, the qualitative picture that emerges from this analysis is that a relatively small fraction of neutrinos are efficiently bound to massive haloes, thereafter strongly clustering, while a majority remain weakly clustered.

Let us emphasize that although the crossing time of a neutrino in a massive halo is significantly shorter than the Hubble time (for the fiducial values of equation 16, Hδt ≈ 1/40 at z = 0), the relevant comparison is not of Hδt with unity but rather with the small ratio v²/ϕ (also ≈1/40 for our adopted fiducial values): even a small relative effect due to the change of the potential on a Hubble time-scale may translate into an order unity effect on the pre-entry velocity and prevent the neutrino from escaping the halo.

The analysis presented here is of course very simplified, for we have used an idealized (and unphysical) profile for the gravitational potential. Moreover, in reality there is no sharp boundary in momentum space between linear behaviour and strong clustering. It is however a robust statement to say that the characteristic momentum of neutrinos efficiently captured in massive haloes is of the order of their temperature. If the critical momentum were much larger than T_ν, the bulk of neutrinos would be captured and the neutrino overdensity in massive clusters would reach values comparable to those for the CDM. If it were much smaller than T_ν, then linear theory would be extremely accurate for neutrinos even in the vicinity of massive clusters. The actual situation is intermediate between these two regimes, as shown in the work of Ringwald & Wong (2004), who used N-one-body simulations to track the evolution of neutrinos around massive clusters. They found that for m_ν = 0.15 eV, linear theory underestimates neutrino clustering by a factor of ∼2 near the centre of 10¹⁴ M_⊙ clusters, and by a factor of ∼4 for 10¹⁵ M_⊙ clusters; for M ≤ 10¹² M_⊙, they found an excellent agreement with linear theory for m_ν = 0.15 eV.

With this caveat in mind, we now proceed to the main thrust of this paper: linear perturbations of neutrinos in a general gravitational potential.

LINEAR EVOLUTION OF HDM PERTURBATIONS IN A GENERAL BACKGROUND

In this section, we consider the linear evolution of an ensemble of non-relativistic HDM particles of mass m in a general gravitational potential ϕ in an expanding universe. We derive an integral equation known as the Gilbert equation (Gilbert 1966), which has been used in various analytic works on neutrinos (Bond & Szalay 1983; Brandenberger, Kaiser & Turok 1987; Singh & Ma 2003; Abazajian et al. 2005). This equation can easily be generalized to the relativistic case (Weinberg 2008; Shoji & Komatsu 2010; de Vega & Sanchez 2012a,b), although it takes on a more complicated form.

Vlasov equation

Derivation for non-relativistic particles in an expanding universe

In principle, the Boltzmann equation in an expanding universe should be derived in a consistent relativistic fashion (see for example Ma & Bertschinger 1995). However, for non-relativistic particles and in the Newtonian (ϕ ≪ 1) and subhorizon (k ≫ aH) limits, which always hold on the scales of interest, we can derive it more simply as follows.

We choose an arbitrary coordinate origin |$\boldsymbol {x}_0$| and denote by |$\boldsymbol r = a(\tau )(\boldsymbol x - \boldsymbol {x}_0)$| the physical coordinate of a particle, where |$\boldsymbol x$| is its comoving coordinate. We also denote by |$\boldsymbol p$| the physical momentum of the particle,

\begin{equation} \boldsymbol p \equiv m \frac{{\rm d} \boldsymbol r}{{\rm d} t} = m H \boldsymbol r + m \frac{{\rm d} \boldsymbol x}{ {\rm d} \tau } \equiv m H \boldsymbol r + a^{-1} \boldsymbol q, \end{equation}

(17)

where the last equality defines the comoving momentum |$\boldsymbol q \equiv a m {\rm d} \, \boldsymbol x/{\rm d} \tau$|⁠. The equation of motion for a test particle is

\begin{equation} \frac{{\rm d} \boldsymbol p}{{\rm d} t} = - m {\nabla }_{\boldsymbol r} \Phi , \ \ \textrm {i.e.} \ \ \frac{{\rm d} \boldsymbol q}{{\rm d} \tau } = - m a {\nabla }_{\boldsymbol x} \Phi - m a \frac{{\rm d}^2 a}{{\rm d} t^2} \boldsymbol r, \end{equation}

(18)

where Φ is the total gravitational potential. We now decompose Φ into a piece due to the background density and a piece sourced by the density perturbations:

\begin{equation} \Phi = \Phi _0 + \phi , \end{equation}

(19)

where

\begin{eqnarray} \Phi _0(\boldsymbol r) &=& -\frac{1}{2a}\frac{{\rm d}^2 a}{{\rm d} t^2} \boldsymbol {r}^2 = -\frac{1}{2} \left(H^2 + \frac{{\rm d} H}{{\rm d}t}\right) \boldsymbol {r}^2\nonumber \\ &=& \frac{2\pi }{3}\left(\overline{\rho } + 3 \overline{p}\right) \boldsymbol {r}^2 \end{eqnarray}

(20)

and ϕ is determined from the Poisson equation

\begin{equation} \nabla _{\boldsymbol r}^2 \phi = 4 \pi \delta \rho . \end{equation}

(21)

Equation (18) for the comoving momentum then simplifies to

\begin{equation} \frac{{\rm d} \boldsymbol q}{{\rm d} \tau } = - m a {\nabla }_{\boldsymbol x} \phi . \end{equation}

(22)

We now define |$f(\boldsymbol x, \boldsymbol q, \tau )$| as the HDM phase-space density, which is the number of particles per unit of physical phase-space (i.e. |${\rm d}N_{\rm hdm} = f(\boldsymbol x, \boldsymbol q, \tau ) {\rm d}^3 r {\rm d}^3p$|⁠). In the absence of collisions, the conservation of phase-space along a particle's trajectory gives the collisionless Boltzmann equation (also known as the Vlasov equation):

\begin{equation} \left.\frac{{\rm d} f}{{\rm d} \tau }\right|_{\rm traj} = \frac{\mathrm{\partial} f}{\mathrm{\partial} \tau } + \frac{{\rm d} \boldsymbol x}{{\rm d} \tau } \cdot {\nabla }_{\boldsymbol x} f + \frac{{\rm d} \boldsymbol q}{{\rm d} \tau }\cdot {\nabla }_{\boldsymbol q} f = 0. \end{equation}

(23)

Finally, using the definition of |$\boldsymbol q$| and equation (22) for its derivative, we find (independently of the chosen origin |$\boldsymbol {x}_0$|⁠)

\begin{equation} \frac{\mathrm{\partial} f}{\mathrm{\partial} \tau } + \frac{\boldsymbol q}{m a} \cdot {\nabla }_{\boldsymbol x} f - m a {\nabla }_{\boldsymbol x} \phi \cdot {\nabla }_{\boldsymbol q} f = 0. \end{equation}

(24)

Linearization and Fourier transform

In a homogeneous universe, |${\nabla }_{\boldsymbol x} \phi = 0$| and |$\nabla _{\boldsymbol x}f = 0$|⁠. As a consequence, the solution to the Vlasov equation must satisfy ∂f₀/∂τ = 0, i.e. be a function of momentum only. Finally, isotropy guarantees that f₀ = f₀(q) is a function of |$q \equiv |\boldsymbol q|$| only.

We now linearize the Vlasov equation about the homogeneous solution: |$f(\boldsymbol x, \boldsymbol q, \tau ) = f_0(q) + \delta f(\boldsymbol x, \boldsymbol q, \tau )$|⁠, assuming that the perturbation δf ≪ f₀ is induced by the presence of the gravitational potential ϕ [our formal expansion parameter is (ma)²ϕ/q²]. The linearized equation becomes (Brandenberger et al. 1987; Ma & Bertschinger 1995)

\begin{equation} \frac{\mathrm{\partial} \delta f}{\mathrm{\partial} \tau } + \frac{\boldsymbol q}{m a} \cdot {\nabla }_{\boldsymbol x} \delta f - \frac{m a}{q} \frac{{\rm d} f_0}{{\rm d} q} \boldsymbol {q} \cdot {\nabla }_{\boldsymbol x} \phi = 0. \end{equation}

(25)

We can now Fourier-transform this equation and obtain

\begin{equation} \frac{\mathrm{\partial} \delta f}{\mathrm{\partial} \tau } + {\rm i} \frac{\boldsymbol {q} \cdot \boldsymbol {k}}{ma} \delta f = {\rm i} \frac{m a}{q} \frac{{\rm d} f_0}{{\rm d} q} (\boldsymbol {q} \cdot \boldsymbol k) \phi , \end{equation}

(26)

Integral solution

For given comoving wavenumber and momentum, equation (26) is a linear first-order ordinary differential equation, and has the explicit solution

\begin{eqnarray} \delta f(\boldsymbol k, \boldsymbol q, \tau ) &=& \textrm {e}^{- {\rm i} \frac{\boldsymbol q \cdot \boldsymbol k}{m}(s-s_{\rm i})} \delta f(\boldsymbol k, \boldsymbol q, \tau _{\rm i}) \nonumber \\ &&+\ {\rm i}\frac{m}{q} \frac{{\rm d} f_0}{{\rm d} q} (\boldsymbol {q} \cdot \boldsymbol k) \int _{\tau _{\rm i}}^{\tau } \textrm {e}^{- {\rm i} \frac{\boldsymbol {q} \cdot \boldsymbol {k}}{m}(s-s^{\prime })} a(\tau ^{\prime }) \phi (\boldsymbol k, \tau ^{\prime }) \, {\rm d} \tau ^{\prime }, \nonumber \\ \end{eqnarray}

(27)

where τ_i is some initial conformal time, and we have used the variable (sometimes referred to as the ‘superconformal time’)

\begin{equation} s(\tau ) \equiv \int _{\tau _{\rm i}}^{\tau } \frac{{\rm d} \tau ^{\prime }}{a(\tau ^{\prime })}. \end{equation}

(28)

We can now evaluate the average of δf over the direction of momentum |$\hat{q}$| (this will be needed shortly):

\begin{eqnarray} &&\!\!\!\!\overline{\delta f}(\boldsymbol k, q, \tau ) \equiv \frac{1}{4 \pi } \int {\rm d}^2 \hat{q} \delta f(\boldsymbol k, \boldsymbol q, \tau ) \end{eqnarray}

(29)

\begin{eqnarray} \hphantom{\overline{\delta f}(\boldsymbol k, q, \tau )}&=& \overline{\delta f}^I(\boldsymbol k, q, \tau _{\rm i}, \tau )\nonumber \\ && + \ m \frac{{\rm d} f_0}{{\rm d} q} k \int _{\tau _{\rm i}}^{\tau } j_1\left(k\frac{q}{m}(s - s^{\prime })\right) a(\tau ^{\prime }) \phi (\boldsymbol k, \tau ^{\prime }){\rm d} \tau ^{\prime }, \nonumber \\ \end{eqnarray}

(30)

where j₁ is the order-one spherical Bessel function of the first kind (see for example Weisstein 2012b), and we have defined

\begin{equation} \overline{\delta f}^I(\boldsymbol k, q, \tau _{\rm i}, \tau ) \equiv \frac{1}{4 \pi } \int {\rm d}^2 \hat{q} \ \textrm {e}^{- {\rm i} \frac{\boldsymbol q \cdot \boldsymbol k}{m}(s - s_{\rm i})} \delta f(\boldsymbol k, \boldsymbol q, \tau _{\rm i}). \end{equation}

(31)

The mass density of the considered particle can be evaluated by integrating the phase-space density over momenta:

\begin{equation} \rho _{\rm hdm}(\boldsymbol x, \tau ) = m a^{-3} \int {\rm d}^3 q f(\boldsymbol x, \boldsymbol q, \tau ). \end{equation}

(32)

We can now obtain the density perturbation in Fourier space,

\begin{eqnarray} \delta \rho _{\rm hdm}(\boldsymbol k, \tau ) &=& m a^{-3}\int {\rm d}^3 q \delta f(\boldsymbol k, \boldsymbol q, \tau )\nonumber \\ &=& m a^{-3}\int 4 \pi q^2 dq \ \overline{\delta f}(\boldsymbol k, q, \tau ). \end{eqnarray}

(33)

Dividing by the mean density |$\overline{\rho }_{\rm hdm}$| (obtained from equation 32 with f = f₀), using equation (30) and integrating the second term by parts, we obtain the following expression for the density contrast |$\delta _{\rm hdm} \equiv \frac{\delta \rho _{\rm hdm}}{\overline{\rho }_{\rm hdm}}$|⁠:

\begin{eqnarray} \!\!\delta _{\rm hdm}(\boldsymbol k, \tau ) &=& \delta ^I_{\rm hdm}(\boldsymbol k, \tau _{\rm i}, \tau )\nonumber \\ && - k^2 \int _{\tau _{\rm i}}^{\tau } (s - s^{\prime }) \ \mathcal {I}\left[\frac{k}{m}(s - s^{\prime })\right] a(\tau ^{\prime })\phi (\boldsymbol k, \tau ^{\prime }) {\rm d} \tau ^{\prime }. \end{eqnarray}

(34)

In the above equation, |$\delta _{\rm hdm}^I(\boldsymbol k, \tau _{\rm i}, \tau )$| is the value of the density contrast evolved from τ_i to τ in the absence of any gravitational potential,

\begin{equation} \delta ^I_{\rm hdm}(\boldsymbol k, \tau _{\rm i}, \tau ) \equiv \frac{\int {\rm d}q q^2 \overline{\delta f}^I(\boldsymbol k, q, \tau _{\rm i}, \tau ) }{\int {\rm d}q \ q^2 f_0(q)}. \end{equation}

(35)

If we expand the phase-space density at the initial time on the basis of Legendre polynomials (see for example Weisstein 2012a),

\begin{equation} \delta f(\boldsymbol k, \boldsymbol q, \tau _{\rm i}) = \sum _{l=0}^{\infty }i^l \delta f_l(\boldsymbol k, q, \tau _{\rm i}) P_l(\hat{k} \cdot \hat{q}), \end{equation}

(36)

we obtain the following expression for δ^I:

\begin{equation} \delta ^I_{\rm hdm}(\boldsymbol k, \tau _{\rm i}, \tau ) = \sum _{l=0}^{\infty }\frac{\int \mathrm{d}q\ q^2 \delta f_l(\boldsymbol k, q, \tau _{\rm i})j_l\left(\frac{k q}{m}(s - s_{\rm i})\right)}{\int {\rm d}q q^2 f_0(q)}. \end{equation}

(37)

The dimensionless function |$\mathcal {I}$| in equation (34) is the Fourier transform of the unperturbed distribution function in momentum space, normalized so that |$\mathcal {I}(0) = 1$| (Brandenberger et al. 1987; Bertschinger & Watts 1988),

\begin{equation} \mathcal {I}[X; f_0] \equiv \frac{\int \mathrm{d}q\ j_0(q X) q^2 f_0(q) }{\int \mathrm{d}q \ q^2 f_0(q)}. \end{equation}

(38)

For neutrinos described by a relativistic Fermi–Dirac distribution, we provide an accurate fitting formula for this function in Appendix C.

Finally, since the gravitational potential is sourced by the total matter density through the Poisson equation

\begin{equation} k^2 \phi = - 4 \pi a^2 \overline{\rho }_{\rm M} \delta _{\rm M} = - \frac{3}{2} H_0^2 \frac{\Omega _{\rm M}}{a} \delta _{\rm M}, \end{equation}

(39)

where Ω_M is the matter fraction at the present day, we can rewrite equation (34) as

\begin{eqnarray} \delta _{\rm hdm}(\boldsymbol k, \tau ) &=& \delta ^I_{\rm hdm}(\boldsymbol k, \tau _{\rm i}, \tau )\nonumber \\ && +\, \frac{3}{2} H_0^2 \Omega _{\rm M} \int _{\tau _{\rm i}}^{\tau } \mathcal {I}\left[\frac{k}{m}(s - s^{\prime })\right] (s - s^{\prime }) \delta _{\rm M}(\boldsymbol k, \tau ^{\prime }) \, \mathrm{d} \tau ^{\prime }.\nonumber \\ \end{eqnarray}

(40)

Note that δ_hdm at time τ depends, in principle, on its value at all prior times through δ_M(τ′ < τ).

Discussion: limiting regimes

The unperturbed phase-space density is in general characterized by a typical comoving momentum q₀, such that f₀(q) decreases rapidly for q ≫ q₀. For neutrinos, for example, q₀ = T_{ν, 0}. We also assume that this is the case for δf; indeed, provided this is true initially, it will remain true as long as linear theory is valid, as can be seen from equation (27).

Large scales

In the limit kq₀(s − s_i)/m ≪ 1, we may Taylor-expand the spherical Bessel functions and obtain, up to terms of |$\mathcal {O}(k^2)$|⁠,

\begin{eqnarray} \delta _{\rm hdm}(k \rightarrow 0, \tau ) &=& \delta _{\rm hdm}(\boldsymbol k, \tau _{\rm i}) - a_i \theta _{\rm hdm}(\boldsymbol k, \tau _{\rm i}) (s - s_{\rm i})\nonumber \\ &&- k^2 \int _{\tau _{\rm i}}^{\tau } a(\tau ^{\prime }) (s - s^{\prime }) \phi (\boldsymbol k, \tau ^{\prime }) {\rm d} \tau ^{\prime }, \end{eqnarray}

(41)

where |$\delta _{\rm hdm}(\boldsymbol k, \tau _{\rm i})$| is the initial overdensity:

\begin{equation} \delta _{\rm hdm}(\boldsymbol k, \tau _{\rm i}) \equiv \frac{\int {\rm d} q q^2 \delta f_0(\boldsymbol k, q, \tau _{\rm i})}{\int {\rm d} q q^2 f_0(q)}\,, \end{equation}

(42)

and |$\theta _{\rm hdm}(\boldsymbol k, \tau _{\rm i})$| is the initial bulk velocity divergence:

\begin{equation} \theta _{\rm hdm}(\boldsymbol k, \tau _{\rm i}) \equiv \frac{1}{3} k \frac{\int {\rm d} q q^2 (q/m a_i) \delta f_1(\boldsymbol k, q, \tau _{\rm i})}{\int {\rm d} q q^2 f_0(q)}\,. \end{equation}

(43)

Equation (41) is nothing but the explicit solution to the (Newtonian, subhorizon) linearized fluid equations in an expanding universe (see for example Ma & Bertschinger 1995):

\begin{eqnarray} \dot{\delta }_{\rm hdm} + \theta _{\rm hdm} &=& 0, \end{eqnarray}

(44)

\begin{eqnarray} \dot{\theta }_{\rm hdm} + \frac{\dot{a}}{a} \theta _{\rm hdm} &=& k^2 \phi . \end{eqnarray}

(45)

We therefore recover the well-known fact that on scales much larger than their free-streaming length, HDM particles behave as a pressureless ideal fluid.

In general, the HDM component should behave like a cold species on scales much larger than its free-streaming length, regardless of whether it is linear. It is therefore important, for our linear approximation to work across all scales, that the CDM is indeed linear on scales larger than the free-streaming length, otherwise we would not get the correct long-wavelength behaviour. This is ensured by the hierarchy of scales k_fs < k_nl.

Small scales

Let us define the function

\begin{equation} E(\tau , \boldsymbol x, \boldsymbol q) \equiv \frac{q^2}{2 m} + m a^2 \phi (\tau , \boldsymbol x). \end{equation}

(46)

The change of E along a particle's trajectory is

\begin{equation} \frac{{{\mathop{\rm d}\nolimits} E}}{{{\mathop{\rm d}\nolimits} \tau }}\left| {_{{\rm{traj}}} \frac{{\partial E}}{{\partial \tau }} + \frac{{{\mathop{\rm d}\nolimits} x}}{{{\mathop{\rm d}\nolimits} \tau }} \cdot \frac{{\partial E}}{{\partial x}} + \frac{{{\mathop{\rm d}\nolimits} q}}{{{\mathop{\rm d}\nolimits} \tau }} \cdot \frac{{\partial E}}{{\partial q}} = \frac{{\partial E}}{{\partial \tau }},} \right. \end{equation}

(47)

\begin{equation} \left| {\frac{{{\mathop{\rm d}\nolimits} E}}{{{\mathop{\rm d}\nolimits} \tau }}\left| {_{{\rm{traj}}} } \right.} \right| \lesssim aHE, \end{equation}

(48)

provided that the potential varies on a Hubble time-scale. This is an upper limit: if ma²|ϕ| ≪ q²/m then the rate of change of E can in fact be much smaller than aHE.

Let us now consider a perturbation with scale small enough that the crossing time is much shorter than the Hubble time (which is itself of the order of or shorter than the time-scale over which E changes). In that case, the quantity E is an adiabatic invariant during the crossing of the perturbation. We assume that any particle within the perturbation can be traced back to a large distance from the perturbation at some earlier time τ_i, where |$m a_i^2 \phi (\tau _{\rm i}, \boldsymbol x_i) \ll q_i^2/m$| and hence E_i = E ≈ q²_i/(2m) (the distance must however be small enough to be crossed in a time-scale short compared to the time for E to change significantly). If furthermore at that initial time the phase-space density is nearly unperturbed, f(τ_i) ≈ f₀(q_i), we deduce, using conservation of phase-space, that, for E > 0,

\begin{eqnarray} f(\tau , \boldsymbol x, \boldsymbol q) &\approx & f_0(\boldsymbol q_i) \approx f_0(\sqrt{2 m E}) \nonumber \\ &=& f_0\left(\sqrt{q^2 + 2 (m a)^2 \phi (\tau , \boldsymbol x)}\right), \end{eqnarray}

(49)

and for E < 0, |$f(\tau , \boldsymbol x, \boldsymbol q) = 0$|⁠. This argument is generic to small-scale perturbations and does not assume linearity.⁴

If we now assume that |ϕ| < q²/(ma)², we can linearize equation (49) in ϕ and obtain

\begin{equation} \delta f(\tau , \boldsymbol x, \boldsymbol q) \approx \frac{(m a)^2}{q} \frac{{\rm d} f_0}{{\rm d} q} \phi (\tau , \boldsymbol x). \end{equation}

(50)

We find the resulting overdensity, after integrating δf by parts over momenta,

\begin{equation} \delta _{\rm hdm} \approx - \overline{v^{-2}} \phi , \end{equation}

(51)

where |$\overline{v^{-2}}$| is the mean inverse velocity squared,

\begin{equation} \overline{v^{-2}} \equiv (m a)^2 \frac{\int {\rm d}q f_0(q)}{\int {\rm d}q q^2 f_0(q)}. \end{equation}

(52)

Note, in passing, that these expressions require that df₀/dq → 0 as q → 0, and would not be valid for a Bose–Einstein distribution, for example.

Therefore, in the regime of shallow potentials, clustering on small scales is proportional to the square of the ratio of the characteristic velocity dispersion of the potential to a characteristic velocity of the HDM particles. Using the Poisson equation (39), we may also rewrite this in Fourier space as

\begin{equation} \delta _{\rm hdm}(k) \approx \frac{3}{2} k^{-2} a^2 H^2 \overline{v^{-2}} \Omega _{\rm M}(a) \delta _{\rm M} \equiv \left( \frac{k}{k_{\rm fs}}\right)^{-2} \delta _{\rm M}(k), \end{equation}

(53)

where we have used the definition of the free-streaming scale, equation (6), and Ω_M(a) is again the relative contribution of the non-relativistic components to the total energy density at scale factor a. This simple dependence is what motivated the exact numerical prefactors used in the definition of k_fs. This result was previously derived by Ringwald & Wong (2004). A similar result was obtained for baryon clustering on scales much smaller than the Jeans scale by Gnedin & Hui (1998) (in that case the Jeans scale plays the role of free-streaming scale, even though physically the absence of clustering is due to the pressure response rather than free-streaming as is the case for neutrinos).

We could also have used our integral expression equation (40) to reach the same results: if k ≫ k_fs, then (i) kq/m(s − s_i) ≫ 1 and the initial conditions become irrelevant, δ^I → 0, and (ii) the function |$\mathcal {I}$| decreases rapidly, such that only the times τ′ ≈ τ matter in the integral. In that case one finds, after a change of variables,

\begin{equation} \delta _{\rm hdm}(k \gg k_{\rm fs}) \approx - (m a)^2 \phi \int _0^{\infty } X \mathcal {I}(X) \, {\rm d} X, \end{equation}

(54)

which can be shown to give precisely equation (51).

APPLICATION TO SIMULATING MASSIVE NEUTRINOS

Equation (40) allows us to compute the density field for HDM from a possibly non-linear dark matter potential. We shall now specialize to the case where the HDM is a neutrino component interacting gravitationally with CDM and baryons, whose non-linear evolution is computed using an N-body code.

Initial conditions

The initial unperturbed distribution function for massive neutrinos which decoupled while relativistic is the Fermi–Dirac distribution,

\begin{equation} f_0(q) = \frac{g}{h^3}\frac{1}{\textrm {e}^{q/T_{\nu ,0}} + 1}, \end{equation}

(55)

where g = 2, the degeneracy of the neutrino species. Only the first two moments of the perturbed distribution function at the initial time are relevant: on small scales the initial conditions are rapidly forgotten, and on large scales the higher order moments scale as v^l and decay rapidly with l as neutrinos are non-relativistic. In principle, one should extract the full momentum information from a Boltzmann code at the startup redshift, and obtain |$\delta f_0(\boldsymbol k, q, \tau _{\rm i})$| and |$\delta f_1(\boldsymbol k, q, \tau _{\rm i})$|⁠. However, on large scales, neutrinos essentially behave as a cold species, with local overdensities and bulk flows roughly independent of their individual momenta, and on small scales, the exact form of initial conditions is irrelevant. These considerations allow us to assume the following simple form for the initial distribution function:

\begin{eqnarray} &&\!\!\!\!f_{\nu }(\boldsymbol x, \boldsymbol q, \tau _{\rm i}) = f_0(|\boldsymbol q-m_{\nu } a_i\boldsymbol v_{\nu }(\boldsymbol x, \tau _{\rm i})|)[1 + \delta _{\nu }(\boldsymbol x, \tau _{\rm i})] \nonumber \\ &&\!\!\!\!\quad \approx f_0(q)\left[1 + \delta _{\nu }(\boldsymbol x, \tau _{\rm i}) - a_i m_{\nu } \boldsymbol v_{\nu }(\boldsymbol x, \tau _{\rm i}) \cdot \hat{q} \frac{{\rm d} \ln f_0}{{\rm d} q}\right], \end{eqnarray}

(56)

where we have assumed that the bulk velocity |$\boldsymbol {v}_{\nu }$| is smaller than the characteristic random velocity. Fourier-transforming this equation, and using |$\boldsymbol {v}_{\nu }(\boldsymbol k) = - \frac{{\rm i}}{k} \hat{k} \theta _{\nu }(\boldsymbol k)$|⁠, we obtain the multipole moments of our approximate initial phase-space density:

\begin{eqnarray} \delta f_0(\boldsymbol k, q, \tau _{\rm i}) &=& f_0(q) \delta _{\nu }(\boldsymbol k, \tau _{\rm i}), \end{eqnarray}

(57)

\begin{eqnarray} \delta f_1(\boldsymbol k, q, \tau _{\rm i}) &=& \frac{{\rm d} f_0}{{\rm d} q} m_{\nu } k^{-1} a_i\theta _{\nu }(\boldsymbol k, \tau _{\rm i}), \end{eqnarray}

(58)

\begin{eqnarray} \delta f_l(\boldsymbol k, q, \tau _{\rm i}) &=& 0, \ \ \ l \ge 2. \end{eqnarray}

(59)

The initial conditions propagate to the current redshift according to equation (37), which gives us, after integrating the second term by parts,

\begin{equation} \delta ^I(\boldsymbol k, \tau _{\rm i}, \tau ) = \mathcal {I}_{s_{\rm i}, s} \left[\delta _{\nu }(\boldsymbol k, \tau _{\rm i}) - a_i \theta _{\nu }(\boldsymbol k, \tau _{\rm i}) (s - s_{\rm i}) \right], \end{equation}

(60)

where |$\mathcal {I}_{s_1, s_2} \equiv \mathcal {I}([s_2 - s_1]k/m)$| and |$\mathcal {I}$| was defined in equation (38).

Phase evolution

Storing the full three-dimensional Fourier transform of the gravitational potential at each time-step in order to perform the integral in equation (40) would require prohibitive amounts of memory. We therefore make some simplifying assumptions so that we need to only store the power spectrum.

Previous works (Brandbyge & Hannestad 2009; Viel et al. 2010) used the fully linear neutrino overdensity (i.e. computed assuming even the CDM is linear) and assumed that the random phases and relative amplitudes of the neutrino density field did not evolve, but were given by the initial conditions. This is true on linear scales, where there is no mode mixing, but not on non-linear scales, where neutrinos, even if linear themselves, evolve in the gravitational potential sourced by non-linear CDM perturbations, for which mode mixing is important.

Here again, we can take advantage of the hierarchy of free-streaming and non-linear scales k_fs(z) < k_nl(z) for currently allowed neutrino masses. For linear scales k < k_nl(z), the CDM phases are approximately constant and equal to the phase at the current time τ for any τ^′ < τ. For scales smaller than the free-streaming scale, the integral in equation (40) is dominated by the most recent times, such that |τ^′ − τ| ≪ τ. It is reasonable to assume that phase evolution and mixing are not too important during a short period of time, and therefore the phases at all relevant times are approximately equal to the phase at the present time. Because of the hierarchy of scales k_fs(z) < k_nl(z), we can therefore assume, for all scales, that the phase and relative amplitude at τ′ < τ is equal to that at the current time τ, and therefore

\begin{equation} \delta _{\rm M}(\boldsymbol k, \tau ^{\prime }) \approx \left(\frac{P_{\rm M}(k, \tau ^{\prime })}{P_{\rm M}(k, \tau )}\right)^{1/2} \delta _{\rm M}(\boldsymbol k, \tau ), \end{equation}

(61)

where P_M is the matter power spectrum. We also assume that the initial condition piece satisfies a similar relation

\begin{equation} \delta _{\nu }(\boldsymbol k, \tau _{\rm i}) \approx \left(\frac{P_{\nu }(k, \tau _{\rm i})}{P_{\nu }(k, \tau )}\right)^{1/2} \delta _{\nu }(\boldsymbol k, \tau ), \end{equation}

(62)

and that the ratio |$\theta _{\nu }(\boldsymbol k, \tau _{\rm i})/\delta _{\nu }(\boldsymbol k, \tau _{\rm i}) \equiv [\theta _{\nu }/\delta _{\nu }]_i(k)$| is a function of k only. The former is justified because initial conditions are only relevant on large (and linear) scales, where equation (62) is indeed correct, and the latter because structure is fully linear at our starting redshift.

Substituting equations (61) and (62) into equation (40), with the initial condition piece given by equation (60), we obtain, first, that the neutrino component is in phase with the CDM, and secondly, that the neutrino power spectrum is given by

\begin{eqnarray} P_{\nu }^{1/2}(k, \tau ) &=& \mathcal {I}_{s_{\rm i}, s} P_{\nu }^{1/2}(k, \tau _{\rm i}) \left\lbrace 1 - (s - s_{\rm i}) a_i [\theta _{\nu }/\delta _{\nu }]_{i}(k)\right\rbrace \nonumber \\ &&+ \frac{3}{2} \Omega _{\rm M} H_0^2 \int _{\tau _{\rm i}}^{\tau } \mathcal {I}_{s^{\prime }, s} P^{1/2}_{\rm M}(k, \tau ^{\prime }) (s - s^{\prime }){\rm d} \tau ^{\prime }, \end{eqnarray}

(63)

where |$\mathcal {I}_{s_1, s_2} \equiv \mathcal {I}([s_2 -s_1]k/m)$|⁠. Because neutrino and CDM overdensities are in phase (following from our assumptions), we can then recover the full neutrino density field from P^1/2_ν(k, τ) and the CDM density field by

\begin{equation} \delta _{\nu }(\boldsymbol k, \tau ) = \left(\frac{P_{\nu }(k, \tau )}{P_{\rm cdm}(k, \tau )}\right)^{1/2} \delta _{\rm cdm}(\boldsymbol k, \tau ). \end{equation}

(64)

Implementation

We implement the effect of massive neutrinos in Fourier space as a modification to the tree-PM-SPH code gadget -3 (Springel 2005; Viel et al. 2010). In order to lower the computational cost, and following Viel et al. (2010), we do not compute the short-range tree force due to and experienced by neutrinos. The size of a PM grid cell for our fiducial simulation is 0.5 h⁻¹Mpc, corresponding to a wavenumber k ≈ 13 h Mpc⁻¹, a scale two orders of magnitude smaller than the neutrino free-streaming length. Neutrino overdensities are therefore completely negligible compared to CDM overdensities below the grid scale, and we are justified in neglecting the tree force due to them. However, neglecting the tree force experienced by neutrinos is not strictly justifiable, and effectively amounts to having a lower resolution for the neutrinos than for the CDM (Viel et al. 2010). As we shall discuss below, we checked that our results are converged with respect to grid size, proving that there is a negligible contribution from scales smaller than our grid resolution, where the tree force would be active.

We evolve CDM particles (and baryons) as well as neutrino overdensities simultaneously as follows.

We generate adiabatic initial conditions using power spectra from camb and setting equal random relative amplitude and phases for the CDM, baryon and neutrino initial overdensities.
Given the total matter overdensity field up to time τ − Δτ, we evolve the CDM and baryons with our tree-PM code up to the next time-step τ. We then evaluate the Fourier-transformed CDM+baryon density, as is required for computing the long-range particle forces, and compute the CDM+baryon power spectrum at time τ.
In order to evaluate the neutrino overdensity at time τ with equation (63), we require the total matter power spectrum up to (and including) the current time τ. Equation (63) is therefore strictly speaking an implicit equation for P_ν(k, τ), and should be solved iteratively. Since the integrand vanishes at τ′ = τ [because of the (s − s′) term], and neutrinos contribute a small fraction of the dark matter, with a good initial guess a single iteration is sufficient. We show this in detail in Appendix B. We therefore approximate P_ν(k, τ) ≈ P_ν(k, τ − Δτ) to evaluate P_ν(k, τ′) inside the integral of equation (63) – we do so by interpolating over our tables of stored power spectra. We checked that P_ν(k, τ) obtained from this first iteration is converged at the level of ∼10⁻⁴ by comparing with the output of a second iteration, which is therefore not required.
We then recover the full phase information for the neutrino overdensity from equation (64), and evaluate the total matter overdensity at the current time-step τ from
\begin{equation} \delta _{\rm M}(\boldsymbol k, \tau ) = (1-f_\nu ) \delta _{\rm cdm, b}(\boldsymbol k, \tau ) + f_\nu \delta _{\nu }(\boldsymbol k, \tau ). \end{equation}
(65)
We finally compute and store the total matter power spectrum at time τ, and reiterate steps (ii) to (iv) until the end of the simulation.

Simulation details

The parameters of our simulations are shown in Table 1. Our cosmological parameters were similar to those in Bird et al. (2012). We set h = 0.7, the total matter fraction Ω_M = 0.3, the scalar spectral index n_s = 1 and the amplitude of the primordial power spectrum, A_s = 2.43 × 10⁻⁹. Our initial conditions were generated, from transfer functions generated by camb, using our own version of N-GenICs modified to use second-order Lagrangian perturbation theory (Scoccimarro 1998) and account correctly for baryons in the initial transfer function.⁵ Where it was important, we used Ω_b = 0.05.

Table 1.

Summary of simulation parameters. m_ν is the mass of one neutrino species. Rows marked with an asterisk were run using both Fourier and particle neutrinos. Cosmological parameters were the same for all simulations and are given in the text. Where neutrino particles were included, the same particle number was used as for the dark matter particles. Simulations with m_ν = 0 included massless neutrinos. Note we held Ω_M constant, making Ω_cdm dependent on Ω_ν. S03NH and S03IH had the same total neutrino mass as S05, but had three species with masses following the normal and inverted hierarchies, respectively.

Name	m_ν (eV)	Box ( Mpc h⁻¹)	N^1/3_CDM	z_i	N^1/3_bar
S00	0	256	512	49
S00P	0	256	1024	49
S05*	0.05	256	512	49
S10*	0.1	256	512	49
S15*	0.15	256	512	49
S20*	0.2	256	512	49
S40*	0.4	256	512	49
F10BA*	0.1	60	512	49	512
S05IC	0.05	256	512	99
S10P*	0.1	256	1024	49
S03NH	0.03, NH	256	512	49
S03IH	0.03, IH	256	512	49

Name	m_ν (eV)	Box ( Mpc h⁻¹)	N^1/3_CDM	z_i	N^1/3_bar
S00	0	256	512	49
S00P	0	256	1024	49
S05*	0.05	256	512	49
S10*	0.1	256	512	49
S15*	0.15	256	512	49
S20*	0.2	256	512	49
S40*	0.4	256	512	49
F10BA*	0.1	60	512	49	512
S05IC	0.05	256	512	99
S10P*	0.1	256	1024	49
S03NH	0.03, NH	256	512	49
S03IH	0.03, IH	256	512	49

Open in new tab

Table 1.

Summary of simulation parameters. m_ν is the mass of one neutrino species. Rows marked with an asterisk were run using both Fourier and particle neutrinos. Cosmological parameters were the same for all simulations and are given in the text. Where neutrino particles were included, the same particle number was used as for the dark matter particles. Simulations with m_ν = 0 included massless neutrinos. Note we held Ω_M constant, making Ω_cdm dependent on Ω_ν. S03NH and S03IH had the same total neutrino mass as S05, but had three species with masses following the normal and inverted hierarchies, respectively.

Name	m_ν (eV)	Box ( Mpc h⁻¹)	N^1/3_CDM	z_i	N^1/3_bar
S00	0	256	512	49
S00P	0	256	1024	49
S05*	0.05	256	512	49
S10*	0.1	256	512	49
S15*	0.15	256	512	49
S20*	0.2	256	512	49
S40*	0.4	256	512	49
F10BA*	0.1	60	512	49	512
S05IC	0.05	256	512	99
S10P*	0.1	256	1024	49
S03NH	0.03, NH	256	512	49
S03IH	0.03, IH	256	512	49

Name	m_ν (eV)	Box ( Mpc h⁻¹)	N^1/3_CDM	z_i	N^1/3_bar
S00	0	256	512	49
S00P	0	256	1024	49
S05*	0.05	256	512	49
S10*	0.1	256	512	49
S15*	0.15	256	512	49
S20*	0.2	256	512	49
S40*	0.4	256	512	49
F10BA*	0.1	60	512	49	512
S05IC	0.05	256	512	99
S10P*	0.1	256	1024	49
S03NH	0.03, NH	256	512	49
S03IH	0.03, IH	256	512	49

Open in new tab

In a previous work (Bird et al. 2012), we found that a significant limitation to the inclusion of massive neutrinos was the restriction on initial redshift due to the slightly relativistic behaviour of the neutrinos at early times. Although our method for computing the growth of neutrino perturbations is non-relativistic, we can at least easily include the correct contribution of neutrinos to the background evolution by computing

\begin{equation} \overline{\rho }_\nu (a) = a^{-3}\int \epsilon (q, a) f_0(q) 4 \pi q^2 {\rm d}q\,, \end{equation}

(66)

where |$\epsilon (q, a) = \sqrt{m^2_\nu + a^{-2} q^2}$| is the total energy of a neutrino with comoving momentum q. Thus, |$\overline{\rho }_\nu (a)$| interpolates between the behaviours of matter and radiation. This is extremely easy to implement with our Fourier method, but somewhat harder in a purely particle simulation, where the energy density is by default computed from the particles’ rest mass and does not account for their kinetic energy. For consistency, we also included the effect of radiation density in the background evolution, and, for CDM only simulations, the effect of massless neutrinos. Our particle-based simulations still have |$\overline{\rho }_{\nu } = \overline{\rho }_{\nu , 0} a^{-3}$|⁠.

CODE VALIDATION

Convergence tests

In this section, we check the convergence of our method with respect to particle number and initial redshift.

Fig. 2 shows the impact of changing the initial redshift from z_i = 49 to 99 for neutrinos of mass m_ν = 0.05 eV. The change at low redshift is of the order of 0.5 per cent and is consistent with slightly more small-scale non-linear growth for the simulation with the higher initial redshift, as expected from a pure CDM simulation.

Figure 2.

Open in new tab Download slide

Change in the total matter power spectrum between S05 (initial redshift z_i = 49) and S05IC (initial redshift z_i = 99). Positive values correspond to more power in S05. Each line represents a different redshift: z = 9 (green dot–dashed), z = 3 (grey dotted), z = 1 (blue dashed) and z = 0 (black solid).

Fig. 3 shows the effect of changing the particle number, comparing S10 to S10P, which has eight times more particles. Changing the number of particles causes the sampling of initial structure to change slightly, introducing sample variance. To minimize this effect, we show the ratio of the suppression due to neutrinos for pairs of simulations with different particle numbers, that is, (P_S10/P_S00)/(P_S10P/P_S00P) − 1. Our results are converged at the 1 per cent level for k < 10 h⁻¹ Mpc, and at the 0.2 per cent level for k < 1 h⁻¹ Mpc. The effect of increasing particle resolution is to increase the overall growth of non-linear structure, and as a consequence to enhance the relative suppression due to massive neutrinos.

Figure 3.

$Effect of the number of CDM particles on the power suppression by neutrinos. What is plotted is the fractional difference between the ratio (PS10/PS00) (both using 5123 CDM particles) and the ratio (PS10P/PS00P) (both using 10243 CDM particles). Positive values correspond to a larger suppression of power in the 10243-particle simulation. Each line represents a different redshift: z = 9 (green dot–dashed), z = 3 (grey dotted), z = 1 (blue dashed) and z = 0 (black solid).$

Open in new tab Download slide

Effect of the number of CDM particles on the power suppression by neutrinos. What is plotted is the fractional difference between the ratio (P_S10/P_S00) (both using 512³ CDM particles) and the ratio (P_S10P/P_S00P) (both using 1024³ CDM particles). Positive values correspond to a larger suppression of power in the 1024³-particle simulation. Each line represents a different redshift: z = 9 (green dot–dashed), z = 3 (grey dotted), z = 1 (blue dashed) and z = 0 (black solid).

Total power spectra

Fig. 4 shows the suppression in the total matter power spectrum, defined as the ratio between the total matter power spectrum including massive neutrinos with M_ν = ∑m_ν = 0.3 eV and the matter power spectrum with massless neutrinos, but unchanged Ω_M. The particle and Fourier neutrino simulation methods agree extremely well, and we see again the enhancement in the suppression due to the delayed onset of non-linear growth, as discussed in Bird et al. (2012).

Figure 4.

The suppression in the total power spectrum at z = 0. Solid lines represent a 256 Mpc h−1 with our Fourier-based method (blue) and the particle method (red). The dashed line represents the prediction of linear theory. The spike at k ≈ 7 h Mpc−1 (corresponding to half the grid frequency) is a numerical artefact.

Open in new tab Download slide

The suppression in the total power spectrum at z = 0. Solid lines represent a 256 Mpc h⁻¹ with our Fourier-based method (blue) and the particle method (red). The dashed line represents the prediction of linear theory. The spike at k ≈ 7 h Mpc⁻¹ (corresponding to half the grid frequency) is a numerical artefact.

Fig. 5 shows the percentage difference in the total power spectrum between the particle and semilinear neutrino simulations, for redshifts between z = 3 and 0 and for a variety of neutrino masses. The agreement is extremely good in all cases. We show those scales where our simulation is resolved to less than 1 per cent; note that the difference between the two methods is generally smaller than the convergence error shown in Figs 2 and 3. The effect of non-linear growth in the neutrino component, neglected in our semilinear method, only begins to become important for M_ν = 1.2 eV, at z < 1. Since neutrinos of this mass are already ruled out by current data, this is unlikely to be a practical limitation.

Figure 5.

Open in new tab Download slide

Percentage differences in the total matter power spectrum between the semilinear Fourier method and particles. ΔP_t = P_Fourier − P_Particle, so positive values correspond to more power in the Fourier method. The simulations used were S05 (top-left), S10 (top-right), S20 (bottom-left) and S40 (bottom-right). Each line represents a different redshift: z = 9 (green dot–dashed), z = 3 (grey dotted), z = 1 (blue dashed) and z = 0 (black solid).

One interesting feature is that on small scales we seem to see (slightly) more power in the Fourier method than the particle method, and the effect increases slightly at high redshift. This is not due to shot noise, which would lead to increased power in the particle method. This appears to be a convergence effect; as we have shown in Section 5.1, these scales are not converged to less than 1 per cent, and increasing the resolution of the simulation leads to a marginal increase of power on small scales. We therefore suspect that our semilinear method has an effective resolution slightly higher than the particle method, although if we computed the tree force for the neutrino particles, the situation would probably be reversed.

It is interesting to compare our results to those obtained using the original fully linear Fourier-space neutrino method of Brandbyge et al. (2008). Fig. 1 of that work is comparable to our Fig. 5. Our semilinear method agrees with the particle-based method slightly better in the linear regime. This is due to our derivation of the neutrino power spectrum from the CDM, which incorporates sample variance in the neutrino component in a way very similar to the particle method. In the non-linear regime, the agreement between the semilinear method and the particle method is significantly better than that between the fully linear and particle methods. The fully linear method shows differences from the particle method of >1 per cent for neutrinos with M_ν = 0.6 eV, and has an error of about 5 per cent at 1.2 eV. We ran test simulations using the fully linear theory power spectrum, but with our improved model for the neutrino phase structure (i.e. assuming that the neutrinos are always in phase with the CDM rather than keeping their initial phases). The result was still in good agreement with Brandbyge et al. (2008), showing that the improved behaviour of the semilinear method in the non-linear regime does indeed result from accounting for the deeper non-linear potential wells sourcing the linear growth of neutrinos.

Bird et al. (2012), using the linear theory neutrino implementation of Viel et al. (2010), found somewhat larger differences between the linear and particle methods than Brandbyge et al. (2008). When preparing this work, we discovered that this was due to a slight inconsistency in the algorithm used; the Fourier-space neutrino implementation was altering the initial CDM transfer function to include neutrinos, adding the neutrino power when calculating long-range forces, but neglecting them when outputting the power spectrum. Once this was corrected, we obtain fully linear Fourier-space neutrino results in good agreement with those of Brandbyge et al. (2008).

Neutrino power spectra

Fig. 6 shows the power spectra for the neutrino component at two different redshifts for the different methods. In the linear regime both codes agree well with linear theory. On non-linear scales our semilinear method shows additional power over the fully linear neutrino power spectrum. On all scales, P_CDM ≫ P_ν and k³P_ν/(2π²) ≪ 1.

Figure 6.

Open in new tab Download slide

The power spectrum of the neutrino component at (left) z = 0 and (right) z = 1. Solid black represents the results of our semilinear Fourier-based method from simulation S10. Dashed red represents the particle method, from S10P, to obtain lower shot noise. Dotted green represents pure linear theory. The vertical dashed grey line represents the approximate non-linear scale for the dark matter.

At z = 0, the particle method produces additional power in the neutrino component. It begins to deviate from our semilinear Fourier-based method at k = k_nl, the non-linear scale for the dark matter. This first becomes apparent for z < 0.5, for all M_ν ≥ 0.1 eV. We have checked that this area of the power spectrum is unchanged with a 512³-particle simulation, verifying that it is not being affected by shot noise. We ascribe this effect to non-linear growth in the neutrinos around the centres of clusters, as described in Ringwald & Wong (2004), Villaescusa-Navarro et al. (2011), and as we discussed in Section 2.4.

This does not appear to have any effect on the total matter power spectrum. Most of the effect of neutrinos on the total matter power spectrum comes from their time-integrated evolution; in particular, the size of the trough is due to their effect on the onset of non-linear growth. Since the highly non-linear effects on the neutrino power spectrum only take place at relatively late times, the impact on the total power spectrum is minimal. Moreover, on scales much smaller than the free-streaming scale, the main effect of including neutrinos is essentially to reduce the matter overdensity by a factor (1 − f_ν), while keeping the same background expansion rate (Lesgourgues & Pastor 2006). As long as |δ_ν| < |δ_cdm|, the exact value of the neutrino overdensity does not matter very much. The same conclusion was reached by Shoji & Komatsu (2009), who compared their third-order perturbation theory calculations (for both CDM and neutrinos, the latter approximated as a perfect fluid with pressure) to a computation where CDM is treated to third order but neutrinos are computed to linear order (without accounting for non-linear growth of potentials), similar to the treatment of Saito et al. (2008). They found that although neutrino overdensities can be underestimated by order unity on non-linear scales, the total matter power spectrum is very accurately obtained, even with the simpler method.⁶

Finally, let us point out that for z > 0.5, subtracting the scale-free shot noise from P_ν in the particle simulation produces very good agreement with the semilinear Fourier method, even at scales where the shot noise dominates. This is further evidence that neutrino shot noise is not having a strong dynamical effect, and is not causing spurious clustering.

Performance

Simulations S05–S20 were consistently faster when using our Fourier method. The speed increase was 13 per cent of the total walltime (which includes time spent reading and writing to disc). Note that the slowest single algorithm in gadget is the tree method for computing short-range forces, which is disabled for neutrinos even for our particle-based simulations, hence a large proportion of the execution time is independent of the method used to simulate neutrinos. More importantly, the total memory usage of gadget was 40 per cent smaller in the Fourier method than with particle neutrinos, essentially identical to the memory usage of a pure dark matter simulation. This is important because memory is often the limiting factor when performing large modern simulations.

The S10P simulation, which had eight times more dark matter particles than S10, took twelve times longer. This scaling is similar to that expected for a pure dark matter simulation, demonstrating that our neutrino method scales well. In fact, the only limit to scalability in the neutrino calculation is the need for interprocess communication when computing the power spectrum.

Overall, our Fourier method appears to have similar performance characteristics to a pure dark matter simulation, as should be expected; the time to compute the neutrino power spectrum is completely negligible compared to the N-body algorithms, and the most costly part of our Fourier algorithm is summing modes on the Fourier-transformed density grid to compute the power spectrum.

APPLICATIONS

Lyman α forest

The Lyman α forest is an indirect probe of the matter power spectrum at small, non-linear scales (k = 0.1-4 h^{− 1} Mpc), and at high redshift, z = 2-4. The power spectrum of the Lyman α flux measures the clustering of the absorption signal from neutral hydrogen in quasar spectra, and can be used to place constraints on the amplitude of primordial perturbations. When combined with constraints from large scales, this can lead to tight constraints on neutrino mass (Seljak 2000; Gratton et al. 2008; Viel et al. 2010). At these high redshifts, shot noise could be an issue for light neutrinos, so it is a natural place to apply our method. In addition, simulations with neutrino particles, dark matter and baryons can become unwieldy.

We used simulation F10BA to simulate the Lyman α forest at z = 2-4. Since the Lyman α forest is not the main focus of our paper, we shall not explain in detail our simulation methodology, and refer the interested reader to Viel et al. (2010). Our simulation calculates the clustering and ionization state of the gas at z = 2 using smooth particle hydrodynamics, together with radiative cooling and a reaction network which assumes optically thin gas in ionization equilibrium. We then simulate the flux in a set of 16 000 quasar spectra by calculating the absorption along random lines of sight. The averaged power spectrum of this flux is the observable quantity. Although constraints are available for k ≲ 4 h⁻¹ Mpc, the Lyman α forest is also sensitive to the collapse scale of hydrogen absorbers k ∼ 60 h⁻¹ Mpc, so very high resolution simulations are required. Fig. 7 shows the change in the flux power spectrum between our two methods; clearly, the differences are quite small. For simulations with the resolution of F10BA and m_ν = 0.1 eV, therefore, shot noise is not having a significant effect on the dynamics, even at z = 4, as also found by Viel et al. (2010). The total effect of neutrinos on the flux power is 5-10 per cent.

Figure 7.

Open in new tab Download slide

The change in the flux power spectrum between our Fourier method and the particle method, using simulation F10BA with m_ν = 0.1 eV. Positive values correspond to more power with the Fourier method. Each line represents a different redshift z = 4 (blue dotted), z = 3 (green dashed), and z = 2 (black solid). SDSS Lyman α data constrain the flux power spectrum for k ≲ 4 h⁻¹ Mpc (McDonald et al. 2006).

Hierarchy

Achieving maximal accuracy for small neutrino masses requires us to account for the neutrino hierarchy. This is difficult in particle-based codes, as two neutrino species doubles the amount of memory required, but has been done (Wagner, Verde & Jimenez 2012). It is, however, essentially trivial using our method. We have performed two simulations (S03IH and S03NH) with a total neutrino mass of M_ν = 0.1 eV using either the normal or inverted hierarchies, both to demonstrate the technique and to examine the effects of the hierarchy. For completeness, the masses of the three neutrino species were 0.022, 0.024 and 0.054 eV for S03NH, and 0, 0.049 and 0.051 eV for S03IH. Our results are shown in Fig. 8, together with the prediction of linear theory. The linear theory change due to the hierarchy is quite small (below our convergence error). Furthermore, the effect in linear theory is mostly due to the change in background evolution caused by one of the possible hierarchies having one massless species. As shown in fig. 16 of Lesgourgues & Pastor (2006), the difference between the hierarchies when all three neutrino species are massive in both cases is much less. Although differences in our cosmological parameters mean that the linear effect is smaller than that found by Wagner et al. (2012), we still find, as they did, an enhancement at non-linear scales, analogous to the enhancement for the total neutrino effect shown in Fig. 4. Determining the hierarchy separately from the total neutrino mass from an effect this subtle is likely to prove challenging. However, obtaining robust constraints with M_ν ≲ 0.1 eV will require accounting for the mass splitting.

Figure 8.

Open in new tab Download slide

The difference in the total matter power spectrum between the inverted and normal neutrino hierarchies. Negative values correspond to more power in the normal scenario. Each line represents a different redshift: z = 3 (grey dotted), z = 1 (blue dashed) and z = 0 (black solid). The orange dot–dashed line represents the linear prediction at z = 0.

21 cm forest

Observations of the spin-flip transition of neutral hydrogen in the 21 cm forest have the potential to probe structure growth at z ∼ 20 (McQuinn et al. 2006). Structure formation at these redshifts is sensitive to the non-linear growth of the first, small haloes. No one has yet examined the impact of neutrinos on the 21 cm forest, and performing the requisite simulations is beyond the scope of this paper. However, our code would be ideal for such a study. Haloes at these redshifts are not yet sufficiently massive for our method to break down in their central region, so one could probe arbitrarily small scales. Furthermore, shot noise would be a severe issue for the particle method at these redshifts.

DISCUSSION AND CONCLUSIONS

We have presented a new semilinear method for simulating the effects of massive neutrinos on cosmic structure, taking advantage of the hierarchy between the neutrino free-streaming scale and the non-linear scale, k_fs < k_nl for currently allowed neutrino masses. The power spectrum of the neutrino component is calculated using perturbation theory, and added to the CDM in Fourier space, as in Brandbyge et al. (2008). However, we improve upon their method by sourcing neutrino perturbations using the gravitational potential obtained from the full non-linear dark matter overdensity rather than from the dark matter power predicted by linear theory. We evolve the CDM, baryons and neutrinos simultaneously and self-consistently.

For observationally relevant neutrino masses, our method gives results for the non-linear power spectrum essentially identical to a fully converged neutrino particle treatment, a significant increase in accuracy over fully linear Fourier-space neutrinos. Our method also has several advantages over a particle treatment; it is faster, uses significantly less memory and does not suffer from shot noise in the neutrino component. Furthermore, it allows for an easy inclusion of the correct background evolution with relativistic corrections, and makes inclusion of the neutrino hierarchy trivial. These properties make it especially suitable for simulating neutrinos at the low end of the currently allowed mass range.

Our treatment is not strictly valid in the inner regions of very massive haloes, where the escape velocity might be larger than the thermal velocity of neutrinos, which may be captured and cluster non-linearly. Accurately describing the distribution of neutrinos in these haloes requires treating them as N-body particles. Because of this effect, we do not accurately recover the non-linear neutrino power spectrum at redshifts z < 0.5. However, this does not affect the total matter power spectrum, because on the one hand most of the neutrino effect is imprinted at earlier times, where our method is extremely accurate, and on the other hand, even with an enhanced clustering due to non-linear growth, neutrino overdensities are very subdominant compared to those of the CDM.

While we have focused on neutrinos, our treatment can also be easily applied to any HDM particle whose characteristic velocity is larger than (or at least of the order of) the escape velocity of massive clusters, and whose free-streaming length today is larger than the non-linear scale. There are many applications of our method beyond the non-linear matter power spectrum; we have briefly discussed the Lyman α and 21 cm forests. It allows the inclusion of neutrinos in an N-body simulation at minimal cost and effort, and should thus be useful for any problem in non-linear structure formation.

We thank Volker Springel for allowing one of us (SB) to use the code gadget-3, Matteo Viel for useful discussions and Daniel Grin for enlightening conversations on various aspects of this work. We are grateful to Anthony Challinor, Julien Lesgourgues and Eiichiro Komatsu for providing detailed comments on the draft of this paper, and acknowledge conversations with Matias Zaldarriaga, Scott Tremaine, Marilena LoVerde and Tobias Heinemann. The authors are supported by the National Science Foundation grants AST-080744 (for YA) and AST-0907969 (for SB).

1

NOVA proposal: http://www-nova.fnal.gov/.

2

The shot noise issue arises for neutrinos for two reasons. First, their large random velocities effectively randomly distribute them over the box, which leads to a Poisson spectrum |$P(k) = 1/\overline{n}$| at all scales. Secondly, their intrinsic clustering is very small on small scales. At large enough redshifts and for small enough scales, the true power may be completely swamped by the Poisson noise.

3

Note that consistently accounting for relativistic corrections in the evolution of neutrinos would also require accounting for CMB inhomogeneities sourcing the gravitational potentials; these are always (rightfully so) neglected in N-body simulations.

4

Interestingly, for a relativistic Fermi–Dirac distribution, the small-scale overdensity obtained from the linearized version of equation (49) is larger than the one obtained from the full non-linear equation.

5

Our initial conditions generator is freely available at http://github.com/sbird/S-GenIC

6

In the notation of Shoji & Komatsu (2009), our method assumes P_tot = f²_cP_{∞, c} + (2f_cf_νg₁ + f²_νg²₁)P_{∞, c}, which is better than the treatment of Saito et al. (2008), in the sense that we use the full non-linear CDM power spectrum as a source for neutrino overdensities.

REFERENCES

Abazajian

K.

,

Switzer

E. R.

,

Dodelson

S.

,

Heitmann

K.

,

Habib

S.

. ,

Phys. Rev. D

,

2005

, vol.

71

pg.

043507

Crossref

Search ADS

Abazajian

K. N.

, et al. ,

Astropart. Phys.

,

2011

, vol.

35

pg.

177

Crossref

Search ADS

Beaulieu

J. P.

, et al.

Coudé du Foresto

V.

,

Gelino

D. M.

,

Ribas

I.

. ,

ASP Conf. Ser. Vol. 430, Pathways Towards Habitable Planets

,

2010

San Francisco

Astron. Soc. Pac.

pg.

266

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Bertschinger

E.

,

Watts

P. N.

. ,

ApJ

,

1988

, vol.

328

pg.

23

Crossref

Search ADS

Bird

S.

,

Viel

M.

,

Haehnelt

M. G.

. ,

MNRAS

,

2012

pg.

420, 2175

Bond

J. R.

,

Szalay

A. S.

. ,

ApJ

,

1983

, vol.

274

pg.

443

Crossref

Search ADS

Bond

J. R.

,

Efstathiou

G.

,

Silk

J.

. ,

Phys. Rev. Lett.

,

1980

, vol.

45

pg.

1980

Crossref

Search ADS

Brandbyge

J.

,

Hannestad

S.

. ,

J. Cosmol. Astropart. Phys.

,

2009

, vol.

5

pg.

2

Crossref

Search ADS

Brandbyge

J.

,

Hannestad

S.

. ,

J. Cosmol. Astropart. Phys.

,

2010

, vol.

1

pg.

21

Crossref

Search ADS

Brandbyge

J.

,

Hannestad

S.

,

Haugbølle

T.

,

Thomsen

B.

. ,

J. Cosmol. Astropart. Phys.

,

2008

, vol.

8

pg.

20

Crossref

Search ADS

Brandbyge

J.

,

Hannestad

S.

,

Haugbølle

T.

,

Wong

Y. Y. Y.

. ,

J. Cosmol. Astropart. Phys.

,

2010

, vol.

9

pg.

14

Crossref

Search ADS

Brandenberger

R.

,

Kaiser

N.

,

Turok

N.

. ,

Phys. Rev. D

,

1987

, vol.

36

pg.

2242

Crossref

Search ADS

Carbone

C.

,

Verde

L.

,

Wang

Y.

,

Cimatti

A.

. ,

J. Cosmol. Astropart. Phys.

,

2011

, vol.

3

pg.

30

Crossref

Search ADS

Cooray

A. R.

. ,

A&A

,

1999

, vol.

348

pg.

31

de Putter

R.

, et al. ,

2012

arXiv:e-prints 1201.1909

de Vega

H. J.

,

Sanchez

N. G.

. ,

Phys. Rev. D

,

2012a

, vol.

85

pg.

043516

Crossref

Search ADS

de Vega

H. J.

,

Sanchez

N. G.

. ,

Phys. Rev. D

,

2012b

, vol.

85

pg.

043517

Crossref

Search ADS

Dehnen

W.

,

Read

J. I.

. ,

Eur. Phys. J. Plus

,

2011

, vol.

126

pg.

55

Crossref

Search ADS

Eitel

K.

. ,

Nucl. Phys. B

,

2005

, vol.

143

pg.

197

Crossref

Search ADS

Fogli

G. L.

,

Lisi

E.

,

Marrone

A.

,

Montanino

D.

,

Palazzo

A.

,

Rotunno

A. M.

. ,

Phys. Rev. D

,

2012

, vol.

86

pg.

013012

Crossref

Search ADS

Gilbert

I. H.

. ,

ApJ

,

1966

, vol.

144

pg.

233

Crossref

Search ADS

Gnedin

N. Y.

,

Hui

L.

. ,

MNRAS

,

1998

, vol.

296

pg.

44

Crossref

Search ADS

Gratton

S.

,

Lewis

A.

,

Efstathiou

G.

. ,

Phys. Rev.

,

2008

, vol.

D77

pg.

083507

Hall

A. C.

,

Challinor

A.

. ,

MNRAS

,

2012

, vol.

425

pg.

3504

Hannestad

S.

,

Haugbølle

T.

,

Schultz

C.

. ,

J. Cosmol. Astropart. Phys.

,

2012

, vol.

2

pg.

45

Crossref

Search ADS

Ichiki

K.

,

Takada

M.

,

Takahashi

T.

. ,

Phys. Rev. D

,

2009

, vol.

79

pg.

023520

Crossref

Search ADS

Jimenez

R.

,

Kitching

T.

,

Peña-Garay

C.

,

Verde

L.

. ,

J. Cosmol. Astropart. Phys.

,

2010

, vol.

5

pg.

35

Crossref

Search ADS

Kaplinghat

M.

,

Knox

L.

,

Song

Y.-S.

. ,

Phys. Rev. Lett.

,

2003

, vol.

91

pg.

241301

Crossref

Search ADS

PubMed

Komatsu

E.

, et al. ,

ApJS

,

2011

, vol.

192

pg.

18

Crossref

Search ADS

Kraus

C.

, et al. ,

Nucl. Phys. B

,

2005

, vol.

143

pg.

499

Crossref

Search ADS

Lesgourgues

J.

,

Pastor

S.

. ,

Phys. Rep.

,

2006

, vol.

429

pg.

307

Crossref

Search ADS

Lesgourgues

J.

,

Tram

T.

. ,

J. Cosmol. Astropart. Phys.

,

2011

, vol.

9

pg.

32

Crossref

Search ADS

Lesgourgues

J.

,

Matarrese

S.

,

Pietroni

M.

,

Riotto

A.

. ,

J. Cosmol. Astropart. Phys.

,

2009

, vol.

6

pg.

17

Crossref

Search ADS

Lewis

A.

,

Challinor

A.

. ,

Phys. Rev. D

,

2002

, vol.

66

pg.

023531

Crossref

Search ADS

Ma

C.-P.

,

Bertschinger

E.

. ,

ApJ

,

1995

, vol.

455

pg.

7

Crossref

Search ADS

Mao

Y.

,

Tegmark

M.

,

McQuinn

M.

,

Zaldarriaga

M.

,

Zahn

O.

. ,

Phys. Rev. D

,

2008

, vol.

78

pg.

023529

Crossref

Search ADS

Marulli

F.

,

Carbone

C.

,

Viel

M.

,

Moscardini

L.

,

Cimatti

A.

. ,

MNRAS

,

2011

, vol.

418

pg.

346

Crossref

Search ADS

McDonald

P.

, et al. ,

ApJS

,

2006

, vol.

163

pg.

80

Crossref

Search ADS

McQuinn

M.

,

Zahn

O.

,

Zaldarriaga

M.

,

Hernquist

L.

,

Furlanetto

S. R.

. ,

ApJ

,

2006

, vol.

653

pg.

815

Crossref

Search ADS

Planck Collaboration

. ,

2005

ESA-SCI

Ringwald

A.

,

Wong

Y. Y. Y.

. ,

J. Cosmol. Astropart. Phys.

,

2004

, vol.

12

pg.

5

Crossref

Search ADS

Saito

S.

,

Takada

M.

,

Taruya

A.

. ,

Phys. Rev. Lett.

,

2008

, vol.

100

pg.

191301

Crossref

Search ADS

PubMed

Saito

S.

,

Takada

M.

,

Taruya

A.

. ,

Phys. Rev. D

,

2009

, vol.

80

pg.

083528

Crossref

Search ADS

Schlegel

D.

,

White

M.

,

Eisenstein

D.

. ,

in Astro2010, A&A Decadal Survey. The Baryon Oscillation Spectroscopic Survey: Precision Measurement of the Absolute Cosmic Distance Scale

,

2009

National Academies Press

pg.

314

p

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Scoccimarro

R.

. ,

MNRAS

,

1998

, vol.

299

pg.

1097

Crossref

Search ADS

Seljak

U.

. ,

MNRAS

,

2000

, vol.

318

pg.

203

Crossref

Search ADS

Seljak

U.

,

Slosar

A.

,

McDonald

P.

. ,

J. Cosmol. Astropart. Phys.

,

2006

, vol.

10

pg.

14

Crossref

Search ADS

Shoji

M.

,

Komatsu

E.

. ,

ApJ

,

2009

, vol.

700

pg.

705

Crossref

Search ADS

Shoji

M.

,

Komatsu

E.

. ,

Phys. Rev. D

,

2010

, vol.

81

pg.

123516

Crossref

Search ADS

Singh

S.

,

Ma

C.-P.

. ,

Phys. Rev. D

,

2003

, vol.

67

pg.

023506

Crossref

Search ADS

Springel

V.

. ,

MNRAS

,

2005

, vol.

364

pg.

1105

Crossref

Search ADS

Takada

M.

,

Komatsu

E.

,

Futamase

T.

. ,

Phys. Rev. D

,

2006

, vol.

73

pg.

083520

Crossref

Search ADS

Vallinotto

A.

,

Das

S.

,

Spergel

D. N.

,

Viel

M.

. ,

Phys. Rev. Lett., 103, 091304

,

2009

Viel

M.

,

Haehnelt

M. G.

,

Springel

V.

. ,

J. Cosmol. Astropart. Phys.

,

2010

, vol.

6

pg.

15

Crossref

Search ADS

Vikhlinin

A.

, et al. ,

ApJ

,

2009

, vol.

692

pg.

1060

Crossref

Search ADS

Villaescusa-Navarro

F.

,

Miralda-Escudé

J.

,

Peña-Garay

C.

,

Quilis

V.

. ,

J. Cosmol. Astropart. Phys.

,

2011

, vol.

6

pg.

27

Crossref

Search ADS

Wagner

C.

,

Verde

L.

,

Jimenez

R.

. ,

ApJ

,

2012

, vol.

752

pg.

L31

Crossref

Search ADS

Wang

S.

,

Haiman

Z.

,

Hu

W.

,

Khoury

J.

,

May

M.

. ,

Phys. Rev. Lett.

,

2005

, vol.

95

pg.

011302

Crossref

Search ADS

PubMed

Weinberg

S.

. ,

Cosmology

,

2008

Oxford Univ. Press

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Weisstein

E. W.

. ,

Legendre Polynomial

,

2012a

MathWorld

A Wolfram Web Resource

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Weisstein

E. W.

. ,

Spherical Bessel Function of the First Kind

,

2012b

MathWorld

A Wolfram Web Resource

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Wong

Y. Y. Y.

. ,

J. Cosmol. Astropart. Phys.

,

2008

, vol.

10

pg.

35

Crossref

Search ADS

Wong

Y. Y. Y.

. ,

Annu. Rev. Nucl. Part. Sci.

,

2011

, vol.

61

pg.

69

Crossref

Search ADS

Xia

J.-Q.

, et al. ,

J. Cosmol. Astropart. Phys.

,

2012

, vol.

6

pg.

10

Crossref

Search ADS

APPENDIX A: ANALYTIC APPROXIMATION FOR LINEAR HDM CLUSTERING

In this section, we derive an approximate expression for δ_hdm valid on all scales during matter domination. We do not use this expression when evaluating neutrino clustering, but it provides some insight. Our main simplification here is to assume that the gravitational potential is nearly constant in time (or equivalently that the matter overdensity grows roughly linearly with the scale factor). During matter domination and for linear evolution, ϕ is in fact strictly constant. With this approximation, we get

\begin{eqnarray} \delta _{\rm hdm}(\boldsymbol k, \tau ) &\approx & \delta ^I_{\rm hdm}(\boldsymbol k, \tau _{\rm i}, \tau )\nonumber \\ && - k^2 \phi (\boldsymbol k, \tau ) \int _{0}^{s-s_{\rm i}} \mathcal {I}\left[\frac{k}{m}\tilde{s}\right] [a(s-\tilde{s})]^2 \tilde{s}\, {\rm d}\tilde{s}, \end{eqnarray}

(A1)

where we made the change of variables |$\tilde{s} \equiv s - s^{\prime }$|⁠. Assuming matter domination, we have a ∝ τ², τ = 2/(aH) and

\begin{equation} a(s - \tilde{s}) = \frac{a(s)}{\left(1 + \frac{1}{2} a^2(s) H(s) \tilde{s} \right)^2}. \end{equation}

(A2)

Changing variables to |$X = (k/m) \tilde{s}$|⁠, we obtain

\begin{eqnarray} &&\!\!\!\! \delta _{\rm hdm}(\boldsymbol k, \tau ) \approx \delta ^I_{\rm hdm}(\boldsymbol k, \tau _{\rm i}, \tau )\nonumber \\ && - (m a)^2 \phi (\boldsymbol k, \tau ) \int _{0}^{\frac{k(s-s_{\rm i})}{m}} \frac{X \mathcal {I}(X)}{\left( 1 + X/X_k \right)^4} {\rm d} X, \end{eqnarray}

(A3)

where X_k ≡ 2k/(ma²H) ∼ q^{− 1}₀k/k_fs(a). Let us now consider scales small enough that initial conditions are ‘forgotten’, kq₀(s − s_i)/m ≫ 1. In this case, δ^I ≈ 0 and the upper limit of the integral in the above equation can be approximated as infinity:

\begin{equation} \delta _{\rm hdm}(\boldsymbol k, \tau ) \approx - (m a)^2 \phi (\boldsymbol k, \tau ) \int _{0}^{+\infty } \frac{X \mathcal {I}(X)}{\left( 1 + X/X_k \right)^4} {\rm d} X. \end{equation}

(A4)

Note that during matter domination,

\begin{equation} q_0 X_k \approx \sqrt{\frac{a_i}{a}} \frac{k q_0(s - s_{\rm i})}{m}, \end{equation}

(A5)

so provided a/a_i is sufficiently large, q₀X_k may be of the order of unity, even if kq₀(s − s_i)/m ≫ 1. Using Poisson equation (39), we arrive at

\begin{equation} \delta _{\rm hdm}(\boldsymbol k, \tau ) \approx \mathcal {F}(k/k_{\rm fs}) \delta _{\rm M}(\boldsymbol k, \tau ), \end{equation}

(A6)

where we defined the dimensionless function

\begin{equation} \mathcal {F}(k/k_{\rm fs}) \equiv \frac{6}{X_k^2}\int _{0}^{+\infty } \frac{X \mathcal {I}(X)}{\left( 1 + X/X_k \right)^4} {\rm d} X. \end{equation}

(A7)

For k ≪ k_fs, q₀X_k → 0 and |$\mathcal {F}(k \ll k_{\rm fs})\rightarrow 1$|⁠, i.e. on large scales, δ_hdm ≈ δ_M. For k ≫ k_fs, q₀X_k → ∞ and one can show that

\begin{equation} \mathcal {F}(k \gg k_{\rm fs}) = \frac{6}{X_k^2}\overline{q^{-2}}, \end{equation}

(A8)

where

\begin{equation} \overline{q^{-2}} \equiv \frac{\int {\rm d}q f_0(q)}{\int {\rm d}q q^2 f_0(q)}. \end{equation}

(A9)

In this limit, we recover the small-scale solution of Section 3.3.2, |$\delta _{\rm hdm} = - \overline{v^{-2}} \phi = \left(k_{\rm fs}/k\right)^2 \delta _{\rm M}$|⁠.

We evaluated the function |$\mathcal {F}$| numerically for massive neutrinos and found that it is approximated to better than 12 per cent for any value of k/k_fs by the very simple form (see also Wong 2008)

\begin{equation} \mathcal {F}(k/k_{\rm fs}) \approx \frac{1}{\left(1 + k/k_{\rm fs}\right)^2}. \end{equation}

(A10)

APPENDIX B: ACCURACY OF THE ITERATIVE SOLUTION FOR THE NEUTRINO POWER SPECTRUM

In this section, we discuss the accuracy of the step where we obtain the neutrino power spectrum at time τ, given the total power spectrum at earlier times P_M(τ′ ≤ τ − Δτ), and the current CDM power spectrum, P_cdm(τ). We start by noticing that, because of our assumption of total correlation of neutrinos and CDM (equat‐ion 64), we have

\begin{equation} P_{\rm M}^{1/2}(k, \tau ) = (1 - f_{\nu }) P_{\rm cdm}^{1/2}(k, \tau ) + f_{\nu } P_{\nu }^{1/2}(k, \tau ). \end{equation}

(B1)

We now split equation (63) in a term that is completely known at the current time-step and a term that depends on the neutrino power spectrum between τ − Δτ and τ:

\begin{equation} P_{\nu }^{1/2}(k, \tau ) = P_{\nu }^{1/2}(k, \tau ; \textrm {known}) + \Delta P_{\nu }^{1/2}(\textrm {implicit}), \end{equation}

(B2)

where

\begin{eqnarray} \Delta P_{\nu }^{1/2}(\textrm {implicit}) \equiv \frac{3}{2} \Omega _{\rm M} H_0^2 f_{\nu } \int _{\tau - \Delta \tau }^{\tau } (s - s^{\prime })\mathcal {I}_{s^{\prime }, s} P^{1/2}_{\nu }(\tau ^{\prime }) \, \mathrm{d} \tau ^{\prime }. \nonumber\\ \end{eqnarray}

(B3)

We see that this renders the equation for P_ν(k, τ) implicit. Now, provided HΔt = aHΔτ ≪ 1, which is obviously the case in an N-body code, the neutrino power spectrum in the integral is approximately |$P_{\nu }(\tau ^{\prime }) = P_{\nu }(\tau )[1 + \mathcal {O}(H \Delta t)]$|⁠. Switching to the integration variable s^′, and using dτ′ = a(τ′)ds′ ≈ a(τ)ds′, we obtain, up to corrections of relative order |$\mathcal {O}(H \Delta t)$|⁠,

\begin{equation} \Delta P_{\nu }^{1/2}(\textrm {implicit}) = \epsilon (\tau ) P_{\nu }^{1/2}(k, \tau ), \end{equation}

(B4)

where we have defined

\begin{equation} \epsilon (\tau ) \equiv \frac{3}{2} \Omega _{\rm M} H_0^2 f_{\nu } a(\tau ) \int _{s-\Delta s}^{s} (s - s^{\prime }) \mathcal {I}_{s^{\prime },s} {\rm d} s^{\prime }, \end{equation}

(B5)

where Δs ≈ Δτ/a ≈ Δt/a². Now |$\mathcal {I} \le 1$| for any value of its argument, and therefore

\begin{eqnarray} \epsilon (\tau ) &\le & \frac{3}{4} \Omega _{\rm M} H_0^2 a^{-3} f_{\nu } (\Delta t)^2 = \frac{3}{4} f_{\nu } \Omega _{\rm M}(a) (H \Delta t)^2\nonumber \\ &<& f_{\nu } (H \Delta t)^2. \end{eqnarray}

(B6)

Therefore, whereas the exact solution is formally

\begin{equation} P_{\nu }^{1/2}(k, \tau ) = \frac{1}{1+ \epsilon (k, \tau )} P_{\nu }^{1/2}(k, \tau ; \textrm {known}), \end{equation}

(B7)

because |$\epsilon = \mathcal {O}\left(f_{\nu }(H \Delta t)^2\right)$| is small, an iterative solution will converge very rapidly, the error after n iterations being of the order of ϵⁿ [where the n = 0 iteration is P_ν = 0 and the (n + 1)-th iteration is obtained by using P_ν(implicit) obtained from the nth iteration]. Instead of initializing P_ν with 0, we moreover initialize it assuming P_ν(τ) ≈ P_ν(τ − Δτ) in the integral. The overall error is therefore extremely small, of the order of |$\mathcal {O}\left(f_{\nu }(H \Delta t)^3\right)$|⁠.

APPENDIX C: FITTING FUNCTION FOR |$\mathcal {I}(X)$|

The function |$\mathcal {I}(X)$| only depends on the dimensionless product x ≡ XT_{ν, 0} (we remind the reader that X was defined as an inverse comoving momentum). For neutrinos described with the relativistic Fermi–Dirac distribution, it can obtained from the following series expansion (Bertschinger & Watts 1988):

\begin{equation} \mathcal {I}(x) = \frac{4}{3 \zeta (3)}\sum _{n=1}^{\infty }(-1)^{n+1} \frac{n}{(n^2 + x^2)^2}. \end{equation}

(C1)

We also find that |$\mathcal {I}(x)$| is approximated to very high accuracy by the fit

\begin{equation} \mathcal {I}_{\rm fit}(x) \equiv \frac{1 + 0.0168 x^2 + 0.0407 x^4}{1 + 2.1734 x^2 + 1.6787 x^{4.1811} + 0.1467 x^8}. \end{equation}

(C2)

This fit gives the correct asymptotic behaviours

\begin{eqnarray} \mathcal {I}(x) &\approx & 1 - \frac{5 \zeta (5)}{2 \zeta (3)} x^2 \approx 1 - 2.1566 \ x^2, \ \ \ x \ll 1, \end{eqnarray}

(C3)

\begin{eqnarray} \mathcal {I}(x) &\approx & \frac{1}{3 \zeta (3) x^4} \approx \frac{0.2773}{x^4}, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x \gg 1. \end{eqnarray}

(C4)

The relative accuracy of the fit is |$|\Delta \mathcal {I}/\mathcal {I}| \lesssim 1$| per cent for 0 ≤ x ≤ 2 and |$|\Delta \mathcal {I}/\mathcal {I}| \lesssim 3$| per cent for 0 ≤ x ≤ +∞; the absolute accuracy is |$|\Delta \mathcal {I}| \lesssim 0.07$| per cent for 0 ≤ x ≤ +∞. We show the function |$\mathcal {I}$| in Fig. C1 We have checked explicitly that modifying our code to use the full function instead of the fitting formula does not affect our results.

Figure C1.

$Function $\mathcal {I}(x)$ that determines the damping of perturbations due to free-streaming of neutrinos (equation 38 with x = XTν, 0).$

Open in new tab Download slide

Function |$\mathcal {I}(x)$| that determines the damping of perturbations due to free-streaming of neutrinos (equation 38 with x = XT_{ν, 0}).

Download all slides

Month:	Total Views:
January 2017	8
February 2017	8
March 2017	12
April 2017	6
May 2017	6
June 2017	6
July 2017	4
August 2017	15
September 2017	5
October 2017	8
November 2017	15
December 2017	14
January 2018	10
February 2018	11
March 2018	8
April 2018	23
May 2018	25
June 2018	15
July 2018	18
August 2018	23
September 2018	4
October 2018	11
November 2018	16
December 2018	10
January 2019	2
February 2019	14
March 2019	19
April 2019	19
May 2019	29
June 2019	12
July 2019	9
August 2019	17
September 2019	9
October 2019	8
November 2019	18
December 2019	10
January 2020	13
February 2020	22
March 2020	19
April 2020	23
May 2020	13
June 2020	25
July 2020	15
August 2020	12
September 2020	10
October 2020	11
November 2020	15
December 2020	4
January 2021	7
February 2021	14
March 2021	8
April 2021	9
May 2021	19
June 2021	18
July 2021	15
August 2021	7
September 2021	13
October 2021	15
November 2021	22
December 2021	8
January 2022	13
February 2022	21
March 2022	6
April 2022	9
May 2022	11
June 2022	15
July 2022	14
August 2022	16
September 2022	24
October 2022	22
November 2022	10
December 2022	17
January 2023	12
February 2023	6
March 2023	11
April 2023	6
May 2023	10
June 2023	13
July 2023	10
August 2023	10
September 2023	12
October 2023	10
November 2023	9
December 2023	9
January 2024	12
February 2024	14
March 2024	20
April 2024	15

Article Contents

An efficient implementation of massive neutrinos in non-linear structure formation simulations

Abstract

INTRODUCTION

NOTATION, CHARACTERISTIC TIME-SCALES AND LENGTH-SCALES

Notation

Relativistic to non-relativistic transition

Free-streaming scale

Neutrino capture in massive haloes

LINEAR EVOLUTION OF HDM PERTURBATIONS IN A GENERAL BACKGROUND

Vlasov equation

Derivation for non-relativistic particles in an expanding universe

Linearization and Fourier transform

Integral solution

Discussion: limiting regimes

Large scales

Small scales

APPLICATION TO SIMULATING MASSIVE NEUTRINOS

Initial conditions

Phase evolution

Implementation

Simulation details

CODE VALIDATION

Convergence tests

Total power spectra

Neutrino power spectra

Performance

APPLICATIONS

Lyman α forest

Hierarchy

21 cm forest

DISCUSSION AND CONCLUSIONS

REFERENCES

APPENDIX A: ANALYTIC APPROXIMATION FOR LINEAR HDM CLUSTERING

APPENDIX B: ACCURACY OF THE ITERATIVE SOLUTION FOR THE NEUTRINO POWER SPECTRUM

APPENDIX C: FITTING FUNCTION FOR |$\mathcal {I}(X)$|

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only