The cosmic neutrino background as a collection of fluids in large-scale structure simulations

Consider a universe containing only a relic population of neutrinos of mass $m_{\nu}$, subdivided into an infinite number of fluids according to some property $\vec{\tau}$. At position $\vec{x}$ and time $s$, each of these neutrino fluids is characterised by a number density $n(s,\vec{x},\vec{\tau})$ and a 4-momentum $P^{\mu}(s,\vec{x},\vec{\tau})$, whose spatial means at a given time $s$ are $\bar{n}(s,\vec{\tau})$ and $\bar{P}^{\mu}(s,\vec{\tau})$ respectively. The mean number density of one such fluid evolves with the scale factor as $\bar{n}(s,\vec{\tau})\propto a^{-3}(s)$ as usual. On the other hand, the spatial part of the mean 4-momentum with upper Lorentz indices scales as $\bar{P}^{i}(s,\vec{\tau})\propto a^{-2}(s)$, while its lower index version $\bar{P}_{i}(s,\vec{\tau})$ is constant in time. The latter leads us to identify $\tau_{i}:=\bar{P}_{i}(s,\vec{\tau})$, and it follows that the mean coordinate 3-velocity of the fluid is $\vec{v}(s,\vec{\tau})=-\vec{\tau}/\xi(s,\tau)$, where $\tau:=|\vec{\tau}|$ and $\xi(s,\tau):=-\sqrt{m_{\nu}^{2}a^{2}(s)+\tau^{2}}$.

At early times when the spatial inhomogeneities are small, we expect $P_{i}(s,\vec{x},\vec{\tau})\rightarrow\tau_{i}$. Then, subdividing the neutrino population into fluids by $\vec{\tau}$ is equivalent to labelling these fluids by their initial comoving 3-momenta. Each fluid thus characterised is unique, and its subsequent evolution under the influence of spatial inhomogeneities can be treated with standard fluid perturbation theory. Reference [36] has derived such a linear perturbation theory for a neutrino fluid characterised by an initial momentum $\vec{\tau}$, which forms the starting point for our analysis. In Fourier space, the equations of motion for a fluid’s number density contrast $\delta(s,\vec{k},\vec{\tau}):=n(s,\vec{k},\vec{\tau})/\bar{n}(s,\vec{\tau})-1$ and the momentum divergence $\theta_{P}(s,\vec{k},\vec{\tau}):=ik^{i}P_{i}(s,\vec{k},\vec{\tau})$ are respectively [36]

	$\displaystyle\delta^{\prime}$	$\displaystyle=$	$\displaystyle i\frac{\vec{k}\cdot\vec{\tau}}{{\mathcal{H}}\xi}\delta+\left(1-\frac{(\vec{k}\cdot\vec{\tau})^{2}}{k^{2}\xi^{2}}\right)\theta_{P}+3\Psi^{\prime}-i\frac{\vec{k}\cdot\vec{\tau}}{{\mathcal{H}}\xi}\left[\left(1+\frac{\tau^{2}}{\xi^{2}}\right)\Psi-\Phi\right],$		(2.4)
	$\displaystyle\theta_{P}^{\prime}$	$\displaystyle=$	$\displaystyle i\frac{\vec{k}\cdot\vec{\tau}}{{\mathcal{H}}\xi}\theta_{P}-\frac{\xi k^{2}}{{\mathcal{H}}}\Phi-\frac{\tau^{2}k^{2}}{\xi{\mathcal{H}}}\Psi.$		(2.5)

The key difference between these equations and standard linearised fluid equations for, e.g., CDM, is that the initial fluid momentum $\vec{\tau}$ now also enters the picture. Thus, this multi-fluid treatment of neutrinos is inherently Lagrangian in momentum space.

