accounting for temperature dependence of the scattering cross section and...

accounting for temperature dependence of the scattering cross section and higher multipole coefficients

Makefile

@@ -38,7 +38,7 @@ OPTFLAG = -O4 -ffast-math #-march=native
 #OPTFLAG = -fast
 
 # your openmp flag (comment for compiling without openmp)
-#OMPFLAG   = -fopenmp
+OMPFLAG   = -fopenmp
 #OMPFLAG   = -mp -mp=nonuma -mp=allcores -g
 #OMPFLAG   = -openmp
 

changelog.log

@@ -2,42 +2,87 @@
 -----------------------------------------
 * this file lists all changes on order to incorporate DM-neutrino scattering
   into CLASS
+* This file lists all changes w.r.t the class_public version, downloaded 2/10/10
+* the commit this code is based on is
+  ----------
+  commit 4be2cb31bb87a5c90d8ae5277ae22252a4732676
+  Merge: de7d9d7 55a99c6
+  Author: lesgourg <lesgourg@physik.rwth-aachen.de>
+  Date:   Mon Sep 10 18:06:10 2018 +0200
+  ----------
 
 1. 'thermodynamics.h'
 ----------------------------------------
-* in struct thermo
-  line 154: double u_nuDM
-	 -> add the coupling strength as a new parameter, the mass is normalized to 100 GeV
+* in struct thermo: add new parameters
+  line 154: double u_nuDM_0
+	 -> interaction cross section normalizedt to a DM mass of 100 GeV,
+	    and the Thomson scattering rate at current time
+  line 155: double n_nuDM;
+	 -> temperature dependence of the DM-neutrino scattering cross section
+	    sigma_nuDM ~ T^n_nuDM
   line 156: short has_coupling_nuDM
 	 -> flag to indicate weather the new coupling should be taken into account
   line 178: int index_th_dmu_nuDM
 	 -> index for the DM-neutrino scattering rate
 
+2. 'perturbations.h'
+----------------------------------------
+* in struct perturb: add new parameters
+  line 103: double * alpha_nuDM;
+	 -> array containing the higer-multipole coefficients for the
+            interaction terms in the 'ur' Boltzmann hierarchy
+
 2. 'input.c'
 ----------------------------------------
+* in function input_init: print information about nu-DM interactions
+  line 495: if has_coupling_nuDM is TRUE print the following parameters
+	 -> u_nuDM_0
+	 -> n_nuDM
+	 -> l_max_ur
+	 -> >ur_fluid_trigger_tau_over_tau_k
+	 -> alpha_nuDM
+
 * in function input_read_parameters
-  line 1297: class_read_double("u_nuDM", pth->u_nuDM)
+  line 1309: class_read_double("u_nuDM_0", pth->u_nuDM_0)
 	  -> read strength of neutrino-DM cooupling
-  line 1298: if the coupling strength is greater than 0 set flag has_coupling_nuDM to TRUE
+  line 1311: if the coupling strength is greater than 0 set flag has_coupling_nuDM to TRUE
+  line 1312: if the coupling strength is greater than 0 read n_nuDM
+  line 2682: if has_coupling_nuDM is TRUE create array for multipole coefficients
+  line 2684: read in alpha_nuDM from the input file
+  line 2691: allocate the array for alpha_nuDM
+	  -> array is of length l_max_ur
+  line 2696: copy values read to alpha_nuDM
+  line 2698: if the # of values read from input is smaller than the l_max_ur copy
+             the largest multipole coefficient read from file to all higher ones
+  line 2702: if no values are given set all entries of alpha_nuDM to default=1.
+
 
-* in function input_default_parameters
-  line 3027: pth->u_nuDM=0.;
-  line 3028: pth->has_coupling_nuDM=_FALSE_;
+* in function input_default_parameters: set default values for all new parameters
+                                        but alpha_nuDM
+  line 3068: pth->u_nuDM_0=0.;
+  line 3069: pth->n_nuDM=0.;
+  line 3070: pth->has_coupling_nuDM=_FALSE_;
 
 3. 'thermodynamics.c'
 ----------------------------------------
 * in function 'thermodynamics_at_z'
