jmsull
commented
3 years ago After too many hours, I found that this was due to a problem in the initial conditions. The ICs as they stand are I believe(?) copied out of Callin06, and even though he includes the massless neutrino hierarchy, the ICs do not consistently account for massless neutrinos.
The reason for this is that when neutrinos are present it is not correct to use both \Theta{0} = \scrN{0} = 1.0/2 and \Phi=1 in the initial conditions (as the code was doing). If \Phi=1, this fixes \Psi through \Phi = -\Psi(1+2/5 f_{\nu}). But the initial value of the photon density perturbation and all the other (adiabatic) perturbations are determined through \Psi, not \Phi. Ignoring this fact is fine when massless neutrinos are not included, and then the only relativistic species is photons. Since photons have negligible quadrupole at early times, turning Neff down appears to fix the problem by sending \Psi to 1 (= \Phi).
It looks like this mistake may have been made in Callin because the discussion around this in Dodelson (the reference he cites for ICs) is not very careful with respect to the neutrinos. The statement that \Theta{0} = \scrN{0} = \Phi/2 is implicitly assuming \Phi=\Psi - i.e. it assumes the non-existence of neutrinos.
If I instead use the MB ICs (eqn 98), which use \Psi to set the density ICs, as expected from the above discussion, we resolve this issue (no more evolution and \Phi' \approx 0):
I will add something describing this in the notes.
In trying to match the CLASS perturbations in #44, we saw before that there was a several percent disagreement pretty much all the time (i.e. even on superhorizon scales / early times). I have looked at this more carefully, and \Phi is clearly doing something strange at very early times (plot attached for k=0.03 h/Mpc) that we were not looking at in that draft PR thread.
The orange dashed line is Bolt and shows that at early times \Phi is evolving in the first 100 steps or so before leveling off until horizon entry. This definitely should not be happening - \Phi (and \Psi) should be constant (and \Phi' negligible) until horizon entry (e.g. MB Fig.2).
You can actually see that initially (very first step) we are in agreement with the initial \Phi value of CLASS, but after this evolution, our \Phi is too high and this affects subsequent evolution, resulting in the percent level differences we see (I am surprised it is not worse, to be honest).
I tried setting Neff to some very small value (1e-4) and this seems to make the problem go away (though late-time evolution won't match CLASS because Neff is wrong):
This is also true if I take out both massive neutrinos and massless neutrinos (by commenting out the lines in the background and einstein equations). The issue is not related to the massive neutrinos, as it persists when only massless neutrinos are included.
I likely made some mistake in the past when setting up the massless neutrinos. I will dig into this in the first few steps to see what about massless neutrinos is making \Phi' non-negligible. I have checked that we get the correct \Psi in the initial conditions, so I guess it is somehow an error in the 00 Einstein eqn (or something less obvious with the neutrinos).