Two observations about equations (2.4) and (2.5) simplify their analysis. Firstly, the direction of $\vec{\tau}$ appears only through its inner product with the Fourier vector, $\mu:=\vec{k}\cdot\vec{\tau}/(k\tau)$. This together with the fact that the initial conditions are also independent of the absolute direction of $\vec{\tau}$ means that we need not consider the full set of allowed $\vec{\tau}$ when representing the neutrino population as a set of discrete fluids. It suffices to deal only with $\tau$ and $\mu$, and modelling relic neutrinos can proceed first by partitioning the population by $\tau$ into $N_{\tau}$ fluids, followed by approximating each these $N_{\tau}$ fluids as having $N_{\mu}$ angular degrees of freedom, in the form of either a discrete set of $\mu$ as in reference [36], or a finite-dimensional subspace of the set of all functions of $\mu$.

Secondly, different fluids interact only through the gravitational potentials $\Phi(s,k)$ and $\Psi(s,k)$. For an infinite number of fluids, each with density contrast $\delta(s,k,\tau,\mu)$ and velocity perturbations $\theta_{P}(s,k,\tau,\mu)$, the potentials in the conformal Newtonian gauge are [37]

	$\displaystyle k^{2}\Psi+3{\mathcal{H}}^{2}(\Psi^{\prime}+\Phi)$	$\displaystyle=$	$\displaystyle-4\pi Ga^{2}\delta\rho,$		(2.6)
	$\displaystyle k^{2}(\Psi^{\prime}+\Phi)$	$\displaystyle=$	$\displaystyle-4\pi Ga^{2}\delta U,$		(2.7)
	$\displaystyle\delta\rho:=-\frac{1}{2}\int_{-1}^{1}{\rm d}\mu\,\delta T^{0}_{0}$	$\displaystyle=$	$\displaystyle\bar{\rho}_{\mathrm{crit}}(\tau)\int_{0}^{\infty}{\rm d}\tau\,\Omega_{\nu}(s,\tau)\int_{-1}^{1}\frac{{\rm d}\mu}{2}\left(\delta-\frac{\tau\mu}{k\xi^{2}}\theta_{P}+\frac{\tau^{2}}{\xi^{2}}\Psi\right),$		(2.8)
	$\displaystyle\delta U:=-\frac{1}{2}\int_{-1}^{1}{\rm d}\mu\frac{ik_{j}}{{\mathcal{H}}}\delta T^{0}_{j}$	$\displaystyle=$	$\displaystyle\bar{\rho}_{\mathrm{crit}}(\tau)\int_{0}^{\infty}{\rm d}\tau\,\Omega_{\nu}(s,\tau)\int_{-1}^{1}\frac{{\rm d}\mu}{2}\left(\frac{\theta_{P}+k\tau\mu\delta}{{\mathcal{H}}\xi}\right),$		(2.9)

where $\Omega_{\nu}(s,\tau)=\Omega_{\nu 0}(\tau)\,{\mathcal{E}}(s,\tau){\mathcal{H}}_{0}^{2}/{\mathcal{H}}^{2}$, with ${\mathcal{E}}(s,\tau)=|\xi|/(m_{\nu}a^{2})$. At a given $\tau$, these potentials depend on the quantities $\delta_{\ell=0}(s,k,\tau)=\frac{1}{2}\int_{-1}^{1}{\rm d}\mu\,\delta$, $\delta_{\ell=1}(s,k,\tau)=\frac{1}{2}\int_{-1}^{1}{\rm d}\mu\,\mu\delta$, $\theta_{P,\ell=0}(s,k,\tau)=\frac{1}{2}\int_{-1}^{1}{\rm d}\mu\,\theta_{P}$, and $\theta_{P,\ell=1}(s,k,\tau)=\frac{1}{2}\int_{-1}^{1}{\rm d}\mu\,\mu\theta_{P}$, which are the Legendre monopoles and dipoles of the fluid perturbations. This motivates us to represent the $N_{\mu}$ angular degrees of freedom using a Legendre polynomial expansion.

We describe first in the following how we sample $N_{\tau}$ fluids following the method of reference [36]. Our Legendre representation of the fluids’ angular degrees of freedom differs from the approach of [36], and we defer its discussion to section 3.1.

