omega_nu_single.h File Reference

#include <gsl/gsl_interp.h>

Include dependency graph for omega_nu_single.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	_rho_nu_single
struct	_hybrid_nu
struct	_omega_nu

Macros
#define	TNUCMB   (pow(4/11.,1/3.)*1.00328)
#define	NUSPECIES   3

Typedefs
typedef struct _rho_nu_single	_rho_nu_single
typedef struct _hybrid_nu	_hybrid_nu
typedef struct _omega_nu	_omega_nu

Functions
void	rho_nu_init (_rho_nu_single *rho_nu_tab, double a0, const double mnu, const double kBtnu)
double	rho_nu (const _rho_nu_single *rho_nu_tab, const double a, const double kT)
void	init_hybrid_nu (_hybrid_nu *const hybnu, const double mnu[], const double vcrit, const double light, const double nu_crit_time, const double kBtnu)
double	particle_nu_fraction (const _hybrid_nu *const hybnu, const double a, int i)
double	nufrac_low (const double qc)
void	init_omega_nu (_omega_nu *const omnu, const double MNu[], const double a0, const double HubbleParam, const double tcmb0)
double	get_omega_nu (const _omega_nu *const omnu, const double a)
double	get_omega_nu_nopart (const _omega_nu *const omnu, const double a)
double	get_omegag (const _omega_nu *const omnu, const double a)
double	omega_nu_single (const _omega_nu *const rho_nu_tab, const double a, const int i)

Detailed Description

Routines for computing the matter density in a single neutrino species

Definition in file omega_nu_single.h.

Macro Definition Documentation

NUSPECIES

#define NUSPECIES   3

Number of massive neutrino species: 3 Could be made configurable at some point Neutrino masses are in eV

Definition at line 20 of file omega_nu_single.h.

TNUCMB

#define TNUCMB   (pow(4/11.,1/3.)*1.00328)

Ratio between the massless neutrino temperature and the CMB temperature. Note there is a slight correction from 4/11 due to the neutrinos being slightly coupled at e+- annihilation. See Mangano et al 2005 (hep-ph/0506164) We use the CLASS default value, chosen so that omega_nu = m_nu / 93.14 h^2 At time of writing this is T_nu / T_gamma = 0.71611. See https://github.com/lesgourg/class_public/blob/master/explanatory.ini

Definition at line 16 of file omega_nu_single.h.

Typedef Documentation

_hybrid_nu

typedef struct _hybrid_nu _hybrid_nu

Definition at line 47 of file omega_nu_single.h.

_omega_nu

typedef struct _omega_nu _omega_nu

Definition at line 89 of file omega_nu_single.h.

_rho_nu_single

typedef struct _rho_nu_single _rho_nu_single

Definition at line 1 of file omega_nu_single.h.

Function Documentation

get_omega_nu()

double get_omega_nu	(	const _omega_nu *const	omnu,
		const double	a
	)

Return the total matter density in neutrinos at scale factor a.

Definition at line 57 of file omega_nu_single.c.

{
        double rhonu=0;
        int mi;
        for(mi=0; mi<NUSPECIES; mi++) {
            if(omnu->nu_degeneracies[mi] > 0){
                 rhonu += omnu->nu_degeneracies[mi] * rho_nu(&omnu->RhoNuTab[mi], a, omnu->kBtnu);
            }
        }
        return rhonu/omnu->rhocrit;
}

References _omega_nu::kBtnu, _omega_nu::nu_degeneracies, NUSPECIES, rho_nu(), _omega_nu::rhocrit, and _omega_nu::RhoNuTab.

Referenced by check_units(), compute_mass(), delta_nu_from_power(), delta_tot_first_init(), get_long_range_timestep_dloga(), get_omega_nu_nopart(), growth(), hubble_function(), init_cosmology(), init_neutrinos_lra(), init_transfer_table(), main(), parse_transfer(), test_get_omega_nu(), and update_delta_tot().

Here is the call graph for this function:

Here is the caller graph for this function:

get_omega_nu_nopart()

double get_omega_nu_nopart	(	const _omega_nu *const	omnu,
		const double	a
	)

Return the total matter density in neutrinos at scale factor a , excluding active particles.

Definition at line 71 of file omega_nu_single.c.

{
    double omega_nu = get_omega_nu(omnu, a);
    double part_nu = get_omega_nu(omnu, 1) * particle_nu_fraction(&omnu->hybnu, a, 0) / (a*a*a);
    return omega_nu - part_nu;
}

References get_omega_nu(), _omega_nu::hybnu, and particle_nu_fraction().

Referenced by check_omega(), delta_nu_from_power(), delta_tot_first_init(), get_delta_nu_combined(), potential_transfer(), test_hybrid_neutrinos(), and update_delta_tot().

Here is the call graph for this function:

