Add massive neutrinos

[jmsull]

Working implementation of massive neutrinos, but some details to sort out before merge.

Summary:

Bolt.jl - add sum of neutrino mass parameters [for now only does a single neutrino with mass given by the sum]

background.jl:

• add non-relativistic neutrino contribution to \Omega_\Lambda
• add background Fermi-Dirac distribution function and lnq derivative
• add integral over comoving momentum q to compute the background massive neutrino density and pressure contributions to H(a) ^this is working well when compared to relativistic and non-relativistic analytic limits, it might be good to add a test case for this

perturbations.jl:

• add massive neutrino initial conditions & hierarchy
• add integrals over perturbed Fermi-Dirac phase space via quadrature
• modify u/du vector to be bigger (n_q*\ell_mnu more points)

utils.jl:

• add a few functions to map unit interval for quadrature

plot_perts_k.jl

• add massive neutrino perturbations, other small edits

New file scripts/def_fd_quadrature.jl to test & compare FastGaussQudrature.jl timings and accuracies.

[xzackli]

I'm very happy to merge this when you think it's ready! It looks good, and we can just PR more when stuff comes up.

Once this is merged, I'll clean up and merge my RECFAST, and also stick helium in the Peebles recombination. I think the temperature evolution should have some effect on P(k), but we have the pieces to get it right.

[xzackli]

Regarding the instability you encountered -- let's turn that into an issue, I imagine it's easy to make an minimal demonstration. My current set of questions is

1. Does this happen if you remove the effect of massive neutrinos on the metric? (i.e. they're just free streaming the entire time)
2. Does this happen if you truncate the massive neutrinos to just the monopole / dipole?
3. Does this happen when you change the q grid?

However, we can deal with it later.

[xzackli]

Here's the sparsity pattern of hierarchy!. After merging, I'll probably move all of the standard model components to the beginning of the state vector. Also, I think we can reduce the bandwidth by reordering the neutrinos so they depend on consecutive elements. I think that will let us use one something like a BlockBandedMatrix.

Does this sparsity pattern seem correct to you?

Here's the script I used to generate this. It takes a while though. You can use spy from Plots to view it.

using Bolt
using BenchmarkTools
using OrdinaryDiffEq

par = CosmoParams()
bg = Background(par)
ih = IonizationHistory(Peebles(), par, bg)

k_grid = quadratic_k(0.1bg.H₀, 1000bg.H₀, 100)
k_i = 10
hierarchy = Hierarchy(BasicNewtonian(), par, bg, ih, k_grid[k_i])
xᵢ = (hierarchy.bg.x_grid)[1]
u₀ = Bolt.initial_conditions(xᵢ, hierarchy)
du_dummy = deepcopy(u₀)

using SparsityDetection, SparseArrays
hf(du, u) = Bolt.hierarchy!(du, u, hierarchy, xᵢ)
sparsity_pattern = jacobian_sparsity(hf, du_dummy, u₀)

[codecov-io]

Codecov Report

Merging #39 (4b918cb) into main (604d535) will decrease coverage by 3.91%. The diff coverage is 0.00%.

@@            Coverage Diff             @@
##             main      #39      +/-   ##
==========================================
- Coverage   18.08%   14.16%   -3.92%     
==========================================
  Files           7        7              
  Lines         271      346      +75     
==========================================
  Hits           49       49              
- Misses        222      297      +75     

Impacted Files	Coverage Δ
src/Bolt.jl	100.00% <ø> (ø)
src/background.jl	0.00% <0.00%> (ø)
src/perturbations.jl	0.00% <0.00%> (ø)
src/util.jl	0.00% <0.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 604d535...4b918cb. Read the comment docs.

[jmsull]

Okay I have never looked at one of these before, but I skimmed the paper linked on the github page so let me see if I have it right (sorry I am sort of doing my learning on the page).

In your call to sparsity_pattern, you are seeing whether a particular row of the output (du) depends on a particular column of the input (u). If it doesn't, those are always zero in the sparsity pattern.

I made a smaller version (decrease the number of multipoles and q to smallest nontrivial case) of this plot to get a feel for this (I added a dashed line to denote the multipole divisions so each dashed line marks nq).

So e.g. in the first row we have that \Theta[0]' depends on \Theta[1] and \Phi', which itself depends on {\Theta[0], \delta, \delta_b, \scrN[0], \scrM[0]->\scrM[nq-1]} and (00 eqn), as well as {\Phi, \Theta[2], \scrN[2], \scrM[2nq]->\scrM[3nq-1]} (through \Psi).

Going to a version of your plot with the realistic size of u (basically your example but \ell\gamma=\ell\nu):

For the massless neutrinos, we have a similar pattern as for photons which looks right just reading off the hierarchy. For massive neutrinos, first looking at the big blocks corresponding to the monopole, dipole, and quadrupole:

• The first of these (top left) show that the monopole (\scrM'[0]->\scrM'[nq-1]) depends on all monopole q-contributions (through \Phi` and therefore the 00-eqn which integrates the monopole over q) as well as the q-diagonals of the dipole contributions, and all q in the quadrupole contributions (similar to the monopole, but through \Phi). We also get the dependence on all monopoles of other species through \Phi' and the relativistic quadrupoles through \Psi.
• The second block shows that the dipole depends on diagonal q-contributions of the monopole and the total integrated quadrupole through (\Psi). We also get the dependence on the relativistic quadrupoles and \Phi through \Psi.

From the quadrupole onward, we get the long skinny band encoding the fact that each multipole at each q only receives contributions from neighboring multipoles at the same q. At the very end we get a shift in the band corresponding to the truncation.

^All of the above behavior lines up with what I would expect from staring at the equations

So the short answer to your question is yes - that looks like what I would expect for the neutrinos.

---

Also, I was looking at the bottom right corner (the \Phi, cdm, baryon subsector) to check my understanding and found a hole. When I look at the MB expression (eqn 67) for the baryon velocity, there is a speed of sound term that is proportional to delta_b that is missing in Callin/Dodelson (and Bolt). Shouldn't we have that? It seems to be in class and camb (some mention at the bottom of this issue).