We sample the initial neutrino fluid momenta $\tau$ from an ultrarelativistic Fermi–Dirac distribution. Defining $x:=\tau/(aT_{\nu})$, where $T_{\nu}(s)=T_{\nu,0}/a(s)$, with $T_{\nu,0}=1.95$ K, is the neutrino temperature, the (time-indepemdent) mean comoving number density of the full neutrino population is given up to a constant multiplicative factor by

	$\displaystyle\bar{n}_{\mathrm{comoving}}$	$\displaystyle\propto\,T^{3}_{\nu,0}\int_{0}^{\infty}{\rm d}x\,\frac{x^{2}}{1+\exp(x)}$		(2.10)
		$\displaystyle=\,\frac{3\zeta(3)}{2}T^{3}_{\nu,0},$

where $\zeta(3)$ is a Riemann zeta function. The cumulative probability ${\mathfrak{P}}(\tau)$ of finding a neutrino with momentum less than $\tau$ is therefore

$${\mathfrak{P}}(\tau)=\frac{2}{3\zeta(3)}\int_{0}^{\tau/T_{\nu,0}}{\rm d}x\,\frac{x^{2}}{1+\exp(x)}.$$

(2.11)

This is by definition a monotonically increasing function. It may therefore be inverted, thereby defining $\tau(\mathfrak{P})$ for $0\leq\mathfrak{P}\leq 1$.

We choose to divide the neutrino population into $N_{\tau}$ bins of equal number densities, implying that all neutrino density fractions $\Omega_{\alpha}(s)$ are equal in the non-relativistic limit. Then, for each neutrino fluid $\alpha=0,\ldots,N_{\tau}-1$ in our collection of $N_{\tau}$ fluids, a convenient choice of $\tau_{\alpha}$ is $\tau_{\alpha}=\tau({\mathfrak{P}}_{\alpha})$ for ${\mathfrak{P}}_{\alpha}=(\alpha+\frac{1}{2})/N_{\tau}$. With this choice, the parameter $N_{\tau}$ determines the momentum resolution within the neutrino population. Thus, although the multi-fluid method allows for an arbitrary binning, which may be tailored for specific applications, we do not investigate this possibility further.

We begin by replacing in equation (2.5) the momentum divergence by the velocity divergence for each neutrino fluid $\alpha$, $\theta_{\alpha}:=\theta_{P,\alpha}/({\mathcal{H}}\xi_{\alpha})$. The resulting equation of motion is

$$\theta_{\alpha}^{\prime}=-\left(1-v_{\alpha}^{2}+\frac{{\mathcal{H}}^{\prime}}{{\mathcal{H}}}\right)\theta_{\alpha}-\frac{ik\mu v_{\alpha}}{{\mathcal{H}}}\theta_{\alpha}-\frac{k^{2}}{{\mathcal{H}}^{2}}(\Phi+v_{\alpha}^{2}\Psi).$$

(3.1)

In the non-relativistic limit we may discard terms of order $v_{\alpha}^{2}$, so that

$$\theta_{\alpha}^{\prime}=-\left(1+\frac{{\mathcal{H}}^{\prime}}{{\mathcal{H}}}\right)\theta_{\alpha}-\frac{ik\mu v_{\alpha}}{{\mathcal{H}}}\theta_{\alpha}-\frac{k^{2}}{{\mathcal{H}}^{2}}\Phi.$$

(3.2)

Similarly, recast in terms of $\theta_{\alpha}$ equation (2.4) becomes

$$\delta_{\alpha}^{\prime}=-\frac{ik\mu v_{\alpha}}{{\mathcal{H}}}\delta_{\alpha}+(1-\mu^{2}v_{\alpha}^{2})\theta_{\alpha}+3\Psi^{\prime}+\frac{ik\mu v_{\alpha}}{{\mathcal{H}}}\left[(1+v_{\alpha}^{2})\Psi-\Phi\right],$$

(3.3)

which reduces further in the subhorizon non-relativistic limit to

