# Neutrino Modelling in Friedmann Equation

I'm trying to model neutrinos in the Friedmann Equation. I've covered the case of the Benchmark Model where we have matter, radiation, curvature, and the cosmological constant, Lambda. I know my coding of the Friedmann equation works because I get the correct plots at different parameters, as you'll see attached below.

Including neutrinos, the Friedmann Equation becomes

$$\begin{eqnarray} H(z)^2 & = & H_0^2 \Big[ (\Omega_c + \Omega_b) ( 1 + z)^3 + \Omega_\gamma ( 1 + z) ^4 \\ & + & \Omega_{DE} ( 1 + z)^{3(1+w)} + \Omega_k ( 1 + z) ^2 + \frac{\rho_{\nu, tot}(z)}{\rho_{crit,0}} \Big]. \end{eqnarray}$$

To solve for the energy density as a function of the scale factor (or redshift), we can solve for the energy density by the following expression for the energy density of a single neutrino species:

$$\rho_\nu (T_\nu) = \frac{g}{(2\pi)^3} \int \frac{\sqrt{p^2 + m^2}}{e^{p/T_\nu} + 1} d^3 p.$$

The energy density critical is $$4870$$ Mev/m$$^3$$. Energy density of the single species can be written as function of scale factor by writing the temperature as a function of scale factor. $$T$$ is simply the expression shown below divided by a: In equation 17, we can write $$d^3p$$ as $$4\pi p^2 dp$$ and $$g=2$$ for a neutrino species. Another thing to notice is that (17) is written in natural units where $$c = h = k = 1$$. I've tried to fix the units and no matter what I do, the density parameter of the neutrino species is always very small (order of $$10^{-9}$$) where it should be between 0.0013 and 0.007 from Ryden, Intro to Cosmology equation (7.54).

I was really hoping someone can help me with the unit conversion from the natural units to the proper units. Everything else I've figured out, I just can't seem to fix the units for equation (17).

Without neutrinos, I get the following plot consisting of various universe models, and they are correct so coding is not the problem. The problem is the unit conversion to proper SI units of (17). Once I get the neutrinos figured out, I want to see how they affect the universe models. Any help is greatly appreciated!

• What is equation (17)? What have you done yourself? I can't see what answer you are expecting. Your equation involves the ratio of the neutrIno density to the critical density so it's unitless and it doesn't matter what units you use so long as they are consistent. What value of $m$ you are using. How have you done the integral? Apr 1 '20 at 9:27
• Sorry, equation 17 is the one for energy density of the single neutrino species; it cut off for some reason. For the values of m, we know that theres a constraint on the sum of neutrino masses, which is roughly 0.06eV/c^2. Just to check, I'm assuming that for a single neutrino species, mass is 0.02eV/c^2. I agree with your statement of the ratio being unitless. However, I'm not sure what the units of (17) are. Since they used natural units, they set c=h=k=1, so units are not consistent in the ratio. I've never encountered natural units so I dont know how to fix (17) such that we have MeV/m^3. Apr 1 '20 at 10:38
• @RobinDhillon Take everything in the same units such as eV. Apr 1 '20 at 18:02

The energy density of a Fermi gas is $$\rho_{\nu}= \int \rho(p)\ dp = \int E(p)F(p)g(p)\ dp$$ $$\rho_{\nu} = \int \left(\sqrt{p^2 c^2 + m^2 c^4}\right)\left(\exp (E/k_BT) + 1\right)^{-1} \left(g_s 4\pi p^2/h^3\right)\ dp$$ in units of energy per unit volume.
Before neutrino decoupling at $$k_B T \sim 1$$ MeV, the neutrinos are ultrarelativistic with $$pc \gg m_{\nu}c^2$$. After decoupling, the shape of the occupation index function $$F(p)$$ does not change - so $$F(p) = \left(pc/k_BT_{\nu} +1\right)^{-1}$$ in subsequent evolution.
Thus $$\rho_{\nu} = \frac{g_s c}{h^3} \int \frac{ \sqrt{p^2 + m_{\nu}^2c^2}} {\exp(pc/k_BT_{\nu}) +1}\ 4\pi p^2\ dp$$
I don't understand where your $$(2\pi)^3$$ comes from, other than to suggest that the unit system is actually $$\hbar = 1$$.