Here is the caller graph for this function:

get_omegag()

double get_omegag	(	const _omega_nu *const	omnu,
		const double	a
	)

Return the photon matter density at scale factor a

Definition at line 79 of file omega_nu_single.c.

{
    const double omegag = 4*STEFAN_BOLTZMANN/(LIGHTCGS*LIGHTCGS*LIGHTCGS)*pow(omnu->tcmb0,4)/omnu->rhocrit;
    return omegag/pow(a,4);
}

References LIGHTCGS, _omega_nu::rhocrit, STEFAN_BOLTZMANN, and _omega_nu::tcmb0.

Referenced by test_get_omegag().

Here is the caller graph for this function:

init_hybrid_nu()

void init_hybrid_nu	(	_hybrid_nu *const	hybnu,
		const double	mnu[],
		const double	vcrit,
		const double	light,
		const double	nu_crit_time,
		const double	kBtnu
	)

Set up parameters for the hybrid neutrinos

Parameters

hybnu	initialised structure
mnu	array of neutrino masses in eV
vcrit	Critical velocity above which to treat neutrinos with particles. Note this is unperturbed velocity TODAY To get velocity at redshift z, multiply by (1+z)
light	speed of light in internal units
nu_crit_time	critical time to make neutrino particles live
kBtnu	Boltzmann constant times neutrino temperature. Dimensionful factor.

Definition at line 235 of file omega_nu_single.c.

{
    hybnu->enabled=1;
    int i;
    hybnu->nu_crit_time = nu_crit_time;
    hybnu->vcrit = vcrit / light;
    for(i=0; i< NUSPECIES; i++) {
        const double qc = mnu[i] * vcrit / light / kBtnu;
        hybnu->nufrac_low[i] = nufrac_low(qc);
    }
}

References _hybrid_nu::enabled, _hybrid_nu::nu_crit_time, nufrac_low(), _hybrid_nu::nufrac_low, NUSPECIES, and _hybrid_nu::vcrit.

Referenced by init_cosmology(), and test_hybrid_neutrinos().

Here is the call graph for this function:

Here is the caller graph for this function:

init_omega_nu()

void init_omega_nu	(	_omega_nu *const	omnu,
		const double	MNu[],
		const double	a0,
		const double	HubbleParam,
		const double	tcmb0
	)

Initialise the neutrino structure, do the time integration and allocate memory for the subclass rho_nu_single

Parameters

omnu	structure to initialise
MNu	array of neutrino masses in eV. Three entries.
a0	initial scale factor.
HubbleParam	Hubble parameter h0, eg, 0.7.
tcmb0	Redshift zero CMB temperature.

Definition at line 20 of file omega_nu_single.c.

{
    int mi;
    /*Explicitly disable hybrid neutrinos*/
    omnu->hybnu.enabled=0;
    /*CMB temperature*/
    omnu->tcmb0 = tcmb0;
    /*Neutrino temperature times k_B*/
    omnu->kBtnu = BOLEVK * TNUCMB * tcmb0;
    /*Store conversion between rho and omega*/
    omnu->rhocrit = (3 * HUBBLE * HubbleParam * HUBBLE * HubbleParam)/ (8 * M_PI * GRAVITY);
    /*First compute which neutrinos are degenerate with each other*/
    for(mi=0; mi<NUSPECIES; mi++){
        int mmi;
        omnu->nu_degeneracies[mi]=0;
        for(mmi=0; mmi<mi; mmi++){
            if(fabs(MNu[mi] -MNu[mmi]) < FLOAT_ACC){
                omnu->nu_degeneracies[mmi]+=1;
                break;
            }
        }
        if(mmi==mi) {
            omnu->nu_degeneracies[mi]=1;
        }
    }
    /*Now allocate a table for the species we want*/
    for(mi=0; mi<NUSPECIES; mi++){
        if(omnu->nu_degeneracies[mi]) {
            rho_nu_init(&omnu->RhoNuTab[mi], a0, MNu[mi], omnu->kBtnu);
        }
        else
            omnu->RhoNuTab[mi].loga = 0;
    }
}

References BOLEVK, _hybrid_nu::enabled, FLOAT_ACC, GRAVITY, HUBBLE, _omega_nu::hybnu, _omega_nu::kBtnu, _rho_nu_single::loga, _omega_nu::nu_degeneracies, NUSPECIES, rho_nu_init(), _omega_nu::rhocrit, _omega_nu::RhoNuTab, _omega_nu::tcmb0, and TNUCMB.

Referenced by init_cosmology(), test_allocate_delta_tot_table(), test_get_omega_nu(), test_get_omegag(), test_hybrid_neutrinos(), test_omega_nu_init_degenerate(), test_omega_nu_init_nondeg(), test_omega_nu_single(), and test_omega_nu_single_exact().