$$\delta_{\alpha}^{\prime}=-\frac{ik\mu v_{\alpha}}{{\mathcal{H}}}\delta_{\alpha}+\theta_{\alpha}$$

(3.4)

if we neglect all terms of order $v_{\alpha}^{2}$ as well as gravitational potentials not multiplied by $k^{2}/{\mathcal{H}}^{2}$.

Departing from reference [36], we now expand $\delta_{\alpha}(k,\mu)$ and $\theta_{\alpha}(k,\mu)$ in Legendre polynomials ${\mathcal{P}}_{\ell}(\mu)$. For an arbitrary function $X(\mu)$, we define the Legendre coefficients $X_{\ell}$ via

	$\displaystyle X(\mu)$	$\displaystyle=\,\sum_{\ell=0}^{\infty}(-i)^{\ell}{\mathcal{P}}_{\ell}(\mu)X_{\ell},$		(3.5)
	$\displaystyle X_{\ell}$	$\displaystyle=\,\frac{i^{\ell}}{2}(2\ell+1)\int_{-1}^{1}d\mu{\mathcal{P}}_{\ell}(\mu)X(\mu).$

This particular choice of definition follows from the structure of our equations of motion (3.2) and (3.4), which couple terms odd and even in $\mu$ only through factors proportional to $i\mu$; the factors of $i$ in the definition (3.5) therefore ensure that the Legendre coefficients $\delta_{\alpha,\ell}(k)$ and $\theta_{\alpha,\ell}(k)$ initialised to real values will remain real under time evolution. The corresponding Legendre-expanded equations of motion are

	$\displaystyle\delta_{\alpha,\ell}^{\prime}$	$\displaystyle=$	$\displaystyle\theta_{\alpha,\ell}+\frac{kv_{\alpha}}{{\mathcal{H}}}\left(\frac{\ell}{2\ell-1}\delta_{\alpha,\ell-1}-\frac{\ell+1}{2\ell+3}\delta_{\alpha,\ell+1}\right),$		(3.6)
	$\displaystyle\theta_{\alpha,\ell}^{\prime}$	$\displaystyle=$	$\displaystyle-\delta^{\mathrm{(K)}}_{\ell 0}\frac{k^{2}}{{\mathcal{H}}^{2}}\Phi-\left(1+\frac{{\mathcal{H}}^{\prime}}{{\mathcal{H}}}\right)\theta_{\alpha,\ell}+\frac{kv_{\alpha}}{{\mathcal{H}}}\left(\frac{\ell}{2\ell-1}\theta_{\alpha,\ell-1}-\frac{\ell+1}{2\ell+3}\theta_{\alpha,\ell+1}\right),$		(3.7)

where $\delta^{\mathrm{(K)}}_{ij}$ is the Kronecker delta function.

This system of equations couples each $\ell\geq 0$ to $\ell+1$, resulting in an infinite system of equations that needs to be truncated in order to limit the decompositions of $\delta_{\alpha}(k,\mu)$ and $\theta_{\alpha}(k,\mu)$ to a finite number $N_{\mu}$ of Legendre polynomials. Reference [39] warns against truncating this system at a finite $\ell_{\mathrm{max}}=N_{\mu}-1$ by merely setting all higher coefficients to zero; doing so would amount to imposing a reflecting boundary condition at $\ell_{\mathrm{max}}$, thereby sending power that ought to flow towards infinite $\ell$ back towards $\ell=0$. Thus, we must approximate $\delta_{\alpha,\ell}$ and $\theta_{\alpha,\ell}$ at $\ell=\ell_{\mathrm{max}}+1$.

