Simulating the Cosmic Neutrino Background Using Collisionless Hydrodynamics

The cosmic neutrino background (CνB) has been robustly detected by observations of the cosmic microwave background perturbations (CMB; Planck Collaboration et al. 2018). At z ∼ 1100 when the primary CMB perturbations are set, the neutrinos are highly relativistic; however, observations of neutrino oscillations (de Salas et al. 2018) suggest that at least two species of neutrino are massive (with the heavier one having a mass of ${m}_{\nu }\gtrsim 0.05$ eV) and so they should be nonrelativistic today. Measuring the properties of the nonrelativistic CνB has become one of the chief goals of upcoming cosmological experiments. Of principle importance is a measurement of the total neutrino mass ${M}_{\nu }=\sum {m}_{\nu }$ that, in conjunction with oscillation experiments, may also resolve the neutrino hierarchy and determine whether any neutrinos are massless. Currently cosmological constraints (e.g., Planck Collaboration et al. 2018; Palanque-Delabrouille et al. 2020) of M_ν are significantly better than terrestrial ones (Aker et al. 2019).

The CνB is distinct in that it becomes nonrelativistic but still remains quite hot and so it does not cluster nearly as much on small scales as the dominant cold dark matter (CDM). The principle effect of this lack of clustering is a modulation of the matter power spectrum in a way that is sensitive to the neutrino energy density ${{\rm{\Omega }}}_{\nu }$. If standard cosmology holds, where neutrinos are the only hot species and have the standard decoupling temperature ${T}_{\nu }={(4/11)}^{1/3}{T}_{\mathrm{CMB}}$, then ${{\rm{\Omega }}}_{\nu }{h}^{2}\,={M}_{\nu }/93.14$ (Mangano et al. 2005). This technique for determining M_ν is a substantial challenge both because the modulation is small and because we principally measure not the underlying matter field but rather biased tracers of it such as galaxies. If M_ν is minimal, observations sensitive to the matter power spectrum by the Dark Energy Spectroscopic Instrument (DESI; Aghamousa et al. 2016), Euclid (Amendola et al. 2018) and CMB-S4 (Abazajian et al. 2016) are forecasted to detect the CνB at around 3σ (Brinckmann et al. 2019). Finding other observables sensitive to neutrinos is therefore critical to improve the robustness of the nonrelativistic CνB detection. Alternatives that have been suggested include void statistics (Massara et al. 2015; Banerjee & Dalal 2016a; Kreisch et al. 2019), the relative velocity effect (Zhu et al. 2014, 2016; Inman et al. 2015, 2017; Zhu & Castorina 2020), the scale-dependent halo bias (LoVerde 2014; Chiang et al. 2018, 2019; Banerjee et al. 2020), differential neutrino condensation (Yu et al. 2017), and via galaxy spins (Yu et al. 2019).

Because some of these new observables may be sensitive to nonlinear behavior of the neutrinos themselves, it is critical to accurately model the dynamical evolution of neutrinos in the large-scale structure. The CνB has been included in simulations in a number of ways. The most straightforward methods invoke linear theory, either interpolating between precomputed results (Brandbyge & Hannestad 2009) or integrating alongside the simulation (Ali-Haïmoud & Bird 2013). These methods are quite fast and lead to the correct modulation of the matter power spectrum but do not accurately resolve neutrino dynamics. The standard method to obtain correct neutrino behavior is to use N-body particles that have large thermal velocities drawn from the Fermi–Dirac distribution (e.g., Brandbyge et al. 2008; Viel et al. 2010; Zennaro et al. 2017). The chief downside of this style of simulation is that the random velocities are so large that they introduce significant randomness into the neutrino density field that substantially erases the true perturbations on small scales. The effects of this Poisson noise can be suppressed through the use of many particles (Emberson et al. 2017) or through the use of hybrid methods where part of the evolution is linear (Brandbyge & Hannestad 2010; Inman et al. 2015; Bird et al. 2018).

There are two methods employed in the literature that are Poisson noise-free. The first is to use neutrino particles but introduce regularity in the thermal velocities such that random clustering does not occur (Banerjee et al. 2018). In practice, this involves utilizing many particles per grid cell, with each cell having the same sampling of random velocities. The second is to instead solve the Boltzmann moment equations on a grid, which has no Poisson noise by construction. The difficulty here is that the Boltzmann moment equations are an infinite hierarchy of equations with no known truncation scheme consistent with the neutrino distribution. One approach is to close the hierarchy by utilizing linear theory for the stress terms in the momentum equation (Dakin et al. 2019), which works well provided that the pressure is always correlated with the density field. The other method is to use N-body particles for closure by using them to estimate the stress terms (instead of the density; Banerjee & Dalal 2016a).