-  line 158: pvecthermo[pth->index_th_dmu_nuDM] = (1.+z)*(1.+z)*pth->u_nuDM*3.*pba->H0*pba->H0/8./_PI_/_G_*pba->Omega0_cdm*pow(_c_,4)*_sigma_/1.e11/_eV_/_Mpc_over_m_; 
-	 -> if has_coupling_nuDM is true compute dmu_nuDM
+  line 158: if has_coupling_nuDM is true compute dmu_nuDM
+	 -> pvecthermo[pth->index_th_dmu_nuDM] = pow(1.+z, 2.+pth->n_nuDM)
+	    *pth->u_nuDM_0*3.*pba->H0*pba->H0/8./_PI_/_G_*pba->Omega0_cdm
+            *pow(_c_,4)*_sigma_/1.e11/_eV_/_Mpc_over_m_; 
 
 * in function 'thermodynamics_init'
   line 379: u_nuDM can't be smaller than zero
   line 383: DM-nu coupling is not implemented with non-flat curvature
   line 632: fill thermodynamics table with values for dmu_nDM
 	 -> only if has_coupling_nuDM is TRUE
+	 -> pth->thermodynamics_table[index_tau*pth->th_size+pth->index_th_dmu_nuDM]
+	    = 3./8./_PI_/_G_*pow(pth->z_table[index_tau]+1.,2.+pth->n_nuDM)
+            *pba->Omega0_cdm*pba->H0*pba->H0*pth->u_nuDM_0*pow(_c_,4)
+            *_sigma_/1.e11/_eV_/_Mpc_over_m_;
 
 * in function 'thermodynamics_indices'
-  line 952: asign index to pth->index_th_dmu_nuDM
+  line 953: asign index to pth->index_th_dmu_nuDM
 	 -> only if has_coupling_nuDM is TRUE
 
 
@@ -49,13 +94,28 @@
 * in function 'perturb_derivs'
   line 6918: double S_nuDM
   -> the following commands are executed only if pth->has_coupling_nuDM is TRUE
-  line 6983: S_nuDM = 4./3. *pvecback[pba->index_bg_rho_ur]/pvecback[pba->index_bg_rho_cdm];
-  line 7238: dy[pv->index_pt_theta_cdm] += S_nuDM*pvecthermo[pth->index_th_dmu_nuDM]*(y[pv->index_pt_theta_ur]-Y[pv->index_pt_theta_cdm]);
-	  -> additional condition ppw->approx[ppw->index_ap_rsa] == (int)rsa_off
-  line 7389: dy[pv->index_pt_theta_ur] += pvecthermo[pth->index_th_dmu_nuDM]*(y[pv->index_pt_theta_cdm]-y[pv->index_pt_theta_ur]);
-  line 7401: dy[pv->index_pt_shear_ur] += -9./10.*pvecthermo[pth->index_th_dmu_nuDM]*y[pv->index_pt_shear_ur];	
-  line 7408: dy[pv->index_pt_l3_ur] -= pvecthermo[pth->index_th_dmu_nuDM]*y[pv->index_pt_l3_ur];
-  line 7415: dy[pv->index_pt_delta_ur+l] -= pvecthermo[pth->index_th_dmu_nuDM]*y[pv->index_pt_delta_ur+l];
-  line 7423: dy[pv->index_pt_delta_ur+l] -= pvecthermo[pth->index_th_dmu_nuDM]*y[pv->index_pt_delta_ur+l];
+  line 6983: S_nuDM = 4./3.*pvecback[pba->index_bg_rho_ur]/pvecback[pba->index_bg_rho_cdm];
+  line 7238: if ppw->approx[ppw->index_ap_rsa] == (int)rsa_off
+	     dy[pv->index_pt_theta_cdm] += S_nuDM*pvecthermo[pth->index_th_dmu_nuDM]
+	     *(y[pv->index_pt_theta_ur]- y[pv->index_pt_theta_cdm]);
+  line 7389: if ppw->approx[ppw->index_ap_rsa] == (int)rsa_off
+	     dy[pv->index_pt_theta_ur] += pvecthermo[pth->index_th_dmu_nuDM]
+	     *(y[pv->index_pt_theta_cdm]-y[pv->index_pt_theta_ur]);
+  line 7401: if ppw->approx[ppw->index_ap_rsa] == (int)rsa_off
+	     and if ppw->approx[ppw->index_ap_ufa] == (int)ufa_off
+             dy[pv->index_pt_shear_ur] += -0.5*ppt->alpha_nuDM[2]
+             *pvecthermo[pth->index_th_dmu_nuDM]*y[pv->index_pt_shear_ur];
+  line 7408: if ppw->approx[ppw->index_ap_rsa] == (int)rsa_off
+	     and if ppw->approx[ppw->index_ap_ufa] == (int)ufa_off
+	     dy[pv->index_pt_l3_ur] -= ppt->alpha_nuDM[l]
+             *pvecthermo[pth->index_th_dmu_nuDM]*y[pv->index_pt_l3_ur];
+  line 7415: if ppw->approx[ppw->index_ap_rsa] == (int)rsa_off
+	     and if ppw->approx[ppw->index_ap_ufa] == (int)ufa_off
+             dy[pv->index_pt_delta_ur+l] -= ppt->alpha_nuDM[l]
+             *pvecthermo[pth->index_th_dmu_nuDM]*y[pv->index_pt_delta_ur+l];
+  line 7423: if ppw->approx[ppw->index_ap_rsa] == (int)rsa_off
+	     and if ppw->approx[ppw->index_ap_ufa] == (int)ufa_off
+             dy[pv->index_pt_delta_ur+l] -= ppt->alpha_nuDM[l]
+             *pvecthermo[pth->index_th_dmu_nuDM]*y[pv->index_pt_delta_ur+l];
 	  -> truncation equation for l_max
   line 7450: UFA equation (CLASS scheme): NO MODIFICATIONS!