In order to do this, note that the terms multiplying $kv_{\alpha}/{\mathcal{H}}$ in equations (3.6) and (3.7) in the limit of large $\ell$ approach finite-difference approximations to the first derivatives of $\delta_{\alpha,\ell}$ and $\theta_{\alpha,\ell}$ with respect to $\ell$. We make this explicit by defining $f_{\ell}(x)=x\delta_{\alpha,x\ell}/(x\ell+1/2)$. Then, the centered second-order finite-difference approximation to the first derivative is $\Delta^{\mathrm{(C)}}_{2,\epsilon}f_{\ell}(x)=[f_{\ell}(x+\epsilon)-f_{\ell}(x-\epsilon)]/(2\epsilon)$, with the actual derivative ${\rm d}f_{\ell}(x)/{\rm d}x=\Delta^{\mathrm{(C)}}_{2,\epsilon}f_{\ell}(x)+{\mathcal{O}}(\epsilon^{2})$. For $x=1$ and $\epsilon=1/\ell$, we find

$$\Delta^{\mathrm{(C)}}_{2,\frac{1}{\ell}}f_{\ell}(1)=\frac{\ell+1}{2\ell+3}\delta_{\alpha,\ell+1}-\frac{\ell-1}{2\ell-1}\delta_{\alpha,\ell-1},$$

(3.8)

which can be easily rearranged to give

$$\frac{\ell}{2\ell-1}\delta_{\alpha,\ell-1}-\frac{\ell+1}{2\ell+3}\delta_{\alpha,\ell+1}=\frac{1}{2\ell-1}\delta_{\alpha,\ell-1}-\Delta^{\mathrm{(C)}}_{2,\frac{1}{\ell}}f_{\ell}(1),$$

(3.9)

i.e., precisely the term in parenthesis in equation (3.6) now written in terms of $\Delta^{\mathrm{(C)}}_{2,1/\ell}f_{\ell}(1)$.

Next, we replace the centered derivative $\Delta^{\mathrm{(C)}}_{2,\epsilon}f_{\ell}(x)$ at $\ell=N_{\mu}-1$ by the backwards derivative $\Delta^{\mathrm{(B)}}_{2,\epsilon}f_{\ell}(x)=[f_{\ell}(x-2\epsilon])-4f_{\ell}(x-\epsilon)+3f_{\ell}(x)]/(2\epsilon)$, which also approximates ${\rm d}f_{\ell}(x)/{\rm d}x$ to ${\mathcal{O}}(\epsilon^{2})\approx{\mathcal{O}}(N_{\mu}^{-2})$. The key is that $\Delta^{\mathrm{(B)}}_{2,1/\ell_{\mathrm{max}}}f_{\ell}(1)$ depends only upon $\delta_{\alpha,\ell}$ for $\ell\leq\ell_{\mathrm{max}}$. This replacement amounts to approximating

$$\frac{\ell_{\mathrm{max}}+1}{2\ell_{\mathrm{max}}+3}\delta_{\alpha,\ell_{\mathrm{max}}+1}\approx\frac{3\ell_{\mathrm{max}}}{2\ell_{\mathrm{max}}+1}\delta_{\alpha,\ell_{\mathrm{max}}}-\frac{3(\ell_{\mathrm{max}}-1)}{2\ell_{\mathrm{max}}-1}\delta_{\alpha,\ell_{\mathrm{max}}-1}+\frac{(\ell_{\mathrm{max}}-2)}{2\ell_{\mathrm{max}}-3}\delta_{\alpha,\ell_{\mathrm{max}}-2}.$$

(3.10)

Similar arguments apply to the truncation of equation (3.7). In practice, we find $N_{\mu}\approx 10$ to be sufficient for percent-level accuracy in $\delta_{\alpha,\ell_{\mathrm{max}}}$ and $\theta_{\alpha,\ell_{\mathrm{max}}}$, which themselves should be subdominant at late times. This conclusion is consistent with that of reference [36], which found that high-accuracy neutrino calculations do not require a high resolution of the $\mu$-dependence of the perturbations.

Finally, in the subhorizon limit ${\mathcal{H}}\ll k$, the Newtonian-gauge potentials are equal, i.e., $\Psi=\Phi$, and Poisson’s equation relates them to the density perturbations of the various constituents of the universe:

$$k^{2}\Phi(s,k)=-\frac{3}{2}{\mathcal{H}}^{2}(s)\left(\Omega_{\rm cb}(s)\delta_{\rm cb}(s,k)+\sum_{\alpha=0}^{N_{\tau}-1}\Omega_{\alpha}(s)\delta_{\alpha,\ell=0}(s,k)\right).$$