In this paper we also solve the Boltzmann moment equations; however, we do so by using a closure scheme that does not rely on either external transfer functions or N-body particles. We do this by decomposing the neutrino phase space into shells of uniform speed, for which a straightforward linear closure scheme exists. It is conceptually quite similar to the multiflow description of Dupuy & Bernardeau (2014, 2015), although we consider particles with the same speed rather than the same velocity. We show slices comparing our method and standard N-body in Figure 1. The lack of Poisson noise is immediately evident. In Section 2 we motivate the closure scheme for the velocity shells. In Section 3 we discuss our numerical implementation of the method into the CUBE code (Yu et al. 2018). We then show results of these simulations, including comparisons to standard N-body methods, in Section 4. We discuss future optimizations and extensions in Sections 5 and 6.

Zoom In Zoom Out Reset image size

2.1. Vlasov Equation and Boltzmann Moment Equations

The phase space, $f(t,{x}^{i},{v}^{i})$, of collisionless particles is described by the Vlasov equation (also known as the collisionless Boltzmann equation):

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial t}+{v}^{i}\displaystyle \frac{\partial f}{\partial {x}^{i}}+{g}^{i}\displaystyle \frac{\partial f}{\partial {v}^{i}}=0,\end{eqnarray} \tag{ 1 }$$

where we employ the Einstein summation convention over the three spatial dimensions. In an expanding universe with scalefactor a, t is the Newtonian time (related to the physical time, t_p, via ${a}^{2}{dt}={{dt}}_{p}$), xⁱ is the comoving coordinate (related to physical coordinates, xⁱ_p, by ${{adx}}^{i}={{dx}}_{p}^{i}$), vⁱ is the conjugate velocity (${v}^{i}={{dx}}^{i}/{dt}$), and gⁱ is the gravitational acceleration field (${g}^{i}=-{a}^{2}\partial \phi /\partial {x}^{i}$ where ϕ is the Newtonian potential; Bertschinger 1995; Inman & Pen 2017). The initial conditions for the CνB is the relativistic Fermi–Dirac distribution:

$$\begin{eqnarray}&&{\bar{f}}_{\nu }=\displaystyle \frac{{\beta }^{3}}{6\pi \zeta (3)}\displaystyle \frac{1}{\exp [\beta v]+1},\end{eqnarray} \tag{ 2 }$$

where $\beta ={m}_{\nu }/{T}_{\nu }$ and the normalization is chosen such that the mean density is unity: $\bar{\rho }=\int {d}^{3}v\bar{f}(v)=1$.

Equation (1) can be converted into a hierarchy of moment equations. We denote coarse-grained fields with the notation $\langle A(t,{x}^{i},{v}^{i})\rangle =\int {d}^{3}{vAf}$. We will concern ourselves with the the density, $\rho =\langle 1\rangle$; the momentum, $\rho {V}^{i}=\langle {v}^{i}\rangle$; the stress–energy tensor $\rho {{\rm{\Sigma }}}^{{ij}}=\langle {v}^{i}{v}^{j}\rangle$; and the third order heat tensor, $\rho {Q}^{{ijk}}=\langle {v}^{i}{v}^{j}{v}^{k}\rangle$. The stress–energy tensor is often expanded as $\rho {{\rm{\Sigma }}}^{{ij}}=\rho {V}^{i}{V}^{j}+\rho {\sigma }^{{ij}}$, where Vⁱ is the velocity field and ${\sigma }^{{ij}}={\rho }^{-1}\langle ({v}^{i}-{V}^{i})({v}^{j}-{V}^{j})\rangle$ is the velocity dispersion. $\rho {Q}^{{ijk}}$ can be expanded as $\rho {Q}^{{ijk}}=\rho {V}^{i}{V}^{j}{V}^{k}+\rho {V}^{i}{\sigma }^{{jk}}$+$\rho {V}^{j}{\sigma }^{{ik}}\,+\rho {V}^{k}{\sigma }^{{ij}}+\rho {q}^{{ijk}}$, where $\rho {q}^{{ijk}}=\langle ({v}^{i}-{V}^{i})({v}^{j}-{V}^{j})({v}^{k}-{V}^{k})\rangle$ is the heat flux tensor with the associated heat flux $(1/2)\rho {q}^{{iik}}$.

By integrating Equation (1) over velocity, we can obtain their equations of motion. These are the continuity equation,

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\rho +\displaystyle \frac{\partial }{\partial {x}^{k}}\rho {V}^{k}=0;\end{eqnarray} \tag{ 3 }$$

the momentum equation,

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\rho {V}^{i}+\displaystyle \frac{\partial }{\partial {x}^{k}}\rho {{\rm{\Sigma }}}^{{ik}}=\rho {g}^{i};\end{eqnarray} \tag{ 4 }$$

and the stress–energy equation,

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\rho {{\rm{\Sigma }}}^{{ij}}+\displaystyle \frac{\partial }{\partial {x}^{k}}\rho {Q}^{{ijk}}=\rho {V}^{i}{g}^{j}+\rho {V}^{j}{g}^{i}.\end{eqnarray} \tag{ 5 }$$

2.2. Hydrodynamic Equations

