including pmf turbulence heating

src/pyhyrec/params.py

@@ -111,6 +111,9 @@ def __init__(self, new_params = None) -> None:
         "decay" : 0.0,
         "on_the_spot" : 0.0,
         "Mpbh" : 1.0,
-        "fpbh" : 0.0,}
+        "fpbh" : 0.0,
+        "sigmaB_PMF" : 0.0,
+        "nB_PMF" : -2.0,
+        "sigmaA_PMF" : 31.78}
 
         super().__init__(new_params)

src/pyhyrec/src/energy_injection.c

@@ -87,15 +87,30 @@ double decay_rate_pmf_turbulences(double z, double tdti, double nB)
   return 3.0*m/2.0 * pow(log(1+tdti), m)/pow(log(1+tdti) + 1.5 * log((1+zi)/(1+z)), m+1);
 }
 
+// Variance of the pmf power spectrum at the Jean's scale
+double sigma_Jeans_pmf(double obh2, double ocbh2)
+{
+  return 2.116 * sqrt(obh2 / 0.02242) * sqrt(ocbh2 / 0.1424);
+}
 
-// Energy injection rate due to ambipolar diffusion
-double gamma_ambipolar(double z, double xp, double B0, double nB)
+/*  Energy injection rate due to ambipolar diffusion (result is in eV / cm^3 / s)*/
+double dEdtdV_heat_turbulences_pmf(double z, double H, double obh2, double ocbh2, double sigmaA, double sigmaB, double nB)
 {
-  return 0.0;
+
+  double zi = 1088;
+  double smooth = (1.0-tanh((z - zi)/50.0))/2.0; // smoothing the introduction of energy injection from PMF
+
+  double en =  1e-25 / (8.0 * M_PI) * 6.241509e+18 * pow(1+z, 4) * sigmaA * sigmaA * pow(sigmaB/sigmaA, 4.0/(5.0+nB)); // in units of eV / cm^3
+  double tdti = sigma_Jeans_pmf(obh2, ocbh2) / sigmaA;
+  return en * decay_rate_pmf_turbulences(z, tdti, nB) * H * smooth;
 }
 
 
 
+
+
+
+
 /***************************************************************************************
 Effect of accreting primordial black holes 
 Since the accuracy is certainly not at the percent level, 
@@ -257,21 +272,41 @@ void update_dEdtdV_dep(double z_out, double dlna, double xe, double Tgas,
 		       double nH, double xH, double H, REC_COSMOPARAMS *params, double *dEdtdV_dep, 
            double *ion, double *exclya, double *dEdtdV_heat) {
 
+  // Injected energy via particle photons and electrons cascades
   double inj  = dEdtdV_inj(z_out, xe, Tgas, params->inj_params);
 
-  if (params->inj_params->on_the_spot == 1) *dEdtdV_dep = inj;
+  if (params->inj_params->on_the_spot == 1) 
+  {
+    *dEdtdV_dep = inj;
+  }
+  else { 
 
-  // Else put in your favorite recipe to translate from injection to deposition
-  // Here I assume injected photon spectrum Compton cools at rate dE/dt = - 0.1 n_h c sigma_T E
-  // This is valid for E_photon ~ MeV or so.
-
-  else { // c sigma_T = 2e-14 (cgs)
+    // Else put in your favorite recipe to translate from injection to deposition
+    // Here I assume injected photon spectrum Compton cools at rate dE/dt = - 0.1 n_h c sigma_T E
+    // This is valid for E_photon ~ MeV or so.
+    
+    // c sigma_T = 2e-14 (cgs)
     *dEdtdV_dep = (*dEdtdV_dep *exp(-7.*dlna) + 2e-15* dlna*nH/H *inj)
                  /(1.+ 2e-15 *dlna*nH/H);                              
   }     
 
   *ion     = *dEdtdV_dep / 3. / nH * xH /EI;
   *exclya  = *ion / 0.75;
-  *dEdtdV_heat    = *dEdtdV_dep * chi_heat(xe); 
+  *dEdtdV_heat = *dEdtdV_dep * chi_heat(xe);
+
+  // add the primordial magnetic field contribution
+  double sigmaB = params->inj_params->sigmaB_PMF;
+
+  if (sigmaB > 0)
+  {
+    double sigmaA = params->inj_params->sigmaA_PMF;
+    double nB = params->inj_params->nB_PMF;
+
+
+    *dEdtdV_heat = *dEdtdV_heat + dEdtdV_heat_turbulences_pmf(z_out, H, params->obh2, params->ocbh2, sigmaA, sigmaB, nB);
+
+    //printf("we are here : %e %e %e %e %e\n ", sigmaB, sigmaA, nB, *dEdtdV_heat, dEdtdV_heat_turbulences_pmf(z_out, H, params->obh2, params->ocbh2, sigmaA, sigmaB, nB));
+
+  }
 }
 

src/pyhyrec/src/energy_injection.h

@@ -30,6 +30,8 @@ typedef struct {
 
   double ion, exclya, dEdtdV_heat;   /* Adding the possibility to have a heating decorrelated from ion and exclya */
 