(3.11)

Here, $\Omega_{\rm cb}(s)$ and $\delta_{\rm cb}(s,k)$ denote respectively the time-dependent energy density fraction and density contrast of the combined CDM+baryon fluid, $\delta_{\alpha,\ell=0}$ is the monopole of the $\alpha$th neutrino fluid density contrast, and the summation is over all $N_{\tau}$ neutrino fluids weighted by $\Omega_{\alpha}(s)$, which we have chosen to be equal for non-relativistic neutrinos.

Then, to summarise, our multi-fluid neutrino perturbation theory uses the subhorizon non-relativistic equations of motion (3.6) and (3.7) truncated as per equation (3.10) at $\ell_{\mathrm{max}}=N_{\mu}-1$, and with a gravitational potential given by Poisson’s equation (3.11).

Closed-form growing-mode solutions to equations (3.6), (3.7), (3.11), (3.12) and (3.13) in a generic $\nu\Lambda$CDM universe are not known. However, since the growing mode is an attractor solution, an approximate solution suffices for the purpose of setting initial conditions. We describe our approximation below. In all cases, we initialise at the scale factor $a_{\mathrm{in}}=10^{-3}$.

Consider first a universe containing only cold matter and radiation. The growing-mode solution for linear CDM+bayron perturbations in such a universe is well known and given by

	$\displaystyle\delta_{\rm cb}(s)$	$\displaystyle=$	$\displaystyle a(s)+\frac{2}{3}a_{\mathrm{eq}},$		(3.14)
	$\displaystyle\theta_{\rm cb}(s)$	$\displaystyle=$	$\displaystyle\delta^{\prime}(s)=a,$		(3.15)

where $a_{\mathrm{eq}}$ is the scale factor at matter–radiation equality. To adapt these expressions into approximate growing-mode initial conditions for $\delta_{\rm cb}$ and $\theta_{\rm cb}$ in a realistic universe that also contains massive neutrinos, we note first of all that at $a_{\rm in}=10^{-3}$ massive neutrinos are typically neither radiation nor cold matter but are transitioning the former to the latter. Then, to approximate this transition effect, we set

$$a_{\mathrm{eq}}=\frac{\Omega_{\gamma 0}+\Omega_{\nu\mathrm{(massless),}0}}{\Omega_{\mathrm{cb}0}},$$

(3.16)

where $\Omega_{\nu\mathrm{(massless),}0}$ is the $z=0$ density fraction corresponding to massless neutrinos with temperature $T_{\nu,0}=1.95$ K. Thus the radiation component counts only those constituents that are truly massless, while only CDM+baryons contribute to the cold matter. In other words, the energy density in massive neutrinos does not feature in the evaluation of $a_{\mathrm{eq}}$ for the purpose of setting initial conditions via equations (3.14) and (3.15).

Meanwhile in the neutrino sector, we note that given a CDM+baryon density contrast $\delta_{\mathrm{cb}}$, the total neutrino monopole density contrast, $\delta_{\nu,\ell=0}:\sum_{\alpha=0}^{N_{\tau}-1}\Omega_{\alpha}\delta_{\alpha,\ell=0}/\sum_{\alpha=0}^{N_{\tau}-1}\Omega_{\alpha}$, can be approximated by [49, 50]

$$\frac{\delta_{\nu,0}(s,k)}{\delta_{\mathrm{cb}}(s,k)}=\frac{k_{\mathrm{FS}}^{2}(1-f_{\nu})}{(k+k_{\mathrm{FS}})^{2}-f_{\nu}k_{\mathrm{FS}}^{2}},$$

(3.17)

where $f_{\nu}:=\Omega_{\rm 0}/\Omega_{\nu 0}$ is the neutrino mass fraction, and