include/perturbations.h

@@ -99,6 +99,9 @@ struct perturbs
    *  parameters and the content of the 'precision', 'background' and
    *  'thermodynamics' structures) */
 
+  /* for the use with neutrino-DM interactions: list of higher order coefficients */
+  double * alpha_nuDM;
+
   //@{
 
   short has_perturbations; /**< do we need to compute perturbations at all ? */
@@ -446,7 +449,6 @@ struct perturb_vector
   int * used_in_sources; /**< boolean array specifying which
                             perturbations enter in the calculation of
                             source functions */
-
 };
 
 

include/thermodynamics.h

@@ -151,8 +151,8 @@ struct thermo
   double annihilation_f_halo; /** takes the contribution of DM annihilation in halos into account*/
   double annihilation_z_halo; /** characteristic redshift for DM annihilation in halos*/
 
-  double u_nuDM;              /** DM-NU scattering cross section normalized to the Thomson cross
-				  section and a DM mass of 100 GeV */
+  double u_nuDM_0;            /** DM-NU scattering cross section normalized to the Thomson crosssection and a DM mass of 100 GeV at T=T_0 */
+  double n_nuDM;             /** temperature dependance of the cross section: sigma ~ T^n_nuDM*/
   short has_coupling_nuDM;   /** flag to indicate weather the scattering is present */
   //@}
 

nuDM.ini

@@ -22,16 +22,17 @@ N_ur = 3.046
 -------------------
 ---> g-cdm coupling
 -------------------
-u_nuDM        = 1.e-5
-gauge         = newtonian
-k_pivot       = 0.05
-sBBN file     = ./bbn/sBBN.dat
+u_nuDM_0   = 1.e-5
+gauge      = newtonian
+k_pivot    = 0.05
+sBBN file  = ./bbn/sBBN.dat
+alpha_nuDM = 0., 1., 1.5    #the firsts to values are for detla and theta and not used!
 
 -----------
 ---> OUTPUT
 -----------
 output               = tCl, pCl, mPk
-root                 = output/2018-10-24_ufa-trigger/nDM_u1em5_eufa0.1_
+root                 = output/test_
 #root                 = output/2018-10-24_mixed-damping/lcdm_
 format               = camb
 k_per_decade_for_pk  = 80
@@ -44,8 +45,8 @@ write thermodynamics = no
 -----------------------
 ---> precision settings
 -----------------------
-ur_fluid_trigger_tau_over_tau_k = 100.
-l_max_ur = 100
+ur_fluid_trigger_tau_over_tau_k = 10.
+l_max_ur = 10
 
 -----------------------
 ---> breake the silence

source/input.c

@@ -491,6 +491,17 @@ int input_init(
     }
   }
 
+  if(pth->has_coupling_nuDM==_TRUE_){
+    printf(" -> Neutrino-Dark Matter coupling enabled with the following parameters:\n");
+    printf("     u_nuDM_0 = %f\n", pth->u_nuDM_0);
+    printf("       n_nuDM = %f\n", pth->n_nuDM);
+    printf("     l_max_ur = %d\n", ppr->l_max_ur);
+    printf("   ufa_triggr = %f\n", ppr->ur_fluid_trigger_tau_over_tau_k);
+    printf("   alpha_nuDM = ");
+    for (i=0; i<ppr->l_max_ur; i++)
+      printf("%f, ", ppt->alpha_nuDM[i]);
+    printf("\n");    
+  }
   return _SUCCESS_;
 
 }