+  double sigmaB_PMF, nB_PMF;      /* adding the possibility for Primordial Magnetic Field Heating (sigmaB_PMF in nG)*/
+  double sigmaA_PMF, sigmaJ_PMF;  /* characteristic amplitude of the PMF on Alfven's scale and Jean's scale */
 } INJ_PARAMS;
 
 
@@ -60,5 +62,6 @@ void update_dEdtdV_dep(double z_out, double dlna, double xe, double Tgas,
 		       double nH, double xH, double H, REC_COSMOPARAMS *params, double *dEdtdV_dep, 
            double *dEdtdV_ion, double *dEdtdV_exclya, double *dEdtdV_heat);
 double decay_rate_pmf_turbulences(double z, double tdti, double nB); /* decay rate of PMF turbulences*/
+double dEdtdV_heat_turbulences_pmf(double z, double H, double obh2, double ocbh2, double sigmaA, double sigmaB, double nB);
 
 #endif

src/pyhyrec/src/history.c

@@ -111,6 +111,9 @@ void init_hyrec(REC_COSMOPARAMS * param, INPUT_COSMOPARAMS cosmo_params, INPUT_I
   param->inj_params->decay       = injection_params.decay;
   param->inj_params->Mpbh        = injection_params.Mpbh;
   param->inj_params->fpbh        = injection_params.fpbh;
+  param->inj_params->sigmaB_PMF  = injection_params.sigmaB_PMF;
+  param->inj_params->nB_PMF      = injection_params.nB_PMF;
+  param->inj_params->sigmaA_PMF  = injection_params.sigmaA_PMF;
 
   param->inj_params->odmh2 = param->ocbh2 - param->obh2;
 
@@ -134,7 +137,7 @@ HYREC_DATA * run_hyrec(INPUT_COSMOPARAMS cosmo_params, INPUT_INJ_PARAMS inj_para
   //sa.sa_flags   = SA_SIGINFO;
 
   //sigaction(SIGSEGV, &sa, NULL);
-  //printf("We are here");
+  //printf("We are here, %e \n", inj_params.fpbh);
 
   HYREC_DATA * data = malloc(sizeof(*data));
 

src/pyhyrec/src/history.h

@@ -85,6 +85,8 @@ typedef struct {
   int on_the_spot;
   double Mpbh, fpbh;
   double decay;
+  double sigmaB_PMF, nB_PMF;
+  double sigmaA_PMF;
 } INPUT_INJ_PARAMS;
 
 void init_hyrec(REC_COSMOPARAMS * param, INPUT_COSMOPARAMS cosmo_params, INPUT_INJ_PARAMS injection_params);

src/pyhyrec/src/hydrogen.h

@@ -41,9 +41,8 @@
 #define PION_MAX  1e-2
 
 
-/****** CONSTANTS IN CGS + EV UNIT SYSTEM *******/
 
-#define EI   13.598286071938324        /* Hydrogen ionization energy in eV, reduced mass, no relativistic corrections */
+
 
 
 /****** CONSTANTS IN CGS + EV UNIT SYSTEM *******/