$$k_{\mathrm{FS}}(s)=\frac{{\mathcal{H}}}{c_{\nu}}\sqrt{\frac{3}{2}\Omega_{\mathrm{m}}(s)}=a{\mathcal{H}}\frac{m_{\nu}}{T_{\nu,0}}\sqrt{\frac{\ln(2)}{\zeta(3)}\Omega_{\mathrm{m}}(s)}$$

(3.18)

is the time-dependent neutrino free-streaming wave number, with $c_{\nu}$ denoting a characteristic thermal speed of the relic neutrino population. Inspired by these approximations, we initialise the monopole ($\ell=0$) of the neutrino fluid $\alpha$ to

	$\displaystyle\delta_{\alpha,0}(s,k)$	$\displaystyle=$	$\displaystyle\frac{k_{\mathrm{FS},\alpha}^{2}(1-f_{\nu})}{(k+k_{\mathrm{FS},\alpha})^{2}-f_{\nu}k_{\mathrm{FS},\alpha}^{2}}\,\delta_{\mathrm{cb}}(s,k),$		(3.19)
	$\displaystyle\theta_{\alpha,0}(s,k)$	$\displaystyle=$	$\displaystyle\delta_{\alpha,0}^{\prime}(s,k),$		(3.20)

where the $\alpha$-dependent free-streaming wave number is now given by

$$k_{\mathrm{FS},\alpha}(s)=\frac{{\mathcal{H}}}{v^{\rm NR}_{\alpha}}\sqrt{\frac{3}{2}\Omega_{\mathrm{m}}(s)}=\frac{m_{\nu}a{\mathcal{H}}}{\tau_{\alpha}}\sqrt{\frac{3}{2}\Omega_{\mathrm{m}}(s)}\,,$$

(3.21)

with the characteristic thermal speed $c_{\nu}$ now replaced with the fluid’s non-relativistic velocity $v_{\alpha}^{\rm NR}:=\tau_{\alpha}/(m_{\nu}a)$. All $\ell>0$ terms of the neutrino density contrast and velocity divergence are initialised to zero at $a_{\mathrm{in}}$.

In summary, we have demonstrated the following about the multi-fluid linear treatment of neutrinos presented above.

• •

  Late-time ($a\geq 0.4$) perturbations converge as the numbers of fluids $N_{\tau}$ and angular modes $N_{\mu}$ are increased, with the CDM+baryon density and velocity perturbations consistent at the $0.1\%$ level for $N_{\tau}=N_{\mu}=10$.

• •

  Neutrino density perturbations agree at the percent level with predictions of the state-of-the-art integral linear response method [27] up to $k\approx 0.1~h/$Mpc for $N_{\tau}=N_{\mu}=10$, and up to $k\approx 2~h/$Mpc for $N_{\tau}=1000$, $N_{\mu}=100$.

• •

  The CDM+baryon growth factor normalised at $z=0$ agrees with the widely-used class Boltzmann code at the $0.1\%$ level for $a\geq 0.4$. Neutrino density perturbations agree to $\approx 2\%$ over that same range up to $k=0.1~h/$Mpc for $N_{\tau}=N_{\mu}=10$, and up to $k=1~h/$Mpc for $N_{\tau}=1000$, $N_{\mu}=100$.

Thus, our perturbation theory is accurate to $0.1\%$ for CDM+baryons at all wavenumbers $k$, and to $\approx 1\%$ for the neutrinos below their free-streaming wave number $k_{\rm FS}$, even for modest numbers of fluids and angular modes, $N_{\tau}=N_{\mu}=10$.

In order to model the time-dependent effects of neutrino clustering on non-linear CDM+ baryon perturbations, we use the full version of Time-RG, rather than the one-loop approximation applied in, e.g., the data analysis of [4]. Full Time-RG and similar perturbation theories are known to suffer from a numerical instability at very small scales [70, 71, 72], which drives the velocity auto-power spectrum to negative values. References [70, 71, 72] trace this instability to the neglect of the velocity field’s vorticity in the CDM+baryon equations of motion, thereby predicting that the instability should arise at length scales much smaller than the velocity dispersion length $\sigma_{v}(\eta):=[\int dkP_{\mathrm{cb}}(k)/(6\pi^{2})]^{1/2}$, or, equivalently, $k\gg k_{v}:=2\pi/\sigma_{v}$. We address this instability by smoothing the CDM+baryon perturbation inputs to the non-linear mode-coupling integrals of the Time-RG equations of motion, via the replacement $P_{\rm cb}(k)\rightarrow P_{\rm cb}(k)/(1+k^{2}/k_{v}^{2})$. This smoothing stabilises Time-RG at $k\gg 1~h/$Mpc.