We wish to re-express these equations in terms of hydrodynamic variables. To that end, we expand the velocity dispersion in terms of primitive variables $\rho {\sigma }^{{ij}}=P{\delta }^{{ij}}+{\pi }^{{ij}}$ where P is the pressure and ${\pi }^{{ij}}$ is the shear stress, which is traceless, ${\pi }^{{ii}}=0$, and symmetric, ${\pi }^{{ij}}={\pi }^{{ji}}$. We now require equations of motion for P and ${\pi }^{{ij}}$, which can be obtained from Equation (5). Instead of using the stresses directly, it is convenient to use the energy, E, and anisotropic stress, ${\tau }^{{ij}}$. The energy is related to the trace of the stress–energy tensor:

$$\begin{eqnarray}\begin{array}{rcl}E & = & \displaystyle \frac{1}{2}\rho {{\rm{\Sigma }}}^{{ii}}\\ & = & \displaystyle \frac{1}{2}\rho {V}^{2}+\displaystyle \frac{3}{2}P\\ & = & \displaystyle \frac{1}{2}\rho {V}^{2}+\displaystyle \frac{P}{\gamma -1},\end{array}\end{eqnarray} \tag{ 6 }$$

which matches the definition of the energy of an ideal fluid with γ = 5/3 (Mitchell et al. 2013). The anisotropic stress is then the traceless component:

$$\begin{eqnarray}\begin{array}{rcl}{\tau }^{{ij}} & = & \rho {{\rm{\Sigma }}}^{{ij}}-\displaystyle \frac{1}{3}\rho {{\rm{\Sigma }}}^{{kk}}{\delta }^{{ij}}\\ & = & \rho {V}^{i}{V}^{j}-\displaystyle \frac{1}{3}\rho {V}^{2}{\delta }^{{ij}}+{\pi }^{{ij}}.\end{array}\end{eqnarray} \tag{ 7 }$$

The equations of motion in terms of the conservative variables are then

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\rho {V}^{i}+\displaystyle \frac{\partial }{\partial {x}^{k}}\left[\rho {V}^{i}{V}^{k}+P{\delta }^{{ik}}+{\pi }^{{ik}}\right]=\rho {g}^{i},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}E+\displaystyle \frac{\partial }{\partial {x}^{k}}\left[{{EV}}^{k}+{{PV}}^{k}+{V}^{l}{\pi }^{{lk}}+\displaystyle \frac{1}{2}\rho {q}^{{iik}}\right]=\rho {V}^{k}{g}^{k},\end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial }{\partial t}{\tau }^{{ij}}+\displaystyle \frac{\partial }{\partial {x}^{k}}\left[{\tau }^{{ij}}{V}^{k}+{{PV}}^{i}{\delta }^{{jk}}+{{PV}}^{j}{\delta }^{{ik}}-\displaystyle \frac{2}{3}{{PV}}^{k}{\delta }^{{ij}}\right.\\ \ \ \left.+{V}^{i}{\pi }^{{jk}}+{V}^{j}{\pi }^{{ik}}-\displaystyle \frac{2}{3}{V}^{l}{\pi }^{{lk}}{\delta }^{{ij}}+\rho {q}^{{ijk}}-\displaystyle \frac{1}{3}\rho {q}^{{llk}}{\delta }^{{ij}}\right]\\ \ \ =\rho {V}^{i}{g}^{j}+\rho {V}^{j}{g}^{i}-\displaystyle \frac{2}{3}\rho {V}^{k}{g}^{k}{\delta }^{{ij}}.\end{array}\end{eqnarray} \tag{ 10 }$$

2.3. Closing the Hierarchy

Solving these equations is exceptionally challenging due to the nontrivial stresses and heat tensor that are present even in linear theory. Strategies that have been employed to estimate the higher moments include linear theory (Dakin et al. 2019), using N-body particles to estimate higher moments (Pueblas & Scoccimarro 2009; Banerjee & Dalal 2016b) or the "Zero Heat Flux" approximation, which sets $\rho {q}^{{ijk}}=0$ (Vorobyov & Theis 2006; Mitchell et al. 2013).

An alternate approach is to try to deduce an (potentially nonlinear) equation of state for the heat flux that depends only on the hydrodynamic variables (see, e.g., Hammett & Perkins 1990; Chust & Belmont 2006; Wang et al. 2015). To do this we will consider two test problems: the free evolution of particles and the linearized evolution of phase space.

2.3.1. Free Evolution of Particles

We first consider the free evolution (gⁱ = 0) of collisionless particles evolving from an initial distribution of $f({t}_{i})$. In this case the exact solution to Equation (1) is given in Fourier space by

$$\begin{eqnarray}&&\mathrm{FT}[f]=\mathrm{FT}[f({t}_{i})]\exp [-{{ik}}^{j}{v}^{j}(t-{t}_{i})].\end{eqnarray} \tag{ 11 }$$

If we now consider an initial distribution consisting solely of particles with the same velocity $f({t}_{i})\propto \rho ({t}_{i},{x}^{i})/(4\pi {v}^{2}){\delta }_{D}(v-u)$, the density field is given by

$$\begin{eqnarray}&&\mathrm{FT}[\rho ]=\mathrm{FT}[\rho ({t}_{i},{x}^{i})]{j}_{0}({ku}(t-{t}_{i})).\end{eqnarray} \tag{ 12 }$$

