Neutrino mass hierarchies#
# import necessary modules
# uncomment to get plots displayed in notebook
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from classy import Class
from scipy.optimize import fsolve
# esthetic definitions for the plots
font = {'size' : 16, 'family':'STIXGeneral'}
axislabelfontsize='large'
matplotlib.rc('font', **font)
matplotlib.rcParams['legend.fontsize']='medium'
# a function returning the three masses given the Delta m^2, the total mass, and the hierarchy (e.g. 'IN' or 'IH')
# taken from a piece of MontePython written by Thejs Brinckmann
def get_masses(delta_m_squared_atm, delta_m_squared_sol, sum_masses, hierarchy):
# any string containing letter 'n' will be considered as refering to normal hierarchy
if 'n' in hierarchy.lower():
# Normal hierarchy massive neutrinos. Calculates the individual
# neutrino masses from M_tot_NH and deletes M_tot_NH
#delta_m_squared_atm=2.45e-3
#delta_m_squared_sol=7.50e-5
m1_func = lambda m1, M_tot, d_m_sq_atm, d_m_sq_sol: M_tot**2. + 0.5*d_m_sq_sol - d_m_sq_atm + m1**2. - 2.*M_tot*m1 - 2.*M_tot*(d_m_sq_sol+m1**2.)**0.5 + 2.*m1*(d_m_sq_sol+m1**2.)**0.5
m1,opt_output,success,output_message = fsolve(m1_func,sum_masses/3.,(sum_masses,delta_m_squared_atm,delta_m_squared_sol),full_output=True)
m1 = m1[0]
m2 = (delta_m_squared_sol + m1**2.)**0.5
m3 = (delta_m_squared_atm + 0.5*(m2**2. + m1**2.))**0.5
return m1,m2,m3
else:
# Inverted hierarchy massive neutrinos. Calculates the individual
# neutrino masses from M_tot_IH and deletes M_tot_IH
#delta_m_squared_atm=-2.45e-3
#delta_m_squared_sol=7.50e-5
delta_m_squared_atm = -delta_m_squared_atm
m1_func = lambda m1, M_tot, d_m_sq_atm, d_m_sq_sol: M_tot**2. + 0.5*d_m_sq_sol - d_m_sq_atm + m1**2. - 2.*M_tot*m1 - 2.*M_tot*(d_m_sq_sol+m1**2.)**0.5 + 2.*m1*(d_m_sq_sol+m1**2.)**0.5
m1,opt_output,success,output_message = fsolve(m1_func,sum_masses/3.,(sum_masses,delta_m_squared_atm,delta_m_squared_sol),full_output=True)
m1 = m1[0]
m2 = (delta_m_squared_sol + m1**2.)**0.5
m3 = (delta_m_squared_atm + 0.5*(m2**2. + m1**2.))**0.5
return m1,m2,m3
# test of this function, returning the 3 masses for total mass of 0.1eV
m1,m2,m3 = get_masses(2.45e-3,7.50e-5,0.1,'NH')
print('NH:',m1,m2,m3,m1+m2+m3)
m1,m2,m3 = get_masses(2.45e-3,7.50e-5,0.1,'IH')
print('IH:',m1,m2,m3,m1+m2+m3)
NH: 0.02193372404432645 0.023581523497277785 0.05448475245839577 0.1
IH: 0.04912714706993182 0.04988463269615915 0.0009882202339108895 0.10000000000000185
# The goal of this cell is to compute the ratio of P(k) for NH and IH with the same total mass
commonsettings = {'N_ur':0,
'N_ncdm':3,
'output':'mPk',
'P_k_max_1/Mpc':3.0,
# The next line should be uncommented fgor higher precision (but significantly slower running)
'ncdm_fluid_approximation':3,
# You may uncomment this line to get more info on the ncdm sector from Class:
'background_verbose':1
}
# array of k values in 1/Mpc
kvec = np.logspace(-4,np.log10(3),100)
# array for storing legend
legarray = []
# loop over total mass values
for sum_masses in [0.1, 0.115, 0.13]:
# normal hierarchy
[m1, m2, m3] = get_masses(2.45e-3,7.50e-5, sum_masses, 'NH')
NH = Class()
NH.set(commonsettings)
NH.set({'m_ncdm':str(m1)+','+str(m2)+','+str(m3)})
NH.compute()
# inverted hierarchy
[m1, m2, m3] = get_masses(2.45e-3,7.50e-5, sum_masses, 'IH')
IH = Class()
IH.set(commonsettings)
IH.set({'m_ncdm':str(m1)+','+str(m2)+','+str(m3)})
IH.compute()
pkNH = []
pkIH = []
for k in kvec:
pkNH.append(NH.pk(k,0.))
pkIH.append(IH.pk(k,0.))
NH.struct_cleanup()
IH.struct_cleanup()
# extract h value to convert k from 1/Mpc to h/Mpc
h = NH.h()
plt.semilogx(kvec/h,1-np.array(pkNH)/np.array(pkIH))
legarray.append(r'$\Sigma m_i = '+str(sum_masses)+'$eV')
plt.axhline(0,color='k')
plt.xlim(kvec[0]/h,kvec[-1]/h)
plt.xlabel(r'$k [h \mathrm{Mpc}^{-1}]$')
plt.ylabel(r'$1-P(k)^\mathrm{NH}/P(k)^\mathrm{IH}$')
plt.legend(legarray)
Running CLASS version 26.0.0
Computing background
-> ncdm species sampled with 9 (resp. 4) points for purpose of background (resp. perturbation) integration. In the relativistic limit it gives Delta N_eff = 1.0132
-> ncdm species sampled with 9 (resp. 4) points for purpose of background (resp. perturbation) integration. In the relativistic limit it gives Delta N_eff = 1.0132
-> ncdm species sampled with 9 (resp. 4) points for purpose of background (resp. perturbation) integration. In the relativistic limit it gives Delta N_eff = 1.0132
-> total N_eff = 3.0396 (sumed over ultra-relativistic species, ncdm and dark radiation)
-> non-cold dark matter species has m_i = 2.193372e-02 eV (so m_i / omega_i = 9.310794e+01 eV)
-> non-cold dark matter species has m_i = 2.358152e-02 eV (so m_i / omega_i = 9.311271e+01 eV)
-> non-cold dark matter species has m_i = 5.448475e-02 eV (so m_i / omega_i = 9.313760e+01 eV)
-> age = 13.768062 Gyr
-> conformal age = 14132.013580 Mpc
-> radiation/matter equality at z = 3407.109348
corresponding to conformal time = 112.715943 Mpc
Running CLASS version 26.0.0
Computing background
-> ncdm species sampled with 9 (resp. 4) points for purpose of background (resp. perturbation) integration. In the relativistic limit it gives Delta N_eff = 1.0132
-> ncdm species sampled with 9 (resp. 4) points for purpose of background (resp. perturbation) integration. In the relativistic limit it gives Delta N_eff = 1.0132
-> ncdm species sampled with 9 (resp. 4) points for purpose of background (resp. perturbation) integration. In the relativistic limit it gives Delta N_eff = 1.0132
-> total N_eff = 3.0396 (sumed over ultra-relativistic species, ncdm and dark radiation)
-> non-cold dark matter species has m_i = 4.912715e-02 eV (so m_i / omega_i = 9.313628e+01 eV)
-> non-cold dark matter species has m_i = 4.988463e-02 eV (so m_i / omega_i = 9.313649e+01 eV)
-> non-cold dark matter species has m_i = 9.882202e-04 eV (so m_i / omega_i = 8.020738e+01 eV)
-> age = 13.767723 Gyr
-> conformal age = 14130.813058 Mpc
-> radiation/matter equality at z = 3407.381684
corresponding to conformal time = 112.707986 Mpc
Running CLASS version 26.0.0
Computing background
-> ncdm species sampled with 9 (resp. 4) points for purpose of background (resp. perturbation) integration. In the relativistic limit it gives Delta N_eff = 1.0132
-> ncdm species sampled with 9 (resp. 4) points for purpose of background (resp. perturbation) integration. In the relativistic limit it gives Delta N_eff = 1.0132
-> ncdm species sampled with 9 (resp. 4) points for purpose of background (resp. perturbation) integration. In the relativistic limit it gives Delta N_eff = 1.0132
-> total N_eff = 3.0396 (sumed over ultra-relativistic species, ncdm and dark radiation)
-> non-cold dark matter species has m_i = 2.820203e-02 eV (so m_i / omega_i = 9.312192e+01 eV)
-> non-cold dark matter species has m_i = 2.950177e-02 eV (so m_i / omega_i = 9.312376e+01 eV)
-> non-cold dark matter species has m_i = 5.729620e-02 eV (so m_i / omega_i = 9.313815e+01 eV)
-> age = 13.763758 Gyr
-> conformal age = 14126.460951 Mpc
-> radiation/matter equality at z = 3407.395810
corresponding to conformal time = 112.707573 Mpc
Running CLASS version 26.0.0
Computing background
-> ncdm species sampled with 9 (resp. 4) points for purpose of background (resp. perturbation) integration. In the relativistic limit it gives Delta N_eff = 1.0132
-> ncdm species sampled with 9 (resp. 4) points for purpose of background (resp. perturbation) integration. In the relativistic limit it gives Delta N_eff = 1.0132
-> ncdm species sampled with 9 (resp. 4) points for purpose of background (resp. perturbation) integration. In the relativistic limit it gives Delta N_eff = 1.0132
-> total N_eff = 3.0396 (sumed over ultra-relativistic species, ncdm and dark radiation)
-> non-cold dark matter species has m_i = 5.074946e-02 eV (so m_i / omega_i = 9.313672e+01 eV)
-> non-cold dark matter species has m_i = 5.148308e-02 eV (so m_i / omega_i = 9.313691e+01 eV)
-> non-col
<matplotlib.legend.Legend at 0x7cdbdb347a10>
plt.savefig('neutrinohierarchy.pdf')
<Figure size 640x480 with 0 Axes>