Observe also that Time-RG has been developed for Newtonian gravity [55, 56]. Thus,to ensure that our computations agree with the outputs of such relativistic linear Boltzmann codes as class and camb on large scales, we use a “scale back” procedure when performing our non-linear perturbative computation: Firstly, we switch off non-linear corrections altogether and run our Time-RG code to $z=0$ using the initial conditions of section 3.2. This exercise yields $\delta_{\mathrm{cb}}(k)$ and $\delta_{\alpha,\ell}(k)$ at $z=0$. Secondly, given the $z=0$ total matter power spectrum $P_{\mathrm{m}}(k)$ from camb or class, we define the $k$-dependent normalisation ${\mathcal{N}}(k)$ via

$$\frac{{\mathcal{N}}(k)}{\Omega_{\mathrm{m}0}}\left[\Omega_{\mathrm{cb}0}\delta_{\mathrm{cb}}(k)+\frac{\Omega_{\nu 0}}{N_{\tau}}\sum_{\alpha=0}^{N_{\tau}-1}\delta_{\alpha,0}(k)\right]=\sqrt{P_{\mathrm{m}}(k)}\,.$$

(5.1)

Finally, we rerun the perturbation theory from the initial scale factor $a_{\mathrm{in}}$, but this time with (i) the non-linear corrections restored, and (ii) the growing-mode initial perturbations of section 3.2 all multiplied by ${\mathcal{N}}(k)$. This procedure ensures that non-linear mode-couplings are computed using the properly-normalised linear power spectrum.

Figure 5 compares the CDM+baryon power spectrum outputs of our Time-RG implementation and of the CosmicEmu $N$-body simulation emulator [73] for the reference cosmology of table 1. Though perturbation theory is not a precision tool at $k\gg 0.1~h/$Mpc, it agrees with the CosmicEmu output up to a multiplicative factor of two or three over the entire $k$-range covered by the emulator, evident in the left panel of figure 5.

The right panel of figure 5, on the other hand, shows the logarithmic derivative of the non-linear CDM+baryon power spectrum with respect to the linear power spectrum normalisation $\sigma_{8}^{2}$. This derivative approaches unity at large scales and rises at smaller scales as non-linearity enhances growth, and is a useful measure of the sensitivity of the non-linear power spectrum to small changes to the linear power spectrum normalisation such as represented by free-streaming neutrinos. Here, we see that while the Time-RG prediction of the logarithmic derivative agrees with that of CosmicEmu at $z=2$, at later times the former overeshoots the latter by a factor of two. We therefore conclude on the basis of figure 5 that Time-RG suffices for estimating the impact of non-linear CDM+baryon clustering at the factor-of-two level.

Equipped with the Time-RG theory for non-linear CDM+baryon perturbations, coupled to a multi-fluid linear theory for neutrinos that respond to the CDM+baryon non-linearities, we are now in a position to answer several important questions in an approximate way:

1. 1.

   By how much is the linear clustering of neutrinos enhanced when the CDM+baryon fluid is allowed to cluster non-linearly?

2. 2.

   What portion of the neutrino population clusters strongly enough that linear response is no longer adequate?

3. 3.

   How much does neutrino clustering, as modelled by linear response calculations, enhance the clustering of the CDM+baryon fluid?

These questions are addressed in three subsections below using our Time-RG plus multi-fluid perturbation theory. We shall revisit them later on in section 6 using more accurate calculations of non-linear CDM+baryon clustering from $N$-body simulations