@@ -1294,9 +1305,11 @@ int input_read_parameters(
     }
   }
   /** parameters for DM-nu scattering */
-  class_read_double("u_nuDM", pth->u_nuDM);
-  if (pth->u_nuDM>0.)
+  class_read_double("u_nuDM_0", pth->u_nuDM_0);
+  if (pth->u_nuDM_0>0.){
     pth->has_coupling_nuDM = _TRUE_;
+    class_read_double("n_nuDM", pth->n_nuDM)
+  }
   
 
   /** (c) define which perturbations and sources should be computed, and down to which scale */
@@ -2663,6 +2676,35 @@ int input_read_parameters(
   class_read_int("l_max_pol_g",ppr->l_max_pol_g);
   class_read_int("l_max_dr",ppr->l_max_dr);
   class_read_int("l_max_ur",ppr->l_max_ur);
+
+  /** this is a bit special: if nu-DM scattering is enabled initialize the list of hihger-order
+      multipole coefficients, set default to 1 and read in vaues */
+  if (pth->has_coupling_nuDM==_TRUE_)
+    {
+      double * read_alpha;
+      parser_read_list_of_doubles(pfc,
+				  "alpha_nuDM",
+				  &entries_read,
+				  &(read_alpha),
+				  &flag2,
+				  errmsg);
+      class_alloc(ppt->alpha_nuDM,
+		  ppr->l_max_ur*sizeof(double),
+		  errmsg);
+      if(flag2==_TRUE_){
+	for(i=0; i<entries_read; i++)
+	  ppt->alpha_nuDM[i] = read_alpha[i];
+	for(i=entries_read; i<ppr->l_max_ur; i++)
+	  ppt->alpha_nuDM[i] = read_alpha[entries_read-1];
+      }
+      else{
+	for(i=0; i<ppr->l_max_ur; i++)
+	  ppt->alpha_nuDM[i] = 1.;
+      }
+      free(read_alpha);
+    }
+  /** end of the nu-DM section */
+      
   if (pba->N_ncdm>0)
     class_read_int("l_max_ncdm",ppr->l_max_ncdm);
   class_read_int("l_max_g_ten",ppr->l_max_g_ten);
@@ -3023,7 +3065,8 @@ int input_default_params(
 
   pth->compute_damping_scale = _FALSE_;
 
-  pth->u_nuDM=0.;
+  pth->u_nuDM_0=0.;
+  pth->n_nuDM=0.;
   pth->has_coupling_nuDM=_FALSE_;
 
   /** - perturbation structure */

source/perturbations.c

@@ -7398,21 +7398,21 @@ int perturb_derivs(double tau,
                  // non-standard term, non-zero if cvis2_ur not 1/3
                  -(1.-ppt->three_cvis2_ur)*(8./15.*(y[pv->index_pt_theta_ur]+metric_shear)));
 	  if(pth->has_coupling_nuDM==_TRUE_)
-	    dy[pv->index_pt_shear_ur] += -9./10.*pvecthermo[pth->index_th_dmu_nuDM]*y[pv->index_pt_shear_ur];
+	    dy[pv->index_pt_shear_ur] += -0.5*ppt->alpha_nuDM[2]*pvecthermo[pth->index_th_dmu_nuDM]*y[pv->index_pt_shear_ur];
 
           /** - -----> Exact ur l=3 */
           l = 3;
           dy[pv->index_pt_l3_ur] = k/(2.*l+1.)*
             (l*2.*s_l[l]*s_l[2]*y[pv->index_pt_shear_ur]-(l+1.)*s_l[l+1]*y[pv->index_pt_l3_ur+1]);
 	  if(pth->has_coupling_nuDM==_TRUE_)
-	    dy[pv->index_pt_l3_ur] -= pvecthermo[pth->index_th_dmu_nuDM]*y[pv->index_pt_l3_ur];
+	    dy[pv->index_pt_l3_ur] -= ppt->alpha_nuDM[l]*pvecthermo[pth->index_th_dmu_nuDM]*y[pv->index_pt_l3_ur];
 
           /** - -----> Exactc ur l>3 */
           for (l = 4; l < pv->l_max_ur; l++) {
             dy[pv->index_pt_delta_ur+l] = k/(2.*l+1)*
               (l*s_l[l]*y[pv->index_pt_delta_ur+l-1]-(l+1.)*s_l[l+1]*y[pv->index_pt_delta_ur+l+1]);
 	    if(pth->has_coupling_nuDM==_TRUE_)
-	      dy[pv->index_pt_delta_ur+l] -= pvecthermo[pth->index_th_dmu_nuDM]*y[pv->index_pt_delta_ur+l];
+	      dy[pv->index_pt_delta_ur+l] -= ppt->alpha_nuDM[l]*pvecthermo[pth->index_th_dmu_nuDM]*y[pv->index_pt_delta_ur+l];
 	  }
 
           /** - -----> Exact ur lmax_ur */
@@ -7420,7 +7420,7 @@ int perturb_derivs(double tau,
           dy[pv->index_pt_delta_ur+l] =
             k*(s_l[l]*y[pv->index_pt_delta_ur+l-1]-(1.+l)*cotKgen*y[pv->index_pt_delta_ur+l]);
 	  if(pth->has_coupling_nuDM==_TRUE_)
-	    dy[pv->index_pt_delta_ur+l] -= pvecthermo[pth->index_th_dmu_nuDM]*y[pv->index_pt_delta_ur+l];
+	    dy[pv->index_pt_delta_ur+l] -= ppt->alpha_nuDM[l]*pvecthermo[pth->index_th_dmu_nuDM]*y[pv->index_pt_delta_ur+l];
         }
 
         else {

source/thermodynamics.c

@@ -153,10 +153,10 @@ int thermodynamics_at_z(
     pvecthermo[pth->index_th_ddg]=0.;
 
     /* calculate the rate for DM-neutrino scattering */
-    /* dmu_nuDM = 3/(8*Pi*G)*sigma_th/10^11eV*c^4/Mpc_over_m*(1+z)^2*u_nuDM*Omega_DM*H_0^2 */
+    /* dmu_nuDM = 3/(8*Pi*G)*sigma_th/10^11eV*c^4/Mpc_over_m*(1+z)^(2+n_nuDM)*u_nuDM_0*Omega_DM*H_0^2 */
     /* actually this expression is exact */
     if (pth->has_coupling_nuDM==_TRUE_)
-      pvecthermo[pth->index_th_dmu_nuDM] = (1.+z)*(1.+z)*pth->u_nuDM*3.*pba->H0*pba->H0/8./_PI_/_G_*pba->Omega0_cdm*pow(_c_,4)*_sigma_/1.e11/_eV_/_Mpc_over_m_; 
+      pvecthermo[pth->index_th_dmu_nuDM] = pow(1.+z, 2.+pth->n_nuDM)*pth->u_nuDM_0*3.*pba->H0*pba->H0/8./_PI_/_G_*pba->Omega0_cdm*pow(_c_,4)*_sigma_/1.e11/_eV_/_Mpc_over_m_; 
 
     /* Calculate Tb */
     pvecthermo[pth->index_th_Tb] = pba->T_cmb*(1.+z);
@@ -376,7 +376,7 @@ int thermodynamics_init(
              "CDM decay effects require the presence of CDM!");
 
   /* test nuDM scattering parameters */
-  class_test((pth->u_nuDM<0.),
+  class_test((pth->u_nuDM_0<0.),
 	     pth->error_message,
 	     "DM-nu coupling strength can't be smaller than zero!");
 
@@ -631,7 +631,9 @@ int thermodynamics_init(
   /* compute dmu_nuDM */
   if(pth->has_coupling_nuDM==_TRUE_){
     for (index_tau=0; index_tau < pth->tt_size; index_tau++) {
-      pth->thermodynamics_table[index_tau*pth->th_size+pth->index_th_dmu_nuDM] = 3./8./_PI_/_G_*(pth->z_table[index_tau]+1.)*(pth->z_table[index_tau]+1.)*pba->Omega0_cdm*pba->H0*pba->H0*pth->u_nuDM*pow(_c_,4)*_sigma_/1.e11/_eV_/_Mpc_over_m_;
+      pth->thermodynamics_table[index_tau*pth->th_size+pth->index_th_dmu_nuDM]
+	= 3./8./_PI_/_G_*pow(pth->z_table[index_tau]+1.,2.+pth->n_nuDM)*pba->Omega0_cdm*pba->H0*pba->H0*pth->u_nuDM_0
+	*pow(_c_,4)*_sigma_/1.e11/_eV_/_Mpc_over_m_;
     }
   }