Here is the call graph for this function:

Here is the caller graph for this function:

nufrac_low()

double nufrac_low

(

const double

qc

)

Integrate the fermi-dirac kernel between 0 and qc to find the fraction of neutrinos that are particles

Definition at line 219 of file omega_nu_single.c.

{
    /*These functions are so smooth that we don't need much space*/
    gsl_integration_workspace * w = gsl_integration_workspace_alloc (100);
    double abserr;
    gsl_function F;
    F.function = &fermi_dirac_kernel;
    F.params = NULL;
    double total_fd;
    gsl_integration_qag (&F, 0, qc, 0, 1e-6,100,GSL_INTEG_GAUSS61, w,&(total_fd), &abserr);
    /*divided by the total F-D probability (which is 3 Zeta(3)/2 ~ 1.8 if MAX_FERMI_DIRAC is large enough*/
    total_fd /= 1.5*1.202056903159594;
    gsl_integration_workspace_free (w);
    return total_fd;
}

References fermi_dirac_kernel().

Referenced by init_hybrid_nu(), Jfrac_high(), test_hybrid_neutrinos(), and test_nufrac_low().

Here is the call graph for this function:

Here is the caller graph for this function:

omega_nu_single()

double omega_nu_single	(	const _omega_nu *const	rho_nu_tab,
		const double	a,
		const int	i
	)

Return the matter density in a single neutrino species

Parameters

rho_nu_tab	structure containing pre-computed matter density values.
i	index of neutrino species we want
a	scale factor desired

Definition at line 267 of file omega_nu_single.c.

{
    /*Deal with case where we want a species degenerate with another one*/
    if(omnu->nu_degeneracies[i] == 0) {
        int j;
        for(j=i; j >=0; j--)
            if(omnu->nu_degeneracies[j]){
                i = j;
                break;
            }
    }
    double omega_nu = rho_nu(&omnu->RhoNuTab[i], a,omnu->kBtnu)/omnu->rhocrit;
    double omega_part = rho_nu(&omnu->RhoNuTab[i], 1,omnu->kBtnu)/omnu->rhocrit;
    omega_part *= particle_nu_fraction(&omnu->hybnu, a, i)/(a*a*a);
    omega_nu -= omega_part;
    return omega_nu;
 
}

References _omega_nu::hybnu, _omega_nu::kBtnu, _omega_nu::nu_degeneracies, particle_nu_fraction(), rho_nu(), _omega_nu::rhocrit, and _omega_nu::RhoNuTab.

Referenced by get_delta_nu_combined(), parse_transfer(), test_get_omega_nu(), test_hybrid_neutrinos(), test_omega_nu_single(), and test_omega_nu_single_exact().

Here is the call graph for this function:

Here is the caller graph for this function:

particle_nu_fraction()

double particle_nu_fraction	(	const _hybrid_nu *const	hybnu,
		const double	a,
		int	i
	)

Get fraction of neutrinos currently followed by particles.

Parameters

hybnu	Structure with hybrid neutrino parameters.
i	index of neutrino species to use.
a	redshift of interest.

Returns

the fraction of neutrinos currently traced by particles. 0 when neutrinos are fully analytic at early times.

Definition at line 251 of file omega_nu_single.c.

{
    /*Return zero if hybrid neutrinos not enabled*/
    if(!hybnu->enabled)
        return 0;
    /*Just use a redshift cut for now. Really we want something more sophisticated,
     * based on the shot noise and average overdensity.*/
    if (a > hybnu->nu_crit_time){
        return hybnu->nufrac_low[i];
    }
    return 0;
}

References _hybrid_nu::enabled, _hybrid_nu::nu_crit_time, and _hybrid_nu::nufrac_low.

Referenced by delta_nu_from_power(), delta_tot_first_init(), get_delta_nu(), get_delta_tot(), get_omega_nu_nopart(), gravpm_force(), omega_nu_single(), test_hybrid_neutrinos(), and update_delta_tot().

Here is the caller graph for this function:

rho_nu()

double rho_nu	(	const _rho_nu_single *	rho_nu_tab,
		const double	a,
		const double	kT
	)

Computes the neutrino density for a single neutrino species at a given redshift, either by looking up in a table, or a simple calculation in the limits, or by direct integration.

Parameters

rho_nu_tab	Pre-computed table of values
a	Redshift desired.
kT	Boltzmann constant times neutrino temperature. Dimensionful factor.

Returns

neutrino density in g cm^-3

Definition at line 180 of file omega_nu_single.c.

