fix for bug with very small but non-zero neutrino masses

fortran/MathUtils.f90

@@ -123,7 +123,7 @@ subroutine brentq(obj,func,ax,bx,tol,xzero,fzero,iflag,fax,fbx)
     real(dl),intent(out)             :: fzero   !! value of `f` at the root (`f(xzero)`)
     integer,intent(out)              :: iflag   !! status flag (`-1`=error, `0`=root found)
     real(dl),intent(in),optional     :: fax     !! if `f(ax)` is already known, it can be input here
-    real(dl),intent(in),optional     :: fbx     !! if `f(ax)` is already known, it can be input here
+    real(dl),intent(in),optional     :: fbx     !! if `f(bx)` is already known, it can be input here
     real(dl), parameter :: one = 1._dl, zero = 0._dl, two =2._dl, three = 3._dl
     real(dl),parameter :: eps   = epsilon(one)  !! original code had d1mach(4)
     real(dl) :: a,b,c,d,e,fa,fb,fc,tol1,xm,p,q,r,s

fortran/massive_neutrinos.f90

@@ -31,11 +31,13 @@ module MassiveNu
         !Quantities for the neutrino background momentum distribution assuming thermal
         real(dl) dam !step in a*m
         real(dl), dimension(:), allocatable ::  r1,p1,dr1,dp1,ddr1
+        real(dl), private :: target_rho
     contains
     procedure :: init => ThermalNuBackground_init
     procedure :: rho_P => ThermalNuBackground_rho_P
     procedure :: rho => ThermalNuBackground_rho
     procedure :: drho => ThermalNuBackground_drho
+    procedure :: find_nu_mass_for_rho => ThermalNuBackground_find_nu_mass_for_rho
     end type TThermalNuBackground
 
     Type(TThermalNuBackground), target :: ThermalNuBackground
@@ -218,8 +220,8 @@ subroutine ThermalNuBackground_rho(this,am,rhonu)
     integer i
 
     !  Compute massive neutrino density in units of the mean
-    !  density of one eigenstate of massless neutrinos.  Use cubic splines to
-    !  interpolate from a table.
+    !  density of one eigenstate of massless neutrinos.  Use series solutions or
+    !  cubic splines to interpolate from a table.
 
     if (am <= am_minp) then
         if (am < 0.01_dl) then
@@ -244,6 +246,58 @@ subroutine ThermalNuBackground_rho(this,am,rhonu)
 
     end subroutine ThermalNuBackground_rho
 
+    function rho_err(this, nu_mass)
+    class(TThermalNuBackground) :: this
+    real(dl) rho_err, nu_mass, rhonu
+
+    call this%rho(nu_mass, rhonu)
+    rho_err = rhonu - this%target_rho
+
+    end function rho_err
+
+    function ThermalNuBackground_find_nu_mass_for_rho(this,rho) result(nu_mass)
+    !  Get eigenstate mass given input density (rho is neutrino density in units of one massless)
+    !  nu_mass=m_n*c**2/(k_B*T_nu0).
+    !  Get number density n of neutrinos from
+    !  rho_massless/n = int q^3/(1+e^q) / int q^2/(1+e^q)=7/180 pi^4/Zeta(3)
+    !  then m = Omega_nu/N_nu rho_crit /n if non-relativistic
+    use MathUtils
+    use config
+    class(TThermalNuBackground) :: this
+    real(dl), intent(in) :: rho
+    real(dl) nu_mass, rhonu, rhonu1, delta
+    real(dl) fzero
+    integer iflag
+
+    if (rho <= 1.001_dl) then
+        !energy density all accounted for by massless result
+        nu_mass=0
+    else
+        !Get mass assuming fully non-relativistic
+        nu_mass=fermi_dirac_const/(1.5d0*zeta3)*rho
+
+        if (nu_mass>4) then
+            !  perturbative correction for velocity when nearly non-relativistic
+            !  Error due to velocity < 1e-5 for mnu~0.06 but can easily correct (assuming non-relativistic today)
+            !  Note that python does not propagate mnu to omnuh2 consistently to the same accuracy; but this makes
+            !  fortran more internally consistent between input and computed Omega_nu h^2
+
+            !Make perturbative correction for the tiny error due to the neutrino velocity
+            call this%rho(nu_mass, rhonu)
+            call this%rho(nu_mass*0.9, rhonu1)
+            delta = rhonu - rho
+            nu_mass = nu_mass*(1 + delta/((rhonu1 - rhonu)/0.1) )
+        else
+            !Directly solve to avoid issues with perturbative result when no longer very relativistic
+            this%target_rho = rho
+            call brentq(this,rho_err,0._dl,nu_mass,0.01_dl,nu_mass,fzero,iflag)
+            if (iflag/=0) call GlobalError('find_nu_mass_for_rho failed to find neutrino mass')
+        end if
+    end if
+
+    end function ThermalNuBackground_find_nu_mass_for_rho
+
+
     function ThermalNuBackground_drho(this,am,adotoa) result (rhonudot)
     !  Compute the time derivative of the mean density in massive neutrinos
     class(TThermalNuBackground) :: this

fortran/results.f90

@@ -453,23 +453,10 @@ subroutine CAMBdata_SetParams(this, P, error, DoReion, call_again)
             !Initialize things for massive neutrinos
             call ThermalNuBackground%Init()
             call this%NuPerturbations%Init(P%Accuracy%AccuracyBoost*P%Accuracy%neutrino_q_boost)
-            !  nu_masses=m_nu(i)*c**2/(k_B*T_nu0).
-            !  Get number density n of neutrinos from
-            !  rho_massless/n = int q^3/(1+e^q) / int q^2/(1+e^q)=7/180 pi^4/Zeta(3)
-            !  then m = Omega_nu/N_nu rho_crit /n
-            !  Error due to velocity < 1e-5 for mnu~0.06 but can easily correct
+            !  nu_masses=m_nu(i)*c**2/(k_B*T_nu0)
             do nu_i=1, this%CP%Nu_mass_eigenstates
-                this%nu_masses(nu_i)=fermi_dirac_const/(1.5d0*zeta3)*this%grhocrit/this%grhor* &
-                    this%CP%omnuh2/h2*this%CP%Nu_mass_fractions(nu_i)/this%CP%Nu_mass_degeneracies(nu_i)
-                block
-                    real(dl) rhonu, rhonu1, delta
-
-                    !Make perturbative correction for the tiny error due to the neutrino velocity
-                    call ThermalNuBackground%rho(this%nu_masses(nu_i), rhonu)
-                    call ThermalNuBackground%rho(this%nu_masses(nu_i)*0.9, rhonu1)
-                    delta = rhonu - this%CP%Nu_mass_fractions(nu_i)*this%grhocrit*this%CP%omnuh2/h2/this%grhormass(nu_i)
-                    this%nu_masses(nu_i) = this%nu_masses(nu_i)*(1 + delta/((rhonu1 - rhonu)/0.1) )
-                end block
+                this%nu_masses(nu_i)= ThermalNuBackground%find_nu_mass_for_rho(this%CP%omnuh2/h2*this%CP%Nu_mass_fractions(nu_i)&
+                    *this%grhocrit/this%grhormass(nu_i))
             end do
         else
             this%nu_masses = 0