diff --git a/doc/RTD/source/ParameterFiles/parameter_description.rst b/doc/RTD/source/ParameterFiles/parameter_description.rst
index cb6248cb66a2e56b1c9d744d5ee78bf5092f6456..b41872aabfef702beaf18308cb4be917ab6553b6 100644
--- a/doc/RTD/source/ParameterFiles/parameter_description.rst
+++ b/doc/RTD/source/ParameterFiles/parameter_description.rst
@@ -171,6 +171,31 @@ w_0 + w_a (1 - a)`. The two parameters in the YAML file are:
If unspecified these parameters default to the default
:math:`\Lambda\rm{CDM}` values of :math:`w_0 = -1` and :math:`w_a = 0`.
+The radiation density :math:`\Omega_r` can also be specified by setting
+an alternative optional parameter:
+
+* The number of ultra-relativistic degrees of freedom :math:`N_\rm{ur}`:
+ ``N_ur``.
+
+The radiation density :math:`\Omega_r` is then automatically inferred from
+:math:`N_\rm{ur}` and the present-day CMB temperature
+:math:`T_{\rm{CMB},0}=2.7255` Kelvin. This parametrization cannot
+be used together with :math:`\Omega_r`. If neither parameter is used, SWIFT
+defaults to :math:`\Omega_r = 0`. Note that :math:`N_\rm{ur}` differs from
+:math:`N_\rm{eff}`, the latter of which also includes massive neutrinos.
+
+Massive neutrinos can be included by specifying the optional parameters:
+
+* The number of massive neutrino species :math:`N_{\nu}`: ``N_nu``,
+* A comma-separated list of neutrino masses in eV: ``M_nu_eV``,
+* A comma-separated list of neutrino degeneracies: ``deg_nu``,
+* The present-day neutrino temperature :math:`T_{\nu,0}`: ``T_nu_0``.
+
+When including massive neutrinos, only ``N_nu`` and ``M_nu_eV`` are necessary.
+By default, SWIFT will assume non-degenerate species and
+:math:`T_{\nu,0}=(4/11)^{1/3}T_{\rm{CMB},0}`. Neutrinos do not contribute to
+:math:`\Omega_m = \Omega_\rm{cdm} + \Omega_b` in our conventions.
+
For a Planck+13 cosmological model (ignoring radiation density as is
commonly done) and running from :math:`z=127` to :math:`z=0`, one would hence
use the following parameters:
diff --git a/examples/main.c b/examples/main.c
index 1b9c21bdb1e9f9effccd29ece728d7bfb60df3bf..620e018132edf809858f31f79e8b9c148d81ec2e 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -1189,13 +1189,17 @@ int main(int argc, char *argv[]) {
const int with_DM_background_particles =
N_total[swift_type_dark_matter_background] > 0;
+ /* Do we have neutrino particles? */
+ const int with_neutrinos = 0; // no for now
+
/* Initialize the space with these data. */
if (myrank == 0) clocks_gettime(&tic);
space_init(&s, params, &cosmo, dim, &hydro_properties, parts, gparts, sinks,
sparts, bparts, Ngas, Ngpart, Nsink, Nspart, Nbpart, periodic,
replicate, remap_ids, generate_gas_in_ics, with_hydro,
with_self_gravity, with_star_formation,
- with_DM_background_particles, talking, dry_run, nr_nodes);
+ with_DM_background_particles, with_neutrinos, talking, dry_run,
+ nr_nodes);
/* Initialise the line of sight properties. */
if (with_line_of_sight) los_init(s.dim, &los_properties, params);
diff --git a/examples/main_fof.c b/examples/main_fof.c
index bed38114d7df9cd46492b1e547295a00c95fcbd3..b40a61b531aa3b17bb281fcc5be2ad52ac06500b 100644
--- a/examples/main_fof.c
+++ b/examples/main_fof.c
@@ -542,14 +542,17 @@ int main(int argc, char *argv[]) {
const int with_DM_background_particles =
N_total[swift_type_dark_matter_background] > 0;
+ /* Do we have neutrino particles? */
+ const int with_neutrinos = 0; // no for now
+
/* Initialize the space with these data. */
if (myrank == 0) clocks_gettime(&tic);
space_init(&s, params, &cosmo, dim, /*hydro_props=*/NULL, parts, gparts,
sinks, sparts, bparts, Ngas, Ngpart, Nsink, Nspart, Nbpart,
periodic, replicate, /*remap_ids=*/0,
/*generate_gas_in_ics=*/0, /*hydro=*/N_total[0] > 0, /*gravity=*/1,
- /*with_star_formation=*/0, with_DM_background_particles, talking,
- /*dry_run=*/0, nr_nodes);
+ /*with_star_formation=*/0, with_DM_background_particles,
+ with_neutrinos, talking, /*dry_run=*/0, nr_nodes);
if (myrank == 0) {
clocks_gettime(&toc);
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index e201f4f05b28226930abd1378e69928c739ae424..075f7a02ec40081fc0738b17aedfb4166686f8f2 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -25,6 +25,11 @@ Cosmology:
Omega_r: 0. # (Optional) Radiation density parameter
w_0: -1.0 # (Optional) Dark-energy equation-of-state parameter at z=0.
w_a: 0. # (Optional) Dark-energy equation-of-state time evolution parameter.
+ T_nu_0: 1.9514 # (Optional) Present-day neutrino temperature in internal units, used for massive neutrinos
+ N_ur: 1.0196 # (Optional) Number of ultra-relativistic species (e.g. massless neutrinos). Assumes instantaneous decoupling T_ur/T_g=(4/11)^(1/3). Cannot be used together with Omega_r.
+ N_nu: 2 # (Optional) Integer number of massive neutrinos. Note that neutrinos do NOT contribute to Omega_m = Omega_cdm + Omega_b in our conventions.
+ M_nu_eV: 0.05, 0.01 # (Optional) Comma-separated list of N_nu nonzero neutrino masses in electron-volts
+ deg_nu: 1.0, 1.0 # (Optional) Comma-separated list of N_nu neutrino degeneracies (default values of 1.0)
# Parameters for the hydrodynamics scheme
SPH:
@@ -128,7 +133,6 @@ Scheduler:
engine_max_parts_per_ghost: 1000 # (Optional) Maximum number of parts per ghost.
engine_max_sparts_per_ghost: 1000 # (Optional) Maximum number of sparts per ghost.
engine_max_parts_per_cooling: 10000 # (Optional) Maximum number of parts per cooling task.
- dependency_graph_frequency: 0 # (Optional) Dumping frequency of the dependency graph. By default, writes only at the first step.
# Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.)
TimeIntegration:
@@ -169,7 +173,7 @@ Logger:
basename: index # Common part of the filenames
initial_buffer_size: 1 # (Optional) Buffer size in GB
buffer_scale: 10 # (Optional) When buffer size is too small, update it with required memory times buffer_scale
-
+
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-2 # Time between statistics output
@@ -234,7 +238,7 @@ DomainDecomposition:
StructureFinding:
config_file_name: stf_input.cfg # Name of the STF config file.
basename: haloes # Common part of the name of output files.
- subdir_per_output: stf # (Optional) Sub-directory (within Snapshots:subdir) to use for each invocation of STF. Defaults to "" (i.e. no sub-directory)
+ subdir_per_output: stf # (Optional) Sub-directory (within Snapshots:subdir) to use for each invocation of STF. Defaults to "" (i.e. no sub-directory)
scale_factor_first: 0.92 # (Optional) Scale-factor of the first structure finding (cosmological run)
time_first: 0.01 # (Optional) Time of the first structure finding output (in internal units).
delta_time: 1.10 # (Optional) Time difference between consecutive structure finding outputs (in internal units) in simulation time intervals.
@@ -283,7 +287,7 @@ EoS:
# Point mass external potentials
PointMassPotential:
- useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
+ useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
position: [50.,50.,50.] # location of external point mass (internal units)
mass: 1e10 # mass of external point mass (internal units)
timestep_mult: 0.03 # Dimensionless pre-factor for the time-step condition
@@ -291,15 +295,15 @@ PointMassPotential:
# Isothermal potential parameters
IsothermalPotential:
- useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
+ useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
position: [100.,100.,100.] # Location of centre of isothermal potential with respect to centre of the box (internal units)
vrot: 200. # Rotation speed of isothermal potential (internal units)
timestep_mult: 0.03 # Dimensionless pre-factor for the time-step condition
epsilon: 0.1 # Softening size (internal units)
-
+
# Hernquist potential parameters
HernquistPotential:
- useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
+ useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
position: [100.,100.,100.] # Location of centre of Henrquist potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units)
idealizeddisk: 0 # (Optional) Whether to run with idealizeddisk or without, 0 used the mass and scalelength as mandatory parameters, while 1 uses more advanced disk dependent paramters
mass: 1e10 # (Optional 0) default parameter, Mass of the Hernquist potential
@@ -314,10 +318,10 @@ HernquistPotential:
bulgefraction: 0.00705852860979 # (Optional 1) Bulge mass fraction (equal to MB in MakeNewDisk and GalIC)
timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration
epsilon: 0.1 # Softening size (internal units)
-
+
# NFW potential parameters
NFWPotential:
- useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
+ useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
position: [0.0,0.0,0.0] # Location of centre of the NFW potential with respect to centre of the box (internal units) if useabspos=0 otherwise with respect to the 0,0,0, coordinates.
concentration: 8. # Concentration of the halo
M_200: 2.0e+12 # Mass of the halo (M_200 in internal units)
@@ -326,7 +330,7 @@ NFWPotential:
# NFW + Miyamoto-Nagai disk potential
NFW_MNPotential:
- useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
+ useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
position: [0.,0.,0.] # Location of centre of isothermal potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units)
timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration
epsilon: 0.01 # Softening size (internal units)
@@ -336,7 +340,7 @@ NFW_MNPotential:
Mdisk: 3.0 # Mass of the disk (internal units)
Rdisk: 4.0 # Effective radius of the disk (internal units)
Zdisk: 0.4704911 # Scale-height of the disk (internal units)
-
+
# Disk-patch potential parameters. The potential is along the x-axis.
DiscPatchPotential:
surface_density: 10. # Surface density of the disc (internal units)
@@ -423,7 +427,7 @@ COLIBRECooling:
He_reion_eV_p_H: 2.0 # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
rapid_cooling_threshold: 0.333333 # Switch to rapid cooling regime for dt / t_cool above this threshold.
delta_logTEOS_subgrid_properties: 0.3 # delta log T above the EOS below which the subgrid properties use Teq assumption
-
+
# Cooling with Grackle 3.0
GrackleCooling:
cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
@@ -458,7 +462,6 @@ EAGLEChemistry:
GEARChemistry:
initial_metallicity: 1 # Initial metallicity of the gas (mass fraction)
scale_initial_metallicity: 1 # Should we scale the initial metallicity with the solar one?
- diffusion_coefficient: 1e-3 # Coefficient for the diffusion (see Shen et al. 2010; differs by gamma^2 due to h^2).
# Parameters related to star formation models -----------------------------------------------
@@ -552,7 +555,7 @@ GEARFeedback:
# Parameters related to AGN models -----------------------------------------------
-
+
# EAGLE AGN model
EAGLEAGN:
subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
@@ -569,7 +572,7 @@ EAGLEAGN:
radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
- use_nibbling: 0 # Continuously transfer small amounts of mass from all gas neighbours to a black hole [1] or stochastically swallow whole gas particles [0]?
+ use_nibbling: 0 # Continuously transfer small amounts of mass from all gas neighbours to a black hole [1] or stochastically swallow whole gas particles [0]?
min_gas_mass_for_nibbling: 9e5 # Minimum mass for a gas particle to be nibbled from [M_Sun]. Only used if use_nibbling is 1.
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
AGN_feedback_model: Random # Feedback modes: Random, Isotropic, MinimumDistance, MinimumDensity
diff --git a/src/cosmology.c b/src/cosmology.c
index 038d87fb6b11d428cebd73cda89cd25efc8f9b61..bc0efcc5116e69ff0d8b891c601f84fb503c3bb8 100644
--- a/src/cosmology.c
+++ b/src/cosmology.c
@@ -43,7 +43,7 @@
#endif
/*! Number of values stored in the cosmological interpolation tables */
-const int cosmology_table_length = 10000;
+const int cosmology_table_length = 30000;
#ifdef HAVE_LIBGSL
/*! Size of the GSL workspace */
@@ -195,8 +195,11 @@ void cosmology_update(struct cosmology *c, const struct phys_const *phys_const,
/* Dark-energy equation of state */
c->w = cosmology_dark_energy_EoS(a, c->w_0, c->w_a);
+ /* Update the neutrino density */
+ c->Omega_nu = cosmology_get_neutrino_density(c, a);
+
/* E(z) */
- const double Omega_r = c->Omega_r;
+ const double Omega_r = c->Omega_r + c->Omega_nu;
const double Omega_m = c->Omega_m;
const double Omega_k = c->Omega_k;
const double Omega_l = c->Omega_lambda;
@@ -237,7 +240,8 @@ void cosmology_update(struct cosmology *c, const struct phys_const *phys_const,
double drift_integrand(double a, void *param) {
const struct cosmology *c = (const struct cosmology *)param;
- const double Omega_r = c->Omega_r;
+ const double Omega_nu = cosmology_get_neutrino_density(c, a);
+ const double Omega_r = c->Omega_r + Omega_nu;
const double Omega_m = c->Omega_m;
const double Omega_k = c->Omega_k;
const double Omega_l = c->Omega_lambda;
@@ -261,7 +265,8 @@ double drift_integrand(double a, void *param) {
double gravity_kick_integrand(double a, void *param) {
const struct cosmology *c = (const struct cosmology *)param;
- const double Omega_r = c->Omega_r;
+ const double Omega_nu = cosmology_get_neutrino_density(c, a);
+ const double Omega_r = c->Omega_r + Omega_nu;
const double Omega_m = c->Omega_m;
const double Omega_k = c->Omega_k;
const double Omega_l = c->Omega_lambda;
@@ -285,7 +290,8 @@ double gravity_kick_integrand(double a, void *param) {
double hydro_kick_integrand(double a, void *param) {
const struct cosmology *c = (const struct cosmology *)param;
- const double Omega_r = c->Omega_r;
+ const double Omega_nu = cosmology_get_neutrino_density(c, a);
+ const double Omega_r = c->Omega_r + Omega_nu;
const double Omega_m = c->Omega_m;
const double Omega_k = c->Omega_k;
const double Omega_l = c->Omega_lambda;
@@ -311,7 +317,8 @@ double hydro_kick_integrand(double a, void *param) {
double hydro_kick_corr_integrand(double a, void *param) {
const struct cosmology *c = (const struct cosmology *)param;
- const double Omega_r = c->Omega_r;
+ const double Omega_nu = cosmology_get_neutrino_density(c, a);
+ const double Omega_r = c->Omega_r + Omega_nu;
const double Omega_m = c->Omega_m;
const double Omega_k = c->Omega_k;
const double Omega_l = c->Omega_lambda;
@@ -334,7 +341,8 @@ double hydro_kick_corr_integrand(double a, void *param) {
double time_integrand(double a, void *param) {
const struct cosmology *c = (const struct cosmology *)param;
- const double Omega_r = c->Omega_r;
+ const double Omega_nu = cosmology_get_neutrino_density(c, a);
+ const double Omega_r = c->Omega_r + Omega_nu;
const double Omega_m = c->Omega_m;
const double Omega_k = c->Omega_k;
const double Omega_l = c->Omega_lambda;
@@ -349,6 +357,184 @@ double time_integrand(double a, void *param) {
return (1. / H) * a_inv;
}
+/**
+ * @brief Evaluates the neutrino density momentum integrand
+ * \f$ x^2 \sqrt{x^2 + y^2} / (1+e^x) \f$, where
+ * \f$ y = a*M_nu/(kb*T) \f$ is the neutrino mass scaled with redshift.
+ *
+ * This is used to evaluate the integral on (0, 1).
+ *
+ * @param x The momentum integration variable
+ * @param param Neutrino mass y scaled by temperature at redshift of interest.
+ * @return The integrand evaluated at x
+ */
+double neutrino_density_integrand(double x, void *param) {
+ double y = *(double *)param;
+ double numerator = x * x * hypot(x, y);
+
+ /* Handle overflows */
+ if (x > 20 + log(numerator)) {
+ return numerator * exp(-x);
+ }
+
+ return numerator / (1.0 + exp(x));
+}
+
+/**
+ * @brief Evaluates the transformed neutrino density momentum integrand
+ * \f$ w^{-4} \sqrt{w^{-2} + y^2} / (1+e^{-w}) \f$, where
+ * \f$ y = a*M_nu/(kb*T) \f$ is the neutrino mass scaled with redshift.
+ *
+ * This is used to evaluate the integral on (1, infinity).
+ *
+ * @param w The transformed momentum integration variable w=1/x
+ * @param param Neutrino mass y scaled by temperature at redshift of interest.
+ * @return The integrand evaluated at w
+ */
+double neutrino_density_integrand_transformed(double w, void *param) {
+ return neutrino_density_integrand(1. / w, param) / (w * w);
+}
+
+/**
+ * @brief Performs the neutrino density momentum integral
+ * \f$ \int_0^\infty x^2 \sqrt{x^2 + y^2} / (1+e^x) dx \f$, where
+ * \f$ y = a*M_nu/(kb*T) \f$ is the neutrino mass scaled with redshift,
+ * without pre-factors.
+ *
+ * @param space The GSL working space
+ * @param y Neutrino mass y scaled by temperature at redshift of interest.
+ * @return The integral evaluated at y
+ */
+double neutrino_density_integrate(gsl_integration_workspace *space, double y) {
+ double intermediate, abserr;
+
+ double result = 0;
+
+ gsl_function F1 = {&neutrino_density_integrand, &y};
+ gsl_function F2 = {&neutrino_density_integrand_transformed, &y};
+ /* Integrate between 0 and 1 */
+ gsl_integration_qag(&F1, 0.0, 1.0, 0, 1.0e-13, GSL_workspace_size,
+ GSL_INTEG_GAUSS61, space, &intermediate, &abserr);
+
+ result += intermediate;
+ /* Integrate between 1 and infinity */
+ gsl_integration_qag(&F2, 0.0, 1.0, 0, 1.0e-13, GSL_workspace_size,
+ GSL_INTEG_GAUSS61, space, &intermediate, &abserr);
+ result += intermediate;
+
+ return result;
+}
+
+/**
+ * @brief Find a time when all neutrinos are still relativistic. Store the
+ * starting, mid, and end points of the neutrino interpolation tables in c.
+ *
+ * @param c The cosmology structure
+ * @param tol Tolerance in density integral
+ * @param space The GSL working space
+ */
+void neutrino_find_relativistic_redshift(struct cosmology *c, double tol,
+ gsl_integration_workspace *space) {
+
+ /* Find the largest neutrino mass */
+ double M_max_eV = c->M_nu_eV[0];
+ for (int i = 1; i < c->N_nu; i++) {
+ M_max_eV = fmax(M_max_eV, c->M_nu_eV[i]);
+ }
+
+ /* A safe starting time when neutrinos are relativistic */
+ double a_safe = 0.5 * tol / M_max_eV;
+
+ /* Dont start the early table later than the simulation */
+ a_safe = fmin(a_safe, 0.9 * c->a_begin);
+
+ /* Start the late table just before the start of the simulation */
+ double a_midpoint = 0.99 * c->a_begin;
+
+ /* End the late table today (a=1) or at a_end, whichever is later */
+ double a_final = fmax(1.0, c->a_end);
+
+ /* Integrate early table on (a_start, a_mid) and late table on (a_mid, a_f) */
+ c->log_a_long_begin = log(a_safe);
+ c->log_a_long_mid = log(a_midpoint);
+ c->log_a_long_end = log(a_final);
+}
+
+/**
+ * @brief Initialise the neutrino density interpolation tables (early and late).
+ */
+void cosmology_init_neutrino_tables(struct cosmology *c) {
+
+ /* Skip if we have no massive neutrinos */
+ if (c->N_nu == 0) return;
+
+#ifdef HAVE_LIBGSL
+
+ /* Allocate memory for longer interpolation tables */
+ if (swift_memalign("cosmo.table", (void **)&c->neutrino_density_early_table,
+ SWIFT_STRUCT_ALIGNMENT,
+ cosmology_table_length * sizeof(double)) != 0)
+ error("Failed to allocate cosmology interpolation table");
+ if (swift_memalign("cosmo.table", (void **)&c->neutrino_density_late_table,
+ SWIFT_STRUCT_ALIGNMENT,
+ cosmology_table_length * sizeof(double)) != 0)
+ error("Failed to allocate cosmology interpolation table");
+
+ /* Initalise the GSL workspace */
+ gsl_integration_workspace *space =
+ gsl_integration_workspace_alloc(GSL_workspace_size);
+
+ /* Find a safe redshift to start the neutrino density interpolation table */
+ neutrino_find_relativistic_redshift(c, 1e-7, space);
+
+ const double pre_factor = 15. * pow(c->T_nu_0 * M_1_PI / c->T_CMB_0, 4);
+ const double early_delta_a =
+ (c->log_a_long_mid - c->log_a_long_begin) / cosmology_table_length;
+ const double late_delta_a =
+ (c->log_a_long_end - c->log_a_long_mid) / cosmology_table_length;
+
+ double result;
+
+ /* Fill the early neutrino density table between (a_long_begin, a_long_mid) */
+ for (int i = 0; i < cosmology_table_length; i++) {
+ double O_nu = 0.;
+ double a = exp(c->log_a_long_begin + early_delta_a * (i + 1));
+
+ /* Integrate the FD distribtution for each species */
+ for (int j = 0; j < c->N_nu; j++) {
+ double y = a * c->M_nu_eV[j] / c->T_nu_0_eV;
+ result = neutrino_density_integrate(space, y);
+ O_nu += c->deg_nu[j] * result * pre_factor * c->Omega_g;
+ }
+
+ c->neutrino_density_early_table[i] = O_nu;
+ }
+
+ /* Fill the late neutrino density table between (a_long_mid, a_long_end) */
+ for (int i = 0; i < cosmology_table_length; i++) {
+ double O_nu = 0.;
+ double a = exp(c->log_a_long_mid + late_delta_a * (i + 1));
+
+ /* Integrate the FD distribtution for each species */
+ for (int j = 0; j < c->N_nu; j++) {
+ double y = a * c->M_nu_eV[j] / c->T_nu_0_eV;
+ result = neutrino_density_integrate(space, y);
+ O_nu += c->deg_nu[j] * result * pre_factor * c->Omega_g;
+ }
+
+ c->neutrino_density_late_table[i] = O_nu;
+ }
+
+ /* Free the workspace */
+ gsl_integration_workspace_free(space);
+
+#else
+
+ error("Code not compiled with GSL. Can't compute cosmology integrals.");
+
+#endif
+}
+
/**
* @brief Initialise the interpolation tables for the integrals.
*/
@@ -525,11 +711,39 @@ void cosmology_init(struct swift_params *params, const struct unit_system *us,
c->Omega_r = parser_get_opt_param_double(params, "Cosmology:Omega_r", 0.);
c->Omega_lambda = parser_get_param_double(params, "Cosmology:Omega_lambda");
c->Omega_b = parser_get_param_double(params, "Cosmology:Omega_b");
- c->Omega_nu = 0.;
c->w_0 = parser_get_opt_param_double(params, "Cosmology:w_0", -1.);
c->w_a = parser_get_opt_param_double(params, "Cosmology:w_a", 0.);
c->h = parser_get_param_double(params, "Cosmology:h");
+ /* Neutrino temperature (inferred from T_CMB_0 if not specified) */
+ c->T_nu_0 = parser_get_opt_param_double(params, "Cosmology:T_nu_0", 0.);
+
+ /* Number of ultra-relativistic (massless) and massive neutrino species */
+ c->N_ur = parser_get_opt_param_double(params, "Cosmology:N_ur", 0.);
+ c->N_nu = parser_get_opt_param_int(params, "Cosmology:N_nu", 0);
+
+ /* Make sure that the cosmological parameters are not overdetermined */
+ if (c->Omega_r != 0. && c->N_ur != 0.) {
+ error("Cannot use both Cosmology:Omega_r and Cosmology:N_ur.");
+ }
+
+ /* If there are massive neutrinos, read the masses and degeneracies */
+ if (c->N_nu > 0) {
+ c->M_nu_eV = (double *)swift_malloc("Mnu", c->N_nu * sizeof(double));
+ c->deg_nu = (double *)swift_malloc("degnu", c->N_nu * sizeof(double));
+
+ /* Set default values */
+ for (int i = 0; i < c->N_nu; i++) {
+ c->M_nu_eV[i] = 0.0;
+ c->deg_nu[i] = 1.0;
+ }
+
+ parser_get_opt_param_double_array(params, "Cosmology:M_nu_eV", c->N_nu,
+ c->M_nu_eV);
+ parser_get_opt_param_double_array(params, "Cosmology:deg_nu", c->N_nu,
+ c->deg_nu);
+ }
+
/* Read the start and end of the simulation */
c->a_begin = parser_get_param_double(params, "Cosmology:a_begin");
c->a_end = parser_get_param_double(params, "Cosmology:a_end");
@@ -546,9 +760,6 @@ void cosmology_init(struct swift_params *params, const struct unit_system *us,
/* Construct derived quantities */
- /* Curvature density (for closure) */
- c->Omega_k = 1. - (c->Omega_m + c->Omega_r + c->Omega_lambda);
-
/* Dark-energy equation of state */
c->w = cosmology_dark_energy_EoS(c->a_begin, c->w_0, c->w_a);
@@ -563,14 +774,86 @@ void cosmology_init(struct swift_params *params, const struct unit_system *us,
c->critical_density_0 =
3. * c->H0 * c->H0 / (8. * M_PI * phys_const->const_newton_G);
- /* CMB temperature in internal units */
+ /* Some constants */
const double cc = phys_const->const_speed_light_c;
- const double T_CMB_0_4 = c->Omega_r * c->critical_density_0 * cc * cc * cc /
- (4. * phys_const->const_stefan_boltzmann);
- c->T_CMB_0 = pow(T_CMB_0_4, 1. / 4.);
+ const double rho_c3_on_4sigma = c->critical_density_0 * cc * cc * cc /
+ (4. * phys_const->const_stefan_boltzmann);
+
+ /* Handle neutrinos only if present */
+ if (c->N_ur == 0. && c->N_nu == 0) {
+ /* Infer T_CMB_0 from Omega_r */
+ c->T_CMB_0 = pow(c->Omega_r * rho_c3_on_4sigma, 1. / 4.);
+ c->T_CMB_0_K =
+ c->T_CMB_0 / units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
+
+ c->Omega_g = c->Omega_r;
+ c->T_nu_0 = 0.;
+ c->T_nu_0_eV = 0.;
+ c->Omega_ur = 0.;
+ c->Omega_nu_0 = 0.;
+ c->Omega_nu = 0.;
+ c->N_eff = 0.;
+
+ c->neutrino_density_early_table = NULL;
+ c->neutrino_density_late_table = NULL;
+ } else {
+ /* Infer T_CMB_0 from Omega_r if the latter is specified */
+ if (c->Omega_r != 0) {
+ c->T_CMB_0 = pow(c->Omega_r * rho_c3_on_4sigma, 1. / 4.);
+ } else {
+ c->T_CMB_0 = phys_const->const_T_CMB_0;
+ }
- c->T_CMB_0_K =
- c->T_CMB_0 / units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
+ /* Approximate the neutrino temperature if unspecified */
+ const double decoupling_factor = cbrt(4. / 11);
+ const double dec_4 = pow(decoupling_factor, 4);
+ if (c->T_nu_0 == 0.) {
+ c->T_nu_0 = c->T_CMB_0 * decoupling_factor;
+ }
+
+ /* Unit conversions */
+ c->T_CMB_0_K =
+ c->T_CMB_0 / units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
+ c->T_nu_0_eV = c->T_nu_0 * phys_const->const_boltzmann_k /
+ phys_const->const_electron_volt;
+
+ /* Set photon density and compute radiation density if necessary */
+ if (c->Omega_r != 0.) {
+ c->Omega_g = c->Omega_r;
+ c->Omega_ur = 0.;
+ } else {
+ /* Infer CMB density from the temperature */
+ c->Omega_g = pow(c->T_CMB_0, 4) / rho_c3_on_4sigma;
+
+ /* Compute the density of ultra-relativistic fermionic species */
+ c->Omega_ur = c->N_ur * (7. / 8.) * dec_4 * c->Omega_g;
+
+ /* Compute the total radiation density */
+ c->Omega_r = c->Omega_g + c->Omega_ur;
+ }
+
+ /* Compute effective number of relativistic species at early times */
+ double N_nu_tot_deg = 0.;
+ for (int i = 0; i < c->N_nu; i++) {
+ N_nu_tot_deg += c->deg_nu[i];
+ }
+ c->N_eff = c->N_ur + N_nu_tot_deg * pow(c->T_nu_0 / c->T_CMB_0, 4) / dec_4;
+
+ /* Initialise the neutrino density interpolation tables if necessary */
+ c->neutrino_density_early_table = NULL;
+ c->neutrino_density_late_table = NULL;
+ cosmology_init_neutrino_tables(c);
+
+ /* Retrieve the present-day total density due to massive neutrinos */
+ c->Omega_nu_0 = cosmology_get_neutrino_density(c, 1);
+ c->Omega_nu = c->Omega_nu_0; // will be updated
+ }
+
+ /* Cold dark matter density */
+ c->Omega_cdm = c->Omega_m - c->Omega_b;
+
+ /* Curvature density (for closure) */
+ c->Omega_k = 1. - (c->Omega_m + c->Omega_r + c->Omega_lambda + c->Omega_nu_0);
/* Initialise the interpolation tables */
c->drift_fac_interp_table = NULL;
@@ -580,7 +863,7 @@ void cosmology_init(struct swift_params *params, const struct unit_system *us,
c->time_interp_table_offset = 0.;
cosmology_init_tables(c);
- /* Set remaining variables to alid values */
+ /* Set remaining variables to valid values */
cosmology_update(c, phys_const, 0);
/* Update the times */
@@ -612,10 +895,22 @@ void cosmology_init_no_cosmo(struct cosmology *c) {
c->h = 1.;
c->w = -1.;
+ c->Omega_cdm = 0.;
+ c->Omega_ur = 0.;
+ c->Omega_g = 0.;
+ c->T_nu_0 = 0.;
+ c->T_nu_0_eV = 0.;
+ c->N_nu = 0;
+ c->N_ur = 0.;
+ c->N_eff = 0.;
+
c->a_begin = 1.;
c->a_end = 1.;
c->log_a_begin = 0.;
c->log_a_end = 0.;
+ c->log_a_long_begin = 0.;
+ c->log_a_long_mid = 0.;
+ c->log_a_long_end = 0.;
c->H = 0.;
c->H0 = 0.;
@@ -656,6 +951,8 @@ void cosmology_init_no_cosmo(struct cosmology *c) {
c->hydro_kick_fac_interp_table = NULL;
c->hydro_kick_corr_interp_table = NULL;
c->time_interp_table = NULL;
+ c->neutrino_density_early_table = NULL;
+ c->neutrino_density_late_table = NULL;
c->time_interp_table_offset = 0.;
c->scale_factor_interp_table = NULL;
@@ -836,6 +1133,30 @@ double cosmology_get_delta_time(const struct cosmology *c,
return t2 - t1;
}
+/**
+ * @brief Compute neutrino density parameter Omega_nu at the given scale-factor
+ * This is the effective present day value, i.e. must be multiplied by (1+z)^4
+ *
+ * @param c The current #cosmology.
+ * @param a The scale factor
+ * @return The density parameter
+ */
+double cosmology_get_neutrino_density(const struct cosmology *c, double a) {
+
+ if (c->N_nu == 0) return 0.;
+
+ const double log_a = log(a);
+
+ if (log_a < c->log_a_long_begin)
+ return c->neutrino_density_early_table[0];
+ else if (log_a < c->log_a_long_mid)
+ return interp_table(c->neutrino_density_early_table, log_a,
+ c->log_a_long_begin, c->log_a_long_mid);
+ else
+ return interp_table(c->neutrino_density_late_table, log_a,
+ c->log_a_long_mid, c->log_a_long_end);
+}
+
/**
* @brief Compute the cosmic time (in internal units) between two scale factors
*
@@ -935,15 +1256,31 @@ double cosmology_get_scale_factor(const struct cosmology *c, double t) {
void cosmology_print(const struct cosmology *c) {
message(
- "Density parameters: [O_m, O_l, O_b, O_k, O_r, O_nu] = [%f, %f, %f, %f, "
- "%f, %f]",
- c->Omega_m, c->Omega_lambda, c->Omega_b, c->Omega_k, c->Omega_r,
- c->Omega_nu);
+ "Density parameters: [O_m, O_l, O_b, O_k, O_r] = [%f, %f, %f, %f, %f]",
+ c->Omega_m, c->Omega_lambda, c->Omega_b, c->Omega_k, c->Omega_r);
+ message(
+ "Additional density parameters: [O_nu_0, O_cdm, O_ur, O_g] = [%f, "
+ "%f, %f, %f]",
+ c->Omega_nu_0, c->Omega_cdm, c->Omega_ur, c->Omega_g);
message("Dark energy equation of state: w_0=%f w_a=%f", c->w_0, c->w_a);
message("Hubble constant: h = %f, H_0 = %e U_t^(-1)", c->h, c->H0);
message("Hubble time: 1/H0 = %e U_t", c->Hubble_time);
- message("CMB temperature at z=0 derived from O_r: T_CMB = %e U_T",
+ message("CMB temperature at z=0 implied by cosmology: T_CMB = %e U_T",
c->T_CMB_0);
+ message("Neutrino temperature at z=0: T_nu = %e U_T", c->T_nu_0);
+ message("Numbers of relatistic species: [N_nu, N_ur, N_eff] = [%d, %f, %f]",
+ c->N_nu, c->N_ur, c->N_eff);
+ /* Print neutrino masses and degeneracies */
+ if (c->N_nu > 0) {
+ char neutrino_mass_string[10 * c->N_nu];
+ char neutrino_deg_string[10 * c->N_nu];
+ for (int i = 0; i < c->N_nu; i++) {
+ sprintf(neutrino_mass_string + i * 10, "%.2e ", c->M_nu_eV[i]);
+ sprintf(neutrino_deg_string + i * 10, "%.2e ", c->deg_nu[i]);
+ }
+ message("Neutrino masses: %seV", neutrino_mass_string);
+ message("Neutrino degeneracies: %s", neutrino_deg_string);
+ }
message("Universe age at present day: %e U_t",
c->universe_age_at_present_day);
}
@@ -956,6 +1293,12 @@ void cosmology_clean(struct cosmology *c) {
swift_free("cosmo.table", c->hydro_kick_corr_interp_table);
swift_free("cosmo.table", c->time_interp_table);
swift_free("cosmo.table", c->scale_factor_interp_table);
+ if (c->N_nu > 0) {
+ swift_free("cosmo.table", c->neutrino_density_early_table);
+ swift_free("cosmo.table", c->neutrino_density_late_table);
+ swift_free("Mnu", c->M_nu_eV);
+ swift_free("degnu", c->deg_nu);
+ }
}
#ifdef HAVE_HDF5
@@ -978,6 +1321,19 @@ void cosmology_write_model(hid_t h_grp, const struct cosmology *c) {
io_write_attribute_d(h_grp, "Omega_k", c->Omega_k);
io_write_attribute_d(h_grp, "Omega_lambda", c->Omega_lambda);
io_write_attribute_d(h_grp, "Omega_nu", c->Omega_nu);
+ io_write_attribute_d(h_grp, "Omega_nu_0", c->Omega_nu_0);
+ io_write_attribute_d(h_grp, "Omega_ur", c->Omega_ur);
+ io_write_attribute_d(h_grp, "Omega_cdm", c->Omega_cdm);
+ io_write_attribute_d(h_grp, "Omega_g", c->Omega_g);
+ io_write_attribute_d(h_grp, "T_nu_0 [internal units]", c->T_nu_0);
+ io_write_attribute_d(h_grp, "T_nu_0 [eV]", c->T_nu_0_eV);
+ io_write_attribute_d(h_grp, "N_eff", c->N_eff);
+ io_write_attribute_d(h_grp, "N_ur", c->N_ur);
+ io_write_attribute_d(h_grp, "N_nu", c->N_nu);
+ if (c->N_nu > 0) {
+ io_write_attribute(h_grp, "M_nu_eV", DOUBLE, c->M_nu_eV, c->N_nu);
+ io_write_attribute(h_grp, "deg_nu", DOUBLE, c->deg_nu, c->N_nu);
+ }
io_write_attribute_d(h_grp, "T_CMB_0 [internal units]", c->T_CMB_0);
io_write_attribute_d(h_grp, "T_CMB_0 [K]", c->T_CMB_0_K);
io_write_attribute_d(h_grp, "w_0", c->w_0);
@@ -999,6 +1355,16 @@ void cosmology_write_model(hid_t h_grp, const struct cosmology *c) {
void cosmology_struct_dump(const struct cosmology *cosmology, FILE *stream) {
restart_write_blocks((void *)cosmology, sizeof(struct cosmology), 1, stream,
"cosmology", "cosmology function");
+
+ /* Also store the neutrino mass and degeneracy arrays if necessary */
+ if (cosmology->N_nu > 0) {
+ restart_write_blocks((double *)cosmology->M_nu_eV, sizeof(double),
+ cosmology->N_nu, stream, "cosmology->M_nu_eV",
+ "neutrino masses eV");
+ restart_write_blocks((double *)cosmology->deg_nu, sizeof(double),
+ cosmology->N_nu, stream, "cosmology->deg_nu",
+ "neutrino degeneracies");
+ }
}
/**
@@ -1014,6 +1380,21 @@ void cosmology_struct_restore(int enabled, struct cosmology *cosmology,
restart_read_blocks((void *)cosmology, sizeof(struct cosmology), 1, stream,
NULL, "cosmology function");
+ /* Restore the neutrino mass and degeneracy arrays if necessary */
+ if (cosmology->N_nu > 0) {
+ cosmology->M_nu_eV =
+ (double *)swift_malloc("Mnu", cosmology->N_nu * sizeof(double));
+ restart_read_blocks((double *)cosmology->M_nu_eV, sizeof(double),
+ cosmology->N_nu, stream, NULL, "neutrino masses eV");
+ cosmology->deg_nu =
+ (double *)swift_malloc("degnu", cosmology->N_nu * sizeof(double));
+ restart_read_blocks((double *)cosmology->deg_nu, sizeof(double),
+ cosmology->N_nu, stream, NULL, "neutrino degeneracies");
+ }
+
/* Re-initialise the tables if using a cosmology. */
- if (enabled) cosmology_init_tables(cosmology);
+ if (enabled) {
+ cosmology_init_neutrino_tables(cosmology);
+ cosmology_init_tables(cosmology);
+ }
}
diff --git a/src/cosmology.h b/src/cosmology.h
index 5a33229a9d7255e47eb79772c7b59202e59e20c8..a72474d16c37403bcc5c25aa127027c63e853a3f 100644
--- a/src/cosmology.h
+++ b/src/cosmology.h
@@ -144,12 +144,24 @@ struct cosmology {
/*! Cosmological constant density parameter */
double Omega_lambda;
- /*! Radiation constant density parameter */
+ /*! Total radiation density parameter (photons and other relics) */
double Omega_r;
- /*! Neutrino density parameter */
+ /*! CMB radiation density parameter (Omega_gamma) */
+ double Omega_g;
+
+ /*! Massive neutrino density parameter */
double Omega_nu;
+ /*! Massive neutrino density parameter at z=0 */
+ double Omega_nu_0;
+
+ /*! Ultra-relativistic species (e.g. massless neutrinos) density parameter */
+ double Omega_ur;
+
+ /*! Cold dark matter density parameter */
+ double Omega_cdm;
+
/*! Curvature density parameter */
double Omega_k;
@@ -159,12 +171,42 @@ struct cosmology {
/*! Dark-energy evolution parameter */
double w_a;
- /*! CMB temperature at z = 0 derived from Omega_r (internal units) */
+ /*! CMB temperature at z = 0 implied by cosmology (internal units) */
double T_CMB_0;
- /*! CMB temperature at z = 0 derived from Omega_r (Kelvins) */
+ /*! CMB temperature at z = 0 implied by cosmology (Kelvins) */
double T_CMB_0_K;
+ /*! Neutrino temperature at z = 0 (internal units) */
+ double T_nu_0;
+
+ /* Neutrino temperature at z = 0 (electron-volts) */
+ double T_nu_0_eV;
+
+ /*! Number of massive neutrino species */
+ int N_nu;
+
+ /*! Number of ultra-relativistic species (excluding massive neutrinos) */
+ double N_ur;
+
+ /*! Effective number of relativistic species (including massive neutrinos) */
+ double N_eff;
+
+ /*! Mass of each massive neutrino species in electron-volts */
+ double *M_nu_eV;
+
+ /*! Degeneracy of each massive neutrino species */
+ double *deg_nu;
+
+ /*! Log of starting expansion factor for neutrino interpolation tables */
+ double log_a_long_begin;
+
+ /*! Log of midpoint expansion factor for neutrino interpolation tables */
+ double log_a_long_mid;
+
+ /*! Log of ending expansion factor for neutrino interpolation tables */
+ double log_a_long_end;
+
/*! Log of starting expansion factor */
double log_a_begin;
@@ -189,6 +231,12 @@ struct cosmology {
/*! Scale factor interpolation table */
double *scale_factor_interp_table;
+ /*! Massive neutrino density interpolation table at early times */
+ double *neutrino_density_early_table;
+
+ /*! Massive neutrino density interpolation table at late times */
+ double *neutrino_density_late_table;
+
/*! Time between Big Bang and first entry in the table */
double time_interp_table_offset;
@@ -217,6 +265,7 @@ double cosmology_get_corr_kick_factor(const struct cosmology *cosmo,
double cosmology_get_delta_time(const struct cosmology *c,
const integertime_t ti_start,
const integertime_t ti_end);
+double cosmology_get_neutrino_density(const struct cosmology *c, double a);
double cosmology_get_delta_time_from_scale_factors(const struct cosmology *c,
const double a_start,
diff --git a/src/engine.c b/src/engine.c
index f2e5c3d4c6629a662488012dcbfa742f38c7b104..2a0492ee2c55abd0d572b9cce100dd3a84c20ac7 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2932,7 +2932,7 @@ void engine_recompute_displacement_constraint(struct engine *e) {
/* Get the cosmological information */
const int with_cosmology = e->policy & engine_policy_cosmology;
const struct cosmology *cosmo = e->cosmology;
- const float Om = cosmo->Omega_m;
+ const float Ocdm = cosmo->Omega_cdm;
const float Ob = cosmo->Omega_b;
const float H0 = cosmo->H0;
const float a = cosmo->a;
@@ -3015,7 +3015,7 @@ void engine_recompute_displacement_constraint(struct engine *e) {
const float min_mass_dm = min_mass[1];
/* Inter-particle sepration for the DM */
- const float d_dm = cbrtf(min_mass_dm / ((Om - Ob) * rho_crit0));
+ const float d_dm = cbrtf(min_mass_dm / (Ocdm * rho_crit0));
/* RMS peculiar motion for the DM */
const float rms_vel_dm = vel_norm_dm / N_dm;
diff --git a/src/space.c b/src/space.c
index 52fbdeb96b9b861c1de23431397258bdddfa95c9..1b60d924742df427f57a9df8350deae953b41e54 100644
--- a/src/space.c
+++ b/src/space.c
@@ -1019,7 +1019,7 @@ void space_init(struct space *s, struct swift_params *params,
size_t Nspart, size_t Nbpart, int periodic, int replicate,
int remap_ids, int generate_gas_in_ics, int hydro,
int self_gravity, int star_formation, int DM_background,
- int verbose, int dry_run, int nr_nodes) {
+ int neutrinos, int verbose, int dry_run, int nr_nodes) {
/* Clean-up everything */
bzero(s, sizeof(struct space));
@@ -1033,6 +1033,7 @@ void space_init(struct space *s, struct swift_params *params,
s->with_hydro = hydro;
s->with_star_formation = star_formation;
s->with_DM_background = DM_background;
+ s->with_neutrinos = neutrinos;
s->nr_parts = Npart;
s->nr_gparts = Ngpart;
s->nr_sparts = Nspart;
@@ -1855,13 +1856,21 @@ void space_check_cosmology(struct space *s, const struct cosmology *cosmo,
const double rho_crit0 = cosmo->critical_density * H0 * H0 / (H * H);
/* Compute the mass density */
- const double Omega_m = (total_mass / volume) / rho_crit0;
+ const double Omega_sim = (total_mass / volume) / rho_crit0;
- if (fabs(Omega_m - cosmo->Omega_m) > 1e-3)
+ /* The density required to match the cosmology */
+ double Omega_cosmo;
+ if (s->with_neutrinos)
+ Omega_cosmo = cosmo->Omega_m + cosmo->Omega_nu_0;
+ else
+ Omega_cosmo = cosmo->Omega_m;
+
+ if (fabs(Omega_sim - Omega_cosmo) > 1e-3)
error(
"The matter content of the simulation does not match the cosmology "
- "in the parameter file cosmo.Omega_m=%e Omega_m=%e",
- cosmo->Omega_m, Omega_m);
+ "in the parameter file cosmo.Omega_m=%e Omega_particles=%e. Are you"
+ "running with neutrinos? Then account for cosmo.Omega_nu_0=%e too.",
+ cosmo->Omega_m, Omega_sim, cosmo->Omega_nu_0);
}
}
diff --git a/src/space.h b/src/space.h
index c726a68d9873f7473c191f649fb3d7fa53b10400..47a8ffe153c3bc4da9101c71a9aa4fbde32bd6fd 100644
--- a/src/space.h
+++ b/src/space.h
@@ -111,6 +111,9 @@ struct space {
/*! Are we running with some DM background particles? */
int with_DM_background;
+ /*! Are we running with neutrino particles? */
+ int with_neutrinos;
+
/*! Width of the top-level cells. */
double width[3];
@@ -346,8 +349,8 @@ void space_init(struct space *s, struct swift_params *params,
struct bpart *bparts, size_t Npart, size_t Ngpart, size_t Nsink,
size_t Nspart, size_t Nbpart, int periodic, int replicate,
int remap_ids, int generate_gas_in_ics, int hydro, int gravity,
- int star_formation, int DM_background, int verbose, int dry_run,
- int nr_nodes);
+ int star_formation, int DM_background, int neutrinos,
+ int verbose, int dry_run, int nr_nodes);
void space_sanitize(struct space *s);
void space_map_cells_pre(struct space *s, int full,
void (*fun)(struct cell *c, void *data), void *data);
diff --git a/src/velociraptor_interface.c b/src/velociraptor_interface.c
index 231349e0d7a10d8920926e876fc09337e8633d63..aabb58906e626f03f4ede727b7462e261d18ea19 100644
--- a/src/velociraptor_interface.c
+++ b/src/velociraptor_interface.c
@@ -66,7 +66,7 @@ struct cosmoinfo {
/*! Radiation density parameter (cosmology.Omega_r) */
double Omega_r;
- /*! Neutrino density parameter (0 in SWIFT) */
+ /*! Neutrino density parameter at z = 0 (cosmology.Omega_nu_0) */
double Omega_nu;
/*! Neutrino density parameter (cosmology.Omega_k) */
@@ -78,7 +78,7 @@ struct cosmoinfo {
/*! Radiation constant density parameter (cosmology.Omega_lambda) */
double Omega_Lambda;
- /*! Dark matter density parameter (cosmology.Omega_m - cosmology.Omega_b) */
+ /*! Dark matter density parameter (cosmology.Omega_cdm) */
double Omega_cdm;
/*! Dark-energy equation of state at the current time (cosmology.w)*/
@@ -774,9 +774,9 @@ void velociraptor_invoke(struct engine *e, const int linked_with_snap) {
cosmo_info.Omega_b = e->cosmology->Omega_b;
cosmo_info.Omega_r = e->cosmology->Omega_r;
cosmo_info.Omega_k = e->cosmology->Omega_k;
- cosmo_info.Omega_nu = 0.;
+ cosmo_info.Omega_nu = e->cosmology->Omega_nu_0;
cosmo_info.Omega_Lambda = e->cosmology->Omega_lambda;
- cosmo_info.Omega_cdm = e->cosmology->Omega_m - e->cosmology->Omega_b;
+ cosmo_info.Omega_cdm = e->cosmology->Omega_cdm;
cosmo_info.w_de = e->cosmology->w;
/* Report the cosmo info we use */
@@ -787,6 +787,7 @@ void velociraptor_invoke(struct engine *e, const int linked_with_snap) {
message("VELOCIraptor conf: Omega_b: %e", cosmo_info.Omega_b);
message("VELOCIraptor conf: Omega_Lambda: %e", cosmo_info.Omega_Lambda);
message("VELOCIraptor conf: Omega_cdm: %e", cosmo_info.Omega_cdm);
+ message("VELOCIraptor conf: Omega_nu: %e", cosmo_info.Omega_nu);
message("VELOCIraptor conf: w_de: %e", cosmo_info.w_de);
}
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 347de263b05836f2579739f3666833736b64968f..c1b3e975705d65a1de4f4a58e0f726e3df2b7522 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -30,7 +30,7 @@ TESTS = testGreetings testMaths testReading.sh testKernel testKernelLongGrav \
testPotentialPair testEOS testUtilities testSelectOutput.sh \
testCbrt testCosmology testOutputList testFormat.sh \
test27cellsStars.sh test27cellsStarsPerturbed.sh testHydroMPIrules \
- testAtomic testGravitySpeed
+ testAtomic testGravitySpeed testNeutrinoCosmology
# List of test programs to compile
check_PROGRAMS = testGreetings testReading testTimeIntegration testKernelLongGrav \
@@ -43,7 +43,7 @@ check_PROGRAMS = testGreetings testReading testTimeIntegration testKernelLongGra
testGravityDerivatives testPotentialSelf testPotentialPair testEOS testUtilities \
testSelectOutput testCbrt testCosmology testOutputList test27cellsStars \
test27cellsStars_subset testCooling testComovingCooling testFeedback testHashmap \
- testAtomic testHydroMPIrules testGravitySpeed
+ testAtomic testHydroMPIrules testGravitySpeed testNeutrinoCosmology
# Rebuild tests when SWIFT is updated.
$(check_PROGRAMS): ../src/.libs/libswiftsim.a
diff --git a/tests/testFFT.c b/tests/testFFT.c
index 1654120b6ca6f059d2be07834ef336be877a86a3..415802d0da276bf7a5b1a7adf219d2d54702796e 100644
--- a/tests/testFFT.c
+++ b/tests/testFFT.c
@@ -84,7 +84,7 @@ int main(int argc, char *argv[]) {
struct space space;
double dim[3] = {1., 1., 1.};
space_init(&space, params, &cosmo, dim, NULL, gparts, NULL, 0, nr_gparts, 0,
- 1, 1, 0, 1, 1, 0);
+ 1, 1, 0, 0, 1, 1, 0);
struct engine engine;
engine.s = &space;
diff --git a/tests/testNeutrinoCosmology.c b/tests/testNeutrinoCosmology.c
new file mode 100644
index 0000000000000000000000000000000000000000..38b2db49bad9a60919630e94c8c0f4d5ef72a7f4
--- /dev/null
+++ b/tests/testNeutrinoCosmology.c
@@ -0,0 +1,264 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2019 Willem Elbers (whe@willemelbers.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/* Some standard headers. */
+#include "../config.h"
+
+/* Some standard headers */
+#include <fenv.h>
+#include <math.h>
+
+/* Includes. */
+#include "swift.h"
+
+#define N_CHECK 20
+#define TOLERANCE 1e-6
+#define EXTERNAL_TOLERANCE 1e-5 // comparing against CLASS
+
+void test_params_init(struct swift_params *params, int testnr) {
+ switch (testnr) {
+ /* One doubly degenerate light neutrino (total 0.10 eV) and one massless */
+ case 0:
+ parser_init("", params);
+ parser_set_param(params, "Cosmology:Omega_m:0.3122688651165");
+ parser_set_param(params, "Cosmology:Omega_lambda:6.853646e-01");
+ parser_set_param(params, "Cosmology:Omega_b:0.049198926683");
+ parser_set_param(params, "Cosmology:h:0.6737");
+ parser_set_param(params, "Cosmology:a_begin:1e-2");
+ parser_set_param(params, "Cosmology:a_end:1.0");
+ parser_set_param(params, "Cosmology:T_nu_0:1.95176");
+ parser_set_param(params, "Cosmology:N_nu:1");
+ parser_set_param(params, "Cosmology:N_ur:1.0196");
+ parser_set_param(params, "Cosmology:M_nu_eV:0.0486");
+ parser_set_param(params, "Cosmology:deg_nu:2");
+ break;
+ /* Three very massive neutrinos */
+ case 1:
+ parser_init("", params);
+ parser_set_param(params, "Cosmology:Omega_m:0.291188231748");
+ parser_set_param(params, "Cosmology:Omega_lambda:0.685365");
+ parser_set_param(params, "Cosmology:Omega_b:0.049199");
+ parser_set_param(params, "Cosmology:h:0.6737");
+ parser_set_param(params, "Cosmology:a_begin:1e-2");
+ parser_set_param(params, "Cosmology:a_end:1.0");
+ parser_set_param(params, "Cosmology:N_nu:3");
+ parser_set_param(params, "Cosmology:M_nu_eV:0.2,0.3,0.4");
+ break;
+ /* One light massive neutrino (0.05 eV) + two massless */
+ case 2:
+ parser_init("", params);
+ parser_set_param(params, "Cosmology:Omega_m:0.3110861814103");
+ parser_set_param(params, "Cosmology:Omega_lambda:0.685365");
+ parser_set_param(params, "Cosmology:Omega_b:0.049199");
+ parser_set_param(params, "Cosmology:h:0.6737");
+ parser_set_param(params, "Cosmology:a_begin:1e-2");
+ parser_set_param(params, "Cosmology:a_end:1.0");
+ parser_set_param(params, "Cosmology:T_nu_0:1.95176");
+ parser_set_param(params, "Cosmology:N_ur:2");
+ parser_set_param(params, "Cosmology:N_nu:1");
+ parser_set_param(params, "Cosmology:M_nu_eV:0.05");
+ break;
+ /* Three very massive neutrinos (same as 1), but starting earlier */
+ case 3:
+ parser_init("", params);
+ parser_set_param(params, "Cosmology:Omega_m:0.291188231748");
+ parser_set_param(params, "Cosmology:Omega_lambda:0.685365");
+ parser_set_param(params, "Cosmology:Omega_b:0.049199");
+ parser_set_param(params, "Cosmology:h:0.6737");
+ parser_set_param(params, "Cosmology:a_begin:1e-3");
+ parser_set_param(params, "Cosmology:a_end:1.0");
+ parser_set_param(params, "Cosmology:N_nu:3");
+ parser_set_param(params, "Cosmology:M_nu_eV:0.2,0.3,0.4");
+ break;
+ }
+
+ parser_set_param(params, "InternalUnitSystem:UnitMass_in_cgs:1.988435e+43");
+ parser_set_param(params,
+ "InternalUnitSystem:UnitLength_in_cgs:3.0856775815e+24");
+ parser_set_param(params,
+ "InternalUnitSystem:UnitVelocity_in_cgs:9.7846194238e+07");
+ parser_set_param(params, "InternalUnitSystem:UnitCurrent_in_cgs:1");
+ parser_set_param(params, "InternalUnitSystem:UnitTemp_in_cgs:1");
+}
+
+int main(int argc, char *argv[]) {
+
+ /* Load a CLASS data table of background quantities for cosmology 0 */
+ const int rows = 51;
+ const int cols = 3;
+ double CLASS_table[rows * cols];
+ FILE *stream = fopen("testNeutrinoCosmology.dat", "r");
+ char line[1024];
+ int row = 0;
+ while (fgets(line, 1024, stream)) {
+ if (line[0] == '#') continue;
+ char *tmp = strdup(line);
+ char *ptr = NULL;
+ CLASS_table[0 + row * 3] = strtod(tmp, &ptr);
+ CLASS_table[1 + row * 3] = strtod(ptr + 1, &ptr);
+ CLASS_table[2 + row * 3] = strtod(ptr + 1, &ptr);
+ row++;
+ free(tmp);
+ }
+ fclose(stream);
+
+ /* Initialize CPU frequency, this also starts time. */
+ unsigned long long cpufreq = 0;
+ clocks_set_cpufreq(cpufreq);
+
+ /* Test a number of different cosmologies */
+ int N_cosmos = 4;
+ double times1[N_CHECK]; // compare t(a) for cosmology 1
+ double times2[N_CHECK]; // with t(a) for cosmology 3
+ for (int testnr = 0; testnr < N_cosmos; testnr++) {
+ message("Initialization of cosmology %i", testnr);
+
+ /* pseudo initialization of params */
+ struct swift_params params;
+ test_params_init(¶ms, testnr);
+
+ /* initialization of unit system */
+ struct unit_system us;
+ units_init_from_params(&us, ¶ms, "InternalUnitSystem");
+
+ /* initialization of phys_const */
+ struct phys_const phys_const;
+ phys_const_init(&us, ¶ms, &phys_const);
+
+ /* initialization of cosmo */
+ struct cosmology cosmo;
+ cosmology_init(¶ms, &us, &phys_const, &cosmo);
+
+ message("Start checking computation...");
+
+ for (int i = 0; i < N_CHECK; i++) {
+ double a = 0.1 + 0.9 * i / (N_CHECK - 1.);
+ /* Compute a(t(a)) and check if same results */
+ double time = cosmology_get_time_since_big_bang(&cosmo, a);
+
+ /* Store the value for cosmologies 1 and 3 to compare later */
+ if (testnr == 1) {
+ times1[i] = time;
+ } else if (testnr == 3) {
+ times2[i] = time;
+ }
+
+ double my_a = cosmology_get_scale_factor(&cosmo, time);
+
+ /* check accuracy */
+ double rel_err = (my_a - a) / a;
+ message("Accuracy of %g at a=%g", rel_err, a);
+ assert(fabs(rel_err) < TOLERANCE);
+ }
+
+ if (testnr == 0) {
+ /* For cosmology 0, compare against independently computed array */
+ message("Relative to CLASS values (a, t, Omega_nu(a))");
+ const double delta_a = (cosmo.log_a_end - cosmo.log_a_begin) / 10000;
+ for (int j = 0; j < 50; j++) {
+ double a1 = exp(cosmo.log_a_begin + delta_a * (j * 200 + 1));
+ double time1 = cosmology_get_time_since_big_bang(&cosmo, a1);
+ // double time1 = cosmo.time_interp_table[j*200] +
+ // cosmo.time_interp_table_offset;
+ double Onu1 = cosmology_get_neutrino_density(&cosmo, a1);
+
+ double a2 = CLASS_table[0 + 3 * j];
+ double time2 = CLASS_table[1 + 3 * j];
+ double Onu2 = CLASS_table[2 + 3 * j];
+
+ /* Class defines years as 365.25 days, we use 365 days. */
+ time2 *= 365.25 / 365.0;
+
+ assert(fabs(a1 - a2) / (a1 + a2) < EXTERNAL_TOLERANCE);
+ assert(fabs(time1 - time2) / (time1 + time2) < EXTERNAL_TOLERANCE);
+ assert(fabs(Onu1 - Onu2) / (Onu1 + Onu2) < EXTERNAL_TOLERANCE);
+
+ message("At a = %e: (%e, %e, %e)", a1, a1 / a2, time1 / time2,
+ Onu1 / Onu2);
+ }
+ }
+
+ message("Everything seems fine with this cosmology.");
+
+ cosmology_clean(&cosmo);
+ }
+
+ /* Compare cosmologies 0 and 3 */
+ for (int j = 0; j < N_CHECK; j++) {
+ /* check accuracy */
+ double err = (times1[j] - times2[j]) / times1[j];
+ message("Agreement of %g at step %i", err, j);
+ assert(fabs(err) < TOLERANCE);
+ }
+
+ /* For the massive neutrino cosmology 3, check that it reproduces
+ * the relativistic (early) and non-relativistic (late) limits */
+
+ message("Initialization...");
+
+ /* pseudo initialization of params */
+ struct swift_params params;
+ test_params_init(¶ms, 3);
+
+ /* initialization of unit system */
+ struct unit_system us;
+ units_init_cgs(&us);
+
+ /* initialization of phys_const */
+ struct phys_const phys_const;
+ phys_const_init(&us, ¶ms, &phys_const);
+
+ /* initialization of cosmo */
+ struct cosmology cosmo;
+ cosmology_init(¶ms, &us, &phys_const, &cosmo);
+
+ message("Start checking the fermion integration...");
+
+ /* Relativistic limit */
+ double Omega_nu_early = cosmology_get_neutrino_density(&cosmo, 1e-10);
+ double Omega_nu_r =
+ cosmo.Omega_g * cosmo.N_eff * (7. / 8.) * pow(4. / 11., 4. / 3.);
+ double err = (Omega_nu_early - Omega_nu_r) / Omega_nu_early;
+ message("Accuracy of %g of the relativistic limit", err);
+ assert(fabs(err) < TOLERANCE);
+
+ /* Zeta function and other constants */
+ const double zeta3 = 1.202056903159594;
+ const double kb = phys_const.const_boltzmann_k;
+ const double hbar = phys_const.const_planck_h / (2 * M_PI);
+ const double cvel = phys_const.const_speed_light_c;
+ const double eV = phys_const.const_electron_volt;
+
+ /* Non-relativistic limit (limit not reached, so lower tolerance is fine)*/
+ double Omega_nu_late = cosmology_get_neutrino_density(&cosmo, 1.0);
+ double Omega_nu_nr = 0.;
+ for (int i = 0; i < cosmo.N_nu; i++) {
+ Omega_nu_nr += 6 * zeta3 / (11 * M_PI * M_PI) * pow(kb * cosmo.T_CMB_0, 3) /
+ pow(cvel * hbar, 3) * cosmo.M_nu_eV[i] * eV /
+ (cosmo.critical_density_0 * cvel * cvel);
+ }
+ double err2 = (Omega_nu_late - Omega_nu_nr) / Omega_nu_late;
+ message("Accuracy of %g of the non-relativistic limit", err2);
+ assert(fabs(err2) < 1e-5);
+
+ message("Neutrino cosmology test successful.");
+
+ cosmology_clean(&cosmo);
+ return 0;
+}
diff --git a/tests/testNeutrinoCosmology.dat b/tests/testNeutrinoCosmology.dat
new file mode 100644
index 0000000000000000000000000000000000000000..351cd2756a246da2ffcd44ccd68917b1eaeef39a
--- /dev/null
+++ b/tests/testNeutrinoCosmology.dat
@@ -0,0 +1,53 @@
+#Computed with CLASS v. 2.9
+#a,proper_time_Gpc,O_nu
+0.01000461,0.016701397553317,3.5163787160972E-05
+0.01096983,0.019228996526818,3.67516040528256E-05
+0.01202818,0.02213429919907,3.85569090708168E-05
+0.01318864,0.025473378968706,4.06036826343241E-05
+0.01446106,0.0293106271508,4.29178039616624E-05
+0.01585623,0.033720161738937,4.55272171628141E-05
+0.01738601,0.038786778537091,4.84619440279236E-05
+0.01906338,0.044607929499909,5.17543596186455E-05
+0.02090259,0.051295577208859,5.54394196800272E-05
+0.02291923,0.058978141111643,5.95548023375683E-05
+0.02513043,0.067803393711092,6.41413614902226E-05
+0.02755497,0.077940926147672,6.92433337546202E-05
+0.03021343,0.08958501664697,7.49086656511501E-05
+0.03312836,0.102958902483026,8.11895180833707E-05
+0.03632453,0.118319082448388,8.8142821279131E-05
+0.03982905,0.135959789513129,9.58306007331672E-05
+0.04367169,0.156219964076042,0.000104320794266
+0.04788506,0.17948690789494,0.000113687563152
+0.05250492,0.206205809197034,0.000124012192172
+0.0575705,0.236887879351767,0.00013538382664
+0.0631248,0.272119744284891,0.000147900195733
+0.06921496,0.312575279195446,0.000161668510521
+0.0758927,0.359027968989762,0.000176806479732
+0.08321469,0.412362655983059,0.000193443160009
+0.09124309,0.473595514035587,0.000211720276868
+0.1000461,0.543891750068211,0.00023179347611
+0.1096983,0.62458472301863,0.000253833120516
+0.1202818,0.717207770565306,0.000278026953631
+0.1318864,0.823506701592335,0.00030458041981
+0.1446106,0.945478283370435,0.000333719041739
+0.1585623,1.0854002866677,0.000365690020889
+0.1738601,1.24586901225677,0.000400765218598
+0.1906338,1.42983264427013,0.000439242312846
+0.2090259,1.64063166511213,0.000481448238577
+0.2291923,1.88201267212233,0.00052774077823
+0.2513043,2.15817979552822,0.000578513184864
+0.2755497,2.47380033173388,0.000634196697354
+0.3021343,2.83399860301076,0.000695263986152
+0.3312836,3.24433010205929,0.000762233064932
+0.3632453,3.71068357665048,0.000835673017592
+0.3982905,4.23914128398719,0.000916206768316
+0.4367169,4.83578207787208,0.00100451843668
+0.4788506,5.50636437013827,0.001101357471308
+0.5250492,6.25594667098269,0.001207545774545
+0.575705,7.08845062088615,0.001323985061712
+0.631248,8.00618086446026,0.001451663744548
+0.6921496,9.00947265502665,0.001591665447889
+0.758927,10.0964700596076,0.001745179121787
+0.8321469,11.2631241396981,0.001913507546324
+0.9124309,12.503470676376,0.002098079745641
+1,13.80342381474,0.002299403771466
diff --git a/theory/Cosmology/bibliography.bib b/theory/Cosmology/bibliography.bib
index 84cec263d2e8195bc831e672184b41d61479fcc2..7467e64a99e1a497f3c550f67b47e0f4a41a37e3 100644
--- a/theory/Cosmology/bibliography.bib
+++ b/theory/Cosmology/bibliography.bib
@@ -176,3 +176,61 @@ archivePrefix = "arXiv",
adsurl = {http://adsabs.harvard.edu/abs/2010MNRAS.401..791S},
adsnote = {Provided by the SAO/NASA Astrophysics Data System}
}
+
+@article{Zennaro2016,
+ title={Initial conditions for accurate N-body simulations of massive neutrino cosmologies},
+ author={Zennaro, Matteo and Bel, Julien and Villaescusa-Navarro, Francisco and Carbone, Carmelita and Sefusatti, Emiliano and Guzzo, Luigi},
+ journal={Monthly Notices of the Royal Astronomical Society},
+ volume={466},
+ number={3},
+ pages={3244--3258},
+ year={2016},
+ publisher={Oxford University Press}
+}
+
+@article{Elbers2020,
+ author = "Elbers, Willem and Frenk, Carlos S. and Jenkins, Adrian and Li, Baojiu and Pascoli, Silvia",
+ title = "{An optimal nonlinear method for simulating relic neutrinos}",
+ eprint = "2010.07321",
+ archivePrefix = "arXiv",
+ primaryClass = "astro-ph.CO",
+ month = "10",
+ year = "2020"
+}
+
+@ARTICLE{Mangano2005,
+ author = {{Mangano}, Gianpiero and {Miele}, Gennaro and {Pastor}, Sergio and {Pinto}, Teguayco and {Pisanti}, Ofelia and {Serpico}, Pasquale D.},
+ title = "{Relic neutrino decoupling including flavour oscillations}",
+ journal = {Nuclear Physics B},
+ keywords = {High Energy Physics - Phenomenology, Astrophysics},
+ year = 2005,
+ month = nov,
+ volume = {729},
+ number = {1-2},
+ pages = {221-234},
+ doi = {10.1016/j.nuclphysb.2005.09.041},
+archivePrefix = {arXiv},
+ eprint = {hep-ph/0506164},
+ primaryClass = {hep-ph},
+ adsurl = {https://ui.adsabs.harvard.edu/abs/2005NuPhB.729..221M},
+ adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+@ARTICLE{Lesgourgues2011,
+ author = {{Lesgourgues}, Julien and {Tram}, Thomas},
+ title = "{The Cosmic Linear Anisotropy Solving System (CLASS) IV: efficient implementation of non-cold relics}",
+ journal = {\jcap},
+ keywords = {Astrophysics - Cosmology and Nongalactic Astrophysics},
+ year = 2011,
+ month = sep,
+ volume = {2011},
+ number = {9},
+ eid = {032},
+ pages = {032},
+ doi = {10.1088/1475-7516/2011/09/032},
+archivePrefix = {arXiv},
+ eprint = {1104.2935},
+ primaryClass = {astro-ph.CO},
+ adsurl = {https://ui.adsabs.harvard.edu/abs/2011JCAP...09..032L},
+ adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
diff --git a/theory/Cosmology/cosmology_standalone.tex b/theory/Cosmology/cosmology_standalone.tex
index 5b5fa228fe4cd1c5cfbd64a5ddb7a7ec466fa7a7..72b802761fad0c87fc12f5895a4f6c887af294bd 100644
--- a/theory/Cosmology/cosmology_standalone.tex
+++ b/theory/Cosmology/cosmology_standalone.tex
@@ -37,6 +37,8 @@ Making cosmology great again.
\input{flrw}
+\input{neutrinos}
+
\input{coordinates}
\input{operators}
diff --git a/theory/Cosmology/flrw.tex b/theory/Cosmology/flrw.tex
index 875d4a4cc3eb8057f62003f90d30432ad879d4dc..d6e8cfcda6ee5eebf53d623afdae00118c47fa2a 100644
--- a/theory/Cosmology/flrw.tex
+++ b/theory/Cosmology/flrw.tex
@@ -19,7 +19,7 @@ and write
\begin{align}
H(a) &\equiv H_0 E(a) \\ E(a) &\equiv\sqrt{\Omega_m a^{-3} + \Omega_r
a^{-4} + \Omega_k a^{-2} + \Omega_\Lambda \exp\left(3\tilde{w}(a)\right)},
-\\
+ \label{eq:Ea} \\
\tilde{w}(a) &= (a-1)w_a - (1+w_0 + w_a)\log\left(a\right),
\label{eq:friedmann}
\end{align}
diff --git a/theory/Cosmology/neutrinos.tex b/theory/Cosmology/neutrinos.tex
new file mode 100644
index 0000000000000000000000000000000000000000..b9a5330f45e4763b115f4b7574ecaf598ee98830
--- /dev/null
+++ b/theory/Cosmology/neutrinos.tex
@@ -0,0 +1,36 @@
+% Willem Elbers, 2019
+
+\subsection{Massive neutrinos}
+
+The non-relativistic transition of massive neutrinos occurs at $a^{-1}\approx 1890 (m_\nu/1\text{ eV})$, which changes the background evolution, affecting the Hubble rate $E(a)$ and therefore most integrated quantities described in this document. Users can include this effect by specifying the number of massive neutrino species $N_\nu$ and their nonzero neutrino masses $m_{\nu,i}$ in eV $(m_{\nu,i}\neq 0, i=1,\dots,N_\nu)$ and (optional) degeneracies $g_i$. In addition, users have the option to specify the present-day neutrino temperature $T_{\nu,0}$\footnote{To match the neutrino density from an accurate calculation of decoupling \citep{Mangano2005}, one can use the value $T_{\nu,0}/T_{\mathrm{CMB},0}=0.71599$ \citep{Lesgourgues2011}.} and and an effective number of ultra-relativistic (massless) species $N_\mathrm{ur}$. Together with the present-day CMB temperature $T_{\mathrm{CMB},0}$, these parameters are used to compute the photon density $\Omega_\gamma$, the ultra-relativistic species density $\Omega_\mathrm{ur}$, and the massive neutrino density $\Omega_\nu(a)$, replacing the total radiation density parameter $\Omega_r$.
+
+In our conventions, the massive neutrino contribution at $a=1$ is \emph{not} included in the present-day matter density $\Omega_m=\Omega_c+\Omega_b$. The radiation term appearing in (\ref{eq:Ea}) is simply replaced by
+\begin{align}
+ \Omega_r a^{-4} &= \left[\Omega_\gamma + \Omega_\mathrm{ur} + \Omega_\nu(a)\right] a^{-4}.
+\end{align}
+The CMB density is constant and given by
+%
+\begin{align}
+ \Omega_\gamma &= \frac{\pi^2}{15}\frac{(k_b T_\text{CMB,0})^4}{(\hbar c)^3}\frac{1}{\rho_{\rm crit}c^2},
+\end{align}
+The ultra-relativistic density is given by
+\begin{align}
+ \Omega_\mathrm{ur} = \frac{7}{8}\left(\frac{4}{11}\right)^{4/3} N_\mathrm{ur}\,\Omega_\gamma.
+\end{align}
+Note that we assume instantaneous decoupling for the ultra-relativistic species. The time-dependent massive neutrino density parameter is \citep{Zennaro2016}
+\begin{align}
+ \Omega_\nu(a) = \Omega_\gamma \sum_{i=1}^{N_\nu}\frac{15}{\pi^4}g_i\left(\frac{T_{\nu,0}}{T_\text{CMB}}\right)^4 \mathcal{F}\left(\frac{a m_{\nu,i}}{k_b T_{\nu,0}}\right), \label{eq:nudensity}
+\end{align}
+where the function $\mathcal{F}$ is given by the momentum integral
+%
+\begin{align}
+ \mathcal{F}(y) = \int_0^{\infty} \frac{x^2\sqrt{x^2+y^2}}{1+e^{x}}\mathrm{d}x.
+\end{align}
+As $\Omega_\nu(a)$ is needed to compute other cosmological integrals, this function should be calculated with greater accuracy. At the start of the simulation, values of (\ref{eq:nudensity}) are tabulated on a piecewise linear grid of $2\times3\times10^4$ values of $a$ spaced between $\log(a_{\nu,\text{begin}})$, $\log(a_{\nu,\text{mid}})$, and $\log(a_{\nu,\text{end}}) = \log(1)=0$. The value of $a_{\nu,\text{begin}}$ is automatically chosen such that the neutrinos are still relativistic at the start of the table. The value of $\log(a_{\nu,\text{mid}})$ is chosen just before the start of the simulation. The integrals $\mathcal{F}(y)$ are evaluated using the 61-points Gauss-Konrod rule implemented in the {\sc gsl} library with a relative error limit of $\epsilon=10^{-13}$. Tabulated values are then linearly interpolated whenever $E(a)$ is computed.
+
+For completeness, we also compute the effective number of relativistic degrees of freedom at high redshift:
+\begin{align}
+ N_\text{eff} = N_\mathrm{ur} + N_\nu\left(\frac{11}{4}\right)^{4/3}\left(\frac{T_{\nu,0}}{T_\text{CMB,0}}\right)^4.
+\end{align}
+
+Neutrino effects also play a role at the perturbation level. These can be included in SWIFT using the fully non-linear $\delta f$ method \citep{Elbers2020}.