{
        double rho_nu_val;
        double amnu=a*rho_nu_tab->mnu;
        const double kTamnu2=(kT*kT/amnu/amnu);
        /*Do it analytically if we are in a regime where we can
         * The next term is 141682 (kT/amnu)^8.
         * At kT/amnu = 8, higher terms are larger and the series stops converging.
         * Don't go lower than 50 here. */
        if(NU_SW*NU_SW*kTamnu2 < 1){
            /*Heavily non-relativistic*/
            rho_nu_val = non_rel_rho_nu(a, kT, amnu, kTamnu2);
        }
        else if(amnu < 1e-6*kT){
            /*Heavily relativistic*/
            rho_nu_val=rel_rho_nu(a, kT);
        }
        else{
            const double loga = log(a);
            /* Deal with early time case. In practice no need to be very accurate
             * so assume relativistic.*/
            if (!rho_nu_tab->loga || loga < rho_nu_tab->loga[0])
                rho_nu_val = rel_rho_nu(a,kT);
            else
                rho_nu_val=gsl_interp_eval(rho_nu_tab->interp,rho_nu_tab->loga,rho_nu_tab->rhonu,loga,rho_nu_tab->acc);
        }
        return rho_nu_val;
}

References _rho_nu_single::acc, _rho_nu_single::interp, _rho_nu_single::loga, _rho_nu_single::mnu, non_rel_rho_nu(), NU_SW, rel_rho_nu(), and _rho_nu_single::rhonu.

Referenced by get_omega_nu(), and omega_nu_single().

Here is the call graph for this function:

Here is the caller graph for this function:

rho_nu_init()

void rho_nu_init	(	_rho_nu_single *	rho_nu_tab,
		double	a0,
		const double	mnu,
		const double	kBtnu
	)

Initialise the tables for matter density in a single neutrino species, by doing numerical integration.

Parameters

rho_nu_tab	Structure to initialise
a0	First scale factor in rho_nu_tab. Do not try and evaluate rho_nu at higher redshift!
mnu	neutrino mass to compute neutrino matter density for in eV
kBtnu	Boltzmann constant times neutrino temperature. Dimensionful factor.

Definition at line 119 of file omega_nu_single.c.

{
     int i;
     double abserr;
     /* Tabulate to 1e-3, unless earlier requested.
      * Need early times for growth function.*/
     if(a0 > 1e-3)
         a0 = 1e-3;
     /* Do not need a table when relativistic*/
     if(a0 * mnu < 1e-6 * kBtnu)
        a0 = 1e-6 * kBtnu / mnu;
     /*Make the table over a slightly wider range than requested, in case there is roundoff error*/
     const double logA0=log(a0)-log(1.2);
     const double logaf=log(NU_SW*kBtnu/mnu)+log(1.2);
     gsl_function F;
     F.function = &rho_nu_int;
     /*Initialise constants*/
     rho_nu_tab->mnu = mnu;
     /*Shortcircuit if we don't need to do the integration*/
     if(mnu < 1e-6*kBtnu || logaf < logA0)
         return;
 
     /*Allocate memory for arrays*/
     rho_nu_tab->loga = (double *) mymalloc("rho_nu_table",2*NRHOTAB*sizeof(double));
     rho_nu_tab->rhonu = rho_nu_tab->loga+NRHOTAB;
     rho_nu_tab->acc = gsl_interp_accel_alloc();
     rho_nu_tab->interp=gsl_interp_alloc(gsl_interp_cspline,NRHOTAB);
     if(!rho_nu_tab->interp || !rho_nu_tab->acc || !rho_nu_tab->loga)
         endrun(2035,"Could not initialise tables for neutrino matter density\n");
 
     gsl_integration_workspace * w = gsl_integration_workspace_alloc (GSL_VAL);
     for(i=0; i< NRHOTAB; i++){
        double param[2];
        rho_nu_tab->loga[i]=logA0+i*(logaf-logA0)/(NRHOTAB-1);
        param[0]=mnu*exp(rho_nu_tab->loga[i]);
        param[1] = kBtnu;
        F.params = &param;
        gsl_integration_qag (&F, 0, 500*kBtnu,0 , 1e-9,GSL_VAL,6,w,&(rho_nu_tab->rhonu[i]), &abserr);
        rho_nu_tab->rhonu[i]=rho_nu_tab->rhonu[i]/pow(exp(rho_nu_tab->loga[i]),4)*get_rho_nu_conversion();
     }
     gsl_integration_workspace_free (w);
     gsl_interp_init(rho_nu_tab->interp,rho_nu_tab->loga,rho_nu_tab->rhonu,NRHOTAB);
     return;
}

References _rho_nu_single::acc, endrun(), get_rho_nu_conversion(), GSL_VAL, _rho_nu_single::interp, _rho_nu_single::loga, _rho_nu_single::mnu, mymalloc, NRHOTAB, NU_SW, rho_nu_int(), and _rho_nu_single::rhonu.

Referenced by init_omega_nu(), and test_rho_nu_init().

Here is the call graph for this function:

Here is the caller graph for this function: