[BUG] neutrino mass incorrectly implemented

[liuxx479]

Describe the bug
Neutrino mass parameter is not properly modeled right now, because astropy takes in m_nu (in unit of eV), but camb uses omnuh2 (no unit).

To Reproduce
Steps to reproduce the behavior:

from astropy import units as u
tr0 = Transfer(cosmo_params={"m_nu":u.eV*np.array([0,0,0]) })
tr1 = Transfer(cosmo_params={"m_nu":u.eV*np.array([0.1, 0.1, 0.1]) } )

tr0.transfer_function - tr1.transfer_function

would show 0, meaning massive and massless neutrinos have identical transfer functions, because m_nu parameter is not passed on to camb

Expected behavior
tr0.transfer_function - tr1.transfer_function should be different

Additional context
I have a quick fix by redefining the CDM parameter omch2:

change line 212 in hmf/density_field/transfer_models.py from:
omch2=(self.cosmo.Om0 - self.cosmo.Ob0) * self.cosmo.h ** 2
to:
omch2=(self.cosmo.Om0 - self.cosmo.Ob0) * self.cosmo.h ** 2 - sum(self.cosmo.m_nu.value)/93.14

But I think there should be a better way to do it.

[steven-murray]

Thanks for the report and the potential fix! When passing on to CAMB, does one also need to specify omnuh2 (as well as taking away its contribution to omch2)? Otherwise, this seems like a perfectly reasonable fix. What do you think could be better? Would be very happy with a PR along these lines :-)

…

On Fri, Dec 11, 2020 at 7:29 PM Jia Liu ***@***.***> wrote: Assigned #117 <#117> to @steven-murray <https://github.com/steven-murray>. — You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <#117 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAJWRXSZAHVXNSSRWLTOTQDSULISHANCNFSM4UXSRDDA> .

[liuxx479]

Hi Steven,
camb requires neutrino hierarchy, number of massive neutrinos, and total neutrino mass and compute the rest automatically:
e.g.
neutrino_hierarchy: Union[str, int] = 'degenerate', num_massive_neutrinos=1,mnu=0.06,
And yes, you need to manually calculate omch2 (in our case, subtracting massive neutrino by hand).

I did a bit editing from my original proposal - instead of manually calculate omega_nu, I use the calculation directly from astropy. I fixed the splitting to degenerate since it should be sufficient for any realistic usage.. I submitted a PR.

A more complicated fix (any neutrino mass values and splitting) can probably be done with more care, following the description under set_cosmology in https://camb.readthedocs.io/en/latest/_modules/camb/model.html#CAMBparams

[liuxx479]

After the first attempt to fix the issue (but fixing the input parameter for camb), my colleague brought to my attention that for neutrinos, the best match of HMF (to simulations) is achieved using P_cb and rho_cb (cb=cdm+baryons) during calculation (though the background cosmology does take into account of massive neutrinos), instead of P_matter and rho_matter (matter = cdm+baryon+neutrinos). See fig. 1 of: https://arxiv.org/abs/1311.1514

Implementation of this is more involved.. I will come back to it later..

[steven-murray]

Note that #280 updated the docs and tests around this issue. However, @liuxx479 's final comment pointing to the paper using cdm instead of cdm + baryons has not been implemented. I do not think of this as a bug, but a potential gain in accuracy.