The inverse Fourier transform of ${j}_{0}({ku}(t-{t}_{i}))$ is ${\delta }_{D}(r\,-u(t-{t}_{i}))/(4\pi {r}^{2})$ and so this can be thought of as the convolution of the initial density field with the inverse square law. In general, Equation (12) satisfies the wavelike spherical Bessel equation in real space:

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}}{\partial {t}^{2}}\rho +\displaystyle \frac{2}{t-{t}_{i}}\displaystyle \frac{\partial }{\partial t}\rho -{u}^{2}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{i}\partial {x}^{i}}\rho =0.\end{eqnarray} \tag{ 13 }$$

After becoming nonrelativistic, the CνB has only very small perturbations, and so no direct convolution needs to be performed.

2.3.2. Linearized Vlasov Equation

We next consider the Vlasov equation linearized about some initial homogeneous and the isotropic velocity distribution $f({t}_{i},x,v)=\bar{f}(v)$:

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial t}+{v}^{i}\displaystyle \frac{\partial f}{\partial {x}^{i}}+{g}^{i}\displaystyle \frac{\partial \bar{f}}{\partial {v}^{i}}=0.\end{eqnarray} \tag{ 14 }$$

We first remark that "linearized" is a misnomer in our particular application: since neutrinos are a small component (${f}_{\nu }\ll 1$), the Vlasov equation is already approximately linear as gⁱ is approximately independent of neutrino density. Instead we are replacing a derivative term with a source term. This is precisely why the "linear response" method does not produce accurate neutrino power spectra despite neutrinos remaining linear (e.g., ${\delta }_{\nu }\ll 1$) at all times. Nonetheless, we will use the standard terminology throughout. The linear response solution to this equation is given in Fourier space by Gilbert (1966) and Ali-Haïmoud & Bird (2013):

$$\begin{eqnarray}&&\mathrm{FT}[f]={\int }_{{t}_{i}}^{t}{dt}^{\prime} (-\mathrm{FT}[{g}^{i}])\displaystyle \frac{\partial \bar{f}}{\partial {v}^{i}}\exp [-{{ik}}^{j}{v}^{j}(t-t^{\prime} )].\end{eqnarray} \tag{ 15 }$$

The moments are then given by

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{FT}[\rho ] & = & {\displaystyle \int }_{0}^{\infty }{dv}4\pi {v}^{2}\bar{f}(v)\\ & & \times \,{\displaystyle \int }_{{t}_{i}}^{t}{dt}^{\prime} (-{{ik}}^{j}\mathrm{FT}[{g}^{j}])(t-t^{\prime} ){j}_{0}({kv}(t-t^{\prime} )),\end{array}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}\begin{array}{rcl}{{ik}}^{i}\mathrm{FT}[\rho {V}^{i}] & = & {\displaystyle \int }_{0}^{\infty }{dv}4\pi {v}^{2}\bar{f}(v)\\ & & \times \,{\displaystyle \int }_{{t}_{i}}^{t}{dt}^{\prime} ({{ik}}^{j}\mathrm{FT}[{g}^{j}]\cos [{kv}(t-t^{\prime} )],\end{array}\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\begin{array}{l}{{ik}}^{i}{{ik}}^{j}\mathrm{FT}[\rho {{\rm{\Sigma }}}^{{ij}}]={\displaystyle \int }_{0}^{\infty }{dv}4\pi {v}^{2}\bar{f}(v){v}^{2}{k}^{2}\\ \ \ \times \,{\displaystyle \int }_{{t}_{i}}^{t}{dt}^{\prime} (-{{ik}}^{j}\mathrm{FT}[{g}^{j}])(t-t^{\prime} ){j}_{0}({kv}(t-t^{\prime} )),\end{array}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\begin{array}{l}{{ik}}^{i}{{ik}}^{j}{{ik}}^{k}\mathrm{FT}[\rho {Q}^{{ijk}}]={\displaystyle \int }_{0}^{\infty }{dv}4\pi {v}^{2}\bar{f}(v){v}^{2}{k}^{2}\\ \quad \times \,{\displaystyle \int }_{{t}_{i}}^{t}{dt}^{\prime} ({{ik}}^{j}\mathrm{FT}[{g}^{j}]\cos [{kv}(t-t^{\prime} )],\end{array}\end{eqnarray} \tag{ 19 }$$

where we have omitted the homogeneous terms. It is immediately clear that density fields can be composed as sums over velocity shells. That is, if the background distribution with ${\bar{f}}_{u}=1/(4\pi {v}^{2}){\delta }_{D}(v-u)$ has moments such as ${\rho }_{u}$, then a more general distribution function, $\bar{f}$, has $\rho ={\int }_{0}^{\infty }{du}4\pi {u}^{2}\bar{f}{\rho }_{u}$ (Inman & Pen 2017). Furthermore, these shells have simple higher moments. In particular, the stress–energy tensor can be written entirely as an isotropic pressure, $\rho {{\rm{\Sigma }}}^{{ij}}=\rho {\sigma }^{{ij}}=P{\delta }^{{ij}}=\rho {u}^{2}{\delta }^{{ij}}$, and so the density satisfies the wave equation:

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}}{\partial {t}^{2}}\rho -{u}^{2}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{i}\partial {x}^{i}}\rho =-\displaystyle \frac{\partial }{\partial {x}^{i}}{g}^{i}.\end{eqnarray} \tag{ 20 }$$

Substituting this relationship and the linearized Equation (3) into the linearized Equation (5) leads to

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial }{\partial {x}^{k}}{Q}^{{ijk}} & = & \displaystyle \frac{\partial }{\partial {x}^{k}}({u}^{2}{V}^{k}{\delta }^{{ij}})\\ & = & \displaystyle \frac{\partial }{\partial {x}^{k}}\left[{u}^{2}({V}^{i}{\delta }^{{jk}}+{V}^{j}{\delta }^{{ik}}+{V}^{k}{\delta }^{{ij}})+{q}^{{ijk}}\right],\end{array}\end{eqnarray} \tag{ 21 }$$

or

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial }{\partial {x}^{k}}{q}^{{ijk}} & = & \displaystyle \frac{\partial }{\partial {x}^{k}}\left[-{u}^{2}{V}^{i}{\delta }^{{jk}}-{u}^{2}{V}^{j}{\delta }^{{ik}}\right]\\ & = & -{u}^{2}\displaystyle \frac{\partial {V}^{i}}{\partial {x}^{j}}-{u}^{2}\displaystyle \frac{\partial {V}^{j}}{\partial {x}^{i}}.\end{array}\end{eqnarray} \tag{ 22 }$$

Linear velocity fields are curl free and can be written in terms of a velocity potential ${V}^{i}=\partial {\phi }_{v}/\partial {x}^{i}$, which leads to

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial {x}^{k}}{q}^{{ijk}}=-2{u}^{2}\displaystyle \frac{{\partial }^{2}{\phi }_{v}}{\partial {x}^{i}\partial {x}^{j}}.\end{eqnarray} \tag{ 23 }$$

The symmetric form of q^ijk satisfying this is then given by

$$\begin{eqnarray}&&{q}^{{ijk}}=-2{u}^{2}\displaystyle \frac{\partial }{\partial {x}^{i}}\displaystyle \frac{\partial }{\partial {x}^{j}}\displaystyle \frac{\partial }{\partial {x}^{k}}{{\rm{\nabla }}}^{-2}{\phi }_{v},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&=-\displaystyle \frac{2}{3}{u}^{2}{{\rm{\nabla }}}^{-2}\left[\displaystyle \frac{\partial }{\partial {x}^{i}}\displaystyle \frac{\partial }{\partial {x}^{j}}{V}^{k}+\displaystyle \frac{\partial }{\partial {x}^{i}}\displaystyle \frac{\partial }{\partial {x}^{k}}{V}^{j}+\displaystyle \frac{\partial }{\partial {x}^{j}}\displaystyle \frac{\partial }{\partial {x}^{k}}{V}^{i}\right],\end{eqnarray} \tag{ 25 }$$

where ∇⁻² denotes the inverse Laplacian.

2.3.3. Shell Equations of Motion

Unfortunately, it is unclear how to best generalize the linear heat flux of a shell once nonlinear evolution occurs. Simply using Equation (22) for the heat tensor leads to an energy flux of $(P-\bar{\rho }{u}^{2}){V}^{k}+{\pi }^{{ij}}$, which could be 0, the nonlinear fluctuation of P, or something else entirely. Lacking a nonlinear scheme, we instead opt for a numerically convenient closure choice that allows for the shell equations to be independent of u:

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial {x}^{k}}\rho {q}^{{ijk}}=\displaystyle \frac{\partial }{\partial {x}^{k}}\left(-\rho {\sigma }^{{ik}}{V}^{j}-\rho {\sigma }^{{jk}}{V}^{i}\right),\end{eqnarray} \tag{ 26 }$$

which yields the much simplified continuity equations for stress energy:

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}E+\displaystyle \frac{\partial }{\partial {x}^{k}}{{EV}}^{k}=\rho {V}^{k}{g}^{k},\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{\tau }^{{ij}}+\displaystyle \frac{\partial }{\partial {x}^{k}}{\tau }^{{ij}}{V}^{k}=\rho {V}^{i}{g}^{j}+\rho {V}^{j}{g}^{i}-\displaystyle \frac{2}{3}\rho {V}^{k}{g}^{k}{\delta }^{{ij}}.\end{eqnarray} \tag{ 28 }$$

The complete set of equations for a shell is therefore Equations (3), (8), (27), and (28). In the absence of gravity, these equations collapse down to wave dynamics. In this work we make one further simplification and neglect the shear stress entirely, ${\pi }^{{ij}}=0$, leading to ideal collisionless hydrodynamics. This simplification is still consistent with linear theory but allows us to save substantial resources by not solving Equation (28). We consider two alternative schemes, the isothermal and ideal gas approximations, in the Appendix.

Let us lastly comment on the difference between this approach and linear response. In linear response, approximations are made to the gravitational terms, i.e., $\rho {g}^{i}\to \bar{\rho }{g}^{i}$. In Equations (27) and (28), we have made approximations to the hydrodynamic fluxes but kept gravitational terms fully intact. We therefore expect this approach to work significantly better as it retains linearity in the neutrino perturbations themselves.

In this section we describe our implementation of the neutrino solver into the CUBE code (Yu et al. 2018). CUBE is a re-factored version of CUBEP³M (Harnois-Déraps et al. 2013) designed to use the minimum amount of memory possible: 6 bytes per particle. It features a three level force decomposition: a "coarse" particle mesh force over n_c cells where the only global fast Fourier transform (FFT) is performed, a "fine" particle mesh force with a resolution of ${n}_{f}=4{n}_{c}$ where local FFTs are performed, and a direct particle–particle (PP) force at the cell level. CUBE is parallelized with both OpenMP and co-array Fortran. Both CUBE and CUBEP³M are publicly available.

3.1. Single-shell Hydrodynamics

3.2. CνB

3.3. Initial Conditions

Since we do not have transfer functions for individual momenta (although in principle these can be extracted from Boltzmann codes), we opt for approximate initial conditions for each velocity shell. To determine an approximate solution for each shell, we first consider their linear equations of motion, Equation (20) (this can also be seen by differentiating Equation (16) twice):

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{\delta }_{u}}{\partial {t}^{2}}+{k}^{2}{u}^{2}{\delta }_{u}={a}^{2}(-{k}^{2}\phi ).\end{eqnarray} \tag{ 29 }$$

On large scales, where k is small, the second term is negligible and ${\delta }_{u}\simeq {\delta }_{\mathrm{CDM}}$. On small scales, the fluid responds instantaneously and the solution is ${\delta }_{u}\simeq {({k}_{u}/k)}^{2}{\delta }_{\mathrm{CDM}}$, where ${k}_{u}=\sqrt{(3/2){{\rm{\Omega }}}_{m}a}{H}_{0}/u$ is the free streaming scale of the shell (Inman & Pen 2017). We utilize the following approximation:

$$\begin{eqnarray}&&{\rho }_{u}=1+{\delta }_{u}=1+{\delta }_{\mathrm{CDM}}\displaystyle \frac{1}{1+\lambda k/{k}_{u}+{(k/{k}_{u})}^{2}}.\end{eqnarray} \tag{ 30 }$$

The corresponding velocity field is obtained via the continuity equation:

$$\begin{eqnarray}&&{V}_{u}^{i}={V}_{\mathrm{CDM}}^{i}\displaystyle \frac{1+1.5\lambda k/{k}_{u}+2{(k/{k}_{u})}^{2}}{{(1+\lambda k/{k}_{u}+{(k/{k}_{u})}^{2})}^{2}}.\end{eqnarray} \tag{ 31 }$$

In the simulations used in this paper we used λ = 1, which is appropriate for z_ν = 0; however, we have subsequently solved for the optimal value of λ and find that λ = 0.75, 0.76, 0.77, and 0.82 for m_ν = 0.05, 0.10, 0.20, and 0.40 at z_ν = 10. These values are are unchanged at z_ν = 5 but increase toward 1 at z_ν = 0. Utilizing a constant λ = 0.75 at either z_ν = 10 or 5 produces neutrino density and velocity transfer functions accurate to around 10% regardless of neutrino mass.

When computing the initial momentum density we opt to include the second order term, i.e., $\rho {V}^{k}=(1+{\delta }_{u}){V}_{u}^{k}$, rather than just V_u^k. Likewise the initial energy is therefore

$$\begin{eqnarray}&&{E}_{u}={\rho }_{u}\left[\displaystyle \frac{1}{2}{V}_{u}^{2}+\displaystyle \frac{{u}^{2}}{\gamma -1}\right]\end{eqnarray} \tag{ 32 }$$

with the kinetic term included.

We run simulations of neutrinos of masses 0.05, 0.10, 0.20, and 0.40 eV using both CUBE and CUBEP³M. We keep the initial Gaussian noise realization equivalent between the two simulations and use the same number of particles but otherwise allow them to run independently. Because the neutrinos are not coupled to the PP force, we run CUBE without the direct force. Importantly this means that the CUBEP³M simulations have better force resolution as they include the pair-wise force on small scales. For cosmological parameters, we match those of Cataneo et al. (2020) so that we can compare to the high-resolution results computed in that work. Specifically we use a flat universe with fixed $1-{{\rm{\Omega }}}_{{\rm{\Lambda }}}={{\rm{\Omega }}}_{m}={{\rm{\Omega }}}_{c}+{{\rm{\Omega }}}_{b}+{{\rm{\Omega }}}_{\nu }=0.291$, Ω_b = 0.047, and h = 0.6898 and implicitly varying Ω_c via ${{\rm{\Omega }}}_{\nu }={m}_{\nu }/93.14/{h}^{2}$. For initial conditions we use ${A}_{s}=2.442\times {10}^{-9}$ and n_s = 0.969. In all cases we use 256³ CDM particles and simulate only a single massive neutrino species. In CUBEP³M the CνB is resolved by 256³ N-body particles, whereas in CUBE they are resolved by three fluids with a grid resolution of 256³. The only exception is for the m_ν = 0.05 eV case where we used 128³ grid cells for the fastest shell. In this simulation we also decreased the Courant–Friedrichs–Lewy (CFL; Courant et al. 1928) condition from 0.7 to 0.5 as we noticed some slight oscillations on small scales in the output. In all simulations for both codes we use z_i = 100 and z_ν = 10. We show slices of CDM and CνB density fields in Figure 1, which demonstrate the lack of Poisson noise in the hydrodynamic simulations well.

We show the ratio of the CνB power spectra to the CDM one in the top panel of Figure 2. Solid lines utilize the collisionless hydrodynamics approach, dashed lines are the particles, dotted lines are linear theory results computed using Code for Anisotropies in the Microwave Background (CAMB; Lewis et al. 2000), and bands are the high-resolution simulations from Cataneo et al. (2020).³ For comparison, we have also included a result from the publicly available Cosmological Massive Neutrino Simulations (MassiveNuS) suite (Liu et al. 2018), which is shown in blue. This simulation included three 0.2 eV neutrinos using the linear response method of Ali-Haïmoud & Bird (2013). This method is designed to use virtually zero memory and processing power on neutrino evolution, which is sufficient to obtain the correct matter power suppression but not the neutrino power spectra directly (unless using the hybrid method described in Bird et al. 2018). It may also miss physical decorrelations between the neutrinos and CDM, like that suggested by Zhu & Castorina (2020). For the lighter neutrino species the hydrodynamic simulations provide significantly better resolution due to lack of Poisson noise, whereas for the heavier ones it is merely comparable (although a lack of power seems preferable to an artificial enhancement). In the bottom panel we show the correlation coefficient between the two codes, defined as the ratio of the cross-power spectrum to the square root of the geometric mean of the autopower spectra. In general, this would go to zero when Poisson noise kicks in. To cancel shot noise in the power spectrum, we divide the neutrino particles randomly into two groups and compute the cross power as described in Inman et al. (2015). For ${m}_{\nu }\lt 0.2$ eV this is very noisy due to the low particle number. However, for the heavier neutrino species, it does do quite well and we find that the hydrodynamics approach is well correlated to at least $k\sim 1$ ${h}^{-1}$Mpc.

Zoom In Zoom Out Reset image size

The modulation of the matter power spectrum is shown in Figure 3. CUBE yields the same result as CUBEP³M, matching linear theory on large scales but with an enhanced dip on nonlinear ones. On nonlinear scales CUBE and CUBEP³M very slightly disagree. We suspect that this is due to the difference in force calculation between the codes, which results in slightly different power spectra.

Zoom In Zoom Out Reset image size

In an ideal simulation method, neutrinos would increase memory and computation time by ${ \mathcal O }({f}_{\nu })$. Our method can satisfy the first criterion easily. As an example, we consider the hypothetical requirements to run a TianNu scale simulation (Emberson et al. 2017) with the CUBE code. TianNu evolved nearly three trillion particles—mostly neutrinos—in a cubic volume of width 1200 ${h}^{-1}$Mpc. Due to Poisson noise, TianNu only resolved neutrino perturbations to $k\sim 1$ h Mpc⁻¹, which generically requires a grid with ∼1024 cells per dimension to resolve. With our new hydrodynamical method, this would require 3 × 5 × 1024³ floating point numbers. This may be compared to the number required in TianNu, 6 × 13824³ floating point numbers, and a savings of approximately ∼10³. The memory usage could be further decreased if faster shells are allowed to be less resolved. CUBE also has a significant CDM memory compression scheme that allows CDM memory reduction to as low as 6 bytes per particle (from 28 in CUBEP³M). In total the TianNu simulation could be run utilizing around 40 times less memory to store particles.

In our current implementation the neutrino computations are quite fast compared to the CDM ones; however, the total computation time is extended well beyond an ${ \mathcal O }({f}_{\nu })$ increase due to the short timestep required by the neutrinos. Nonetheless, we are optimistic that this can be resolved by individual timestepping, where neutrinos take many more short timesteps than CDM. To obtain an idea of how much improvement can be obtained, we extract various timing quantities at each timestep: the time taken in the CDM subroutines (t_c—note that this includes the force calculation, which has some neutrino operations as well), the time taken in the neutrino hydrodynamic routines (t_ν), the required CDM timestep (${\rm{\Delta }}{t}_{c}$), and the required neutrino timestep (${\rm{\Delta }}{t}_{\nu }$) utilizing CFL = 0.5. The number of neutrino steps that can be taken per CDM one is then given by N_ν = Δt_c/Δt_ν.

In Figure 4 we show the ratio of the time necessary for neutrino calculations, ${N}_{\nu }{t}_{\nu }$, compared to the time needed for CDM ones, t_c, for the most challenging 0.05 eV case. At high redshifts, the neutrino computations are substantial; however, at lower redshifts, they quickly become subdominant. Bird et al. (2018) demonstrated that neutrinos with velocities $\gtrsim 100$ km s⁻¹ can be treated with linear response down to $z\sim 1$ and so we expect that we can always set z_ν such that neutrinos are subdominant. Indeed when the subgrid pair-wise force, which requires a much smaller timestep as well, is also used the neutrinos may even become a negligible fraction. We ran a single CDM-only simulation using CUBE with the PP force enabled with an extended range of one fine cell. The dashed–dotted line in Figure 4 shows the ratio of the timestep required for the PP force compared to other non-PP timesteps, giving an estimate of the scaling when subgrid forces are included. CUBE is currently undergoing optimizations for upcoming large-scale simulations, and we hope to include individual timestepping in a subsequent update.

Zoom In Zoom Out Reset image size

We have solved the Vlasov equation for the CνB using a novel technique based on decomposing the homogeneous phase space into velocity shells, which are then evolved via a closed set of Boltzmann moment equations. We have found this approach to accurately model the neutrino perturbations well into the nonlinear regime. This method does not have Poisson noise, which plagues the more common particle based methods. Because of the hydrodynamical method's accuracy and low computational requirements, it is ideally suited for studying novel observables that depend on neutrino perturbations.

To model the CνB we have been interested in subsonic shells (V ≪ u). A natural extension is to look for approximate closure schemes in the supersonic regime, which could provide a useful model for warm dark matter. While the equations we employed are consistent with warm dark matter linear evolution,⁴ it would be quite surprising if they remained valid in highly nonlinear regimes. Nonetheless other closure relationships, perhaps based on higher orders of perturbation theory, may exist and be more valid in nonlinear regimes. Even an approximate scheme may be beneficial, as warm dark matter is known to exhibit artificial fragmentation when simulated using N-body particles (Wang & White 2007), which could potentially be ameliorated through the use of hydrodynamic equations.

We are extremely grateful to Andrew MacFadyen for providing notes on implementing hydrodynamics. We thank Yacine Ali-Haïmoud for valuable discussions, J.D. Emberson for assistance with the power spectra of Cataneo et al. (2020), and Jia Liu for help with the MassiveNuS (Liu et al. 2018) power spectra. We thank the Columbia Lensing group (http://columbialensing.org) for making their suite of simulated maps available and NSF for supporting the creation of those maps through grant AST-1210877 and XSEDE allocation AST-140041. We thank New Mexico State University (USA) and Instituto de Astrofisica de Andalucia CSIC (Spain) for hosting the Skies & Universes site for cosmological simulation products. H.R.Y. acknowledges National Science Foundation of China No. 11903021. This work was supported in part through the NYU IT High Performance Computing resources, services, and staff expertise. This research has made use of NASA's Astrophysics Data System Bibliographic Services.

Software: NumPy v1.13.3 (van der Walt et al. 2011), Matplotlib v2.1.0(Hunter 2007), SciPy v1.0.0rc1 (Jones et al. 2001).

In this section we test to what extent the isothermal and ideal gas approximations can be used to model the shells. For the isothermal gas, we assume that P = ρu² at all times, with the sound speed of c_s = u. We show the results in Figure A1. In the top panel the solid colored lines show the isothermal approximation, whereas the solid black lines show the collisionless hydrodynamics result. The bottom panel shows the difference in the correlation coefficient between the isothermal gas and the collisionless hydrodynamics, i.e., negative values indicate collisionless hydrodynamics is better correlated. We find that the isothermal approximation appears to overpredict neutrino perturbations but remains well correlated with the particle field.

Zoom In Zoom Out Reset image size

For the ideal gas approximation, we employ the ideal gas energy equation:

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}E+\displaystyle \frac{\partial }{\partial {x}^{k}}\left[(E+P){V}^{k}\right]=\rho {V}^{k}{g}^{k}.\end{eqnarray} \tag{ A1 }$$

Because the flux $(E+P){V}^{k}$ is simply proportional to V^k upon linearization—the same as in collisionless hydrodynamics—we expect that we can reproduce linear evolution by simply changing the initial energy to

$$\begin{eqnarray}&&{E}_{u}={\rho }_{u}\left[\displaystyle \frac{1}{2}{V}_{u}^{2}+\displaystyle \frac{{u}^{2}}{\gamma (\gamma -1)}\right]\end{eqnarray} \tag{ A2 }$$

and the sound speed is the standard $\sqrt{\gamma P/\rho }$. We show the result in Figure A2, which shows that the ideal gas approximation slightly underpredicts neutrino perturbations. Again, this does not result in a substantially worsened correlation coefficient.

Zoom In Zoom Out Reset image size

We conclude that both approximations are sufficient to model the CνB, although perhaps slightly less accurate than collisionless hydrodynamics. As both isothermal and ideal gases are commonly included in cosmological codes, it is therefore straightforward for these codes to include Poisson noise-free neutrino perturbations provided they can simulate multiple fluids simultaneously.