diff --git a/hi_class.ini b/hi_class.ini index c32de95e..41b00aa6 100644 --- a/hi_class.ini +++ b/hi_class.ini @@ -79,18 +79,29 @@ expansion_smg = 0.5 #this value will be overwritten using the closure equation. ii) "propto_scale": all the alphas are proportional to the scale factor. Then, you have to provide a vector containg the factor of proportionality of each alpha and the initial value of the Planck Mass - iii) "planck_linear": this parametrization depends on one single parameter. It was used - by the Planck collaboration in 1502.01590 (this reference provides all the details - about this parametrization) - iv) "planck_exponential": this parametrization depends on two parameters. It was used - by the Planck collaboration in 1502.01590 (this reference provides all the details - about this parametrization) + iii) "eft_alphas_power_law": delta_M, alpha_{K, B, T} are proportional to the scale factor to + some power. Then, you have to provide a vector containg eight parameters, the four + proportionality constants and the four exponents of the scale factor. This generalize + the "propto_scale" parametrization with the remarkable exception that in that case + alpha_M was proportional to the scale factor, while here it is the Planck mass that + scales with the scale factor + iv) "eft_gammas_power_law": this parametrization uses the basis of the eft framework + (Omega, gamma_i) and assumes that they are power laws of the scale factor + (g1_0*pow(a,g1_exp)). You have to provide a vector containing eight parameters, + the four proportionality constants and the four exponents of the scale factor + v) "eft_gammas_exponential": this parametrization uses the basis of the eft framework + (Omega, gamma_i) and assumes that they are exponentials in the scale factor. In a + schematic way it reads: exp(gi_0*pow(a,gi_exp))-1. You have to provide a vector + containing eight parameters, the four proportionality constants and the four + exponents of the scale factor The parametrization you want to used is stored in the variable "gravity_model", while the value of the parameters is stored in the vector "parameters_smg" i) "propto_omega" -> x_k, x_b, x_m, x_t, M*^2_ini (default) ii) "propto_scale" -> x_k, x_b, x_m, x_t, M*^2_ini - iii) "planck_linear" -> Omega_0 - iv) "planck_exponential" -> alpha_M0, beta + iii) "eft_alphas_power_law" -> delta_M_0, x_k, x_b, x_t, delta_M_0_exp, x_k_exp, x_b_exp, x_t_exp + iv) "eft_gammas_power_law" -> Om_0, g_1, g_2, g_3, Om_0_exp, g_1_exp, g_2_exp, g_3_exp + v) "eft_gammas_exponential" -> Om_0, g_1, g_2, g_3, Om_0_exp, g_1_exp, g_2_exp, g_3_exp + gravity_model = propto_omega parameters_smg = 1., 0., 0., 0., 1. diff --git a/include/background.h b/include/background.h index f00901bf..7c9f33a2 100755 --- a/include/background.h +++ b/include/background.h @@ -11,7 +11,7 @@ #include "parser.h" enum spatial_curvature {flat,open,closed}; -enum gravity_model {propto_omega, propto_scale, planck_linear, planck_exponential}; //write here the different models +enum gravity_model {propto_omega, propto_scale, eft_alphas_power_law, eft_gammas_power_law, eft_gammas_exponential}; //write here the different models // initial conditions for the perturbations enum pert_initial_conditions {single_clock, zero}; diff --git a/source/background.c b/source/background.c index 46597a7b..7c5808bd 100755 --- a/source/background.c +++ b/source/background.c @@ -2195,16 +2195,17 @@ int background_initial_conditions( pvecback_integration[pba->index_bi_M_pl_smg] = pba->parameters_2_smg[4]; break; - case planck_linear: - //NOTE: The Planck collaboration decided to consider models with M_pl^2=1+\Omega. Here, even if we have an additional integration parameter (i.e. M_pl_ini), we decided to fix it to M_pl^2=1+a*Omega_0 in order to be coherent with the choice of the Planck collaboration. - pvecback_integration[pba->index_bi_M_pl_smg] = 1+a*pba->parameters_2_smg[0]; + case eft_alphas_power_law: + pvecback_integration[pba->index_bi_M_pl_smg] = 1. + pba->parameters_2_smg[0]*pow(a, pba->parameters_2_smg[4]); break; - case planck_exponential: - //NOTE: The Planck collaboration decided to consider models with M_pl^2=1+\Omega. Here, even if we have an additional integration parameter (i.e. M_pl_ini), we decided to fix M_pl^2=exp(alpha_M0*pow(a, beta)/beta) in order to be coherent with the choice of the Planck collaboration. - pvecback_integration[pba->index_bi_M_pl_smg] = exp(pba->parameters_2_smg[0]*pow(a, pba->parameters_2_smg[1])/pba->parameters_2_smg[1]); + case eft_gammas_power_law: + pvecback_integration[pba->index_bi_M_pl_smg] = 1. + pba->parameters_2_smg[0]*pow(a,pba->parameters_2_smg[4]) + pba->parameters_2_smg[3]*pow(a,pba->parameters_2_smg[7]); break; - + case eft_gammas_exponential: + pvecback_integration[pba->index_bi_M_pl_smg] = exp(pba->parameters_2_smg[0]*pow(a,pba->parameters_2_smg[4])) + exp(pba->parameters_2_smg[3]*pow(a,pba->parameters_2_smg[7])) -1.; + break; + } if (pba->M_pl_evolution_smg == _TRUE_) @@ -2558,29 +2559,66 @@ int background_gravity_functions( pvecback[pba->index_bg_mpl_running_smg] = c_m*a; pvecback[pba->index_bg_M2_smg] = M_pl; } - else if (pba->gravity_model_smg == planck_linear) { - //NOTE: With this parametrization every function it is expressed analytically. Then, it is possible to choose both to take the derivative of M_pl to obtain alpha_M or to integrate alpha_M to obtain M_pl. Even if the two results are undistinguishable, we choose the latter option, since in Class integrals are more stable numerically. + else if (pba->gravity_model_smg == eft_alphas_power_law) { - double Omega = a*pba->parameters_2_smg[0]; + double M_0 = pba->parameters_2_smg[0]; + double c_k = pba->parameters_2_smg[1]; + double c_b = pba->parameters_2_smg[2]; + double c_t = pba->parameters_2_smg[3]; + double M_0_exp = pba->parameters_2_smg[4]; + double c_k_exp = pba->parameters_2_smg[5]; + double c_b_exp = pba->parameters_2_smg[6]; + double c_t_exp = pba->parameters_2_smg[7]; - pvecback[pba->index_bg_tensor_excess_smg] = 0.; - pvecback[pba->index_bg_mpl_running_smg] = Omega/(1.+Omega); - pvecback[pba->index_bg_braiding_smg] = -pvecback[pba->index_bg_mpl_running_smg]; - pvecback[pba->index_bg_kineticity_smg] = 3.*((2.+3.*Omega)*(pvecback[pba->index_bg_rho_smg]+pvecback[pba->index_bg_p_smg])+Omega*(pvecback[pba->index_bg_rho_tot_wo_smg]+pvecback[pba->index_bg_p_tot_wo_smg]))/2./(1.+Omega)/rho_tot; + pvecback[pba->index_bg_kineticity_smg] = c_k*pow(a, c_k_exp); + pvecback[pba->index_bg_braiding_smg] = c_b*pow(a, c_b_exp); + pvecback[pba->index_bg_tensor_excess_smg] = c_t*pow(a, c_t_exp); + pvecback[pba->index_bg_mpl_running_smg] = M_0*M_0_exp*pow(a, M_0_exp)/(1. + M_0*pow(a, M_0_exp)); pvecback[pba->index_bg_M2_smg] = M_pl; } - else if (pba->gravity_model_smg == planck_exponential) { - //NOTE: With this parametrization every function it is expressed analytically. Then, it is possible to choose both to take the derivative of M_pl to obtain alpha_M or to integrate alpha_M to obtain M_pl. Even if the two results are undistinguishable, we choose the latter option, since in Class integrals are more stable numerically. + else if ((pba->gravity_model_smg == eft_gammas_power_law) || (pba->gravity_model_smg == eft_gammas_exponential)) { + + double Omega=0., g1=0., g2=0., g3=0., Omega_p=0., Omega_pp=0., g3_p=0.; + double Omega_0=0., g1_0=0., g2_0=0., g3_0=0.; + double Omega_exp=0., g1_exp=0., g2_exp=0., g3_exp=0.; + + Omega_0 = pba->parameters_2_smg[0]; + g1_0 = pba->parameters_2_smg[1]; + g2_0 = pba->parameters_2_smg[2]; + g3_0 = pba->parameters_2_smg[3]; + Omega_exp = pba->parameters_2_smg[4]; + g1_exp = pba->parameters_2_smg[5]; + g2_exp = pba->parameters_2_smg[6]; + g3_exp = pba->parameters_2_smg[7]; + + if (pba->gravity_model_smg == eft_gammas_power_law) { + Omega = Omega_0*pow(a,Omega_exp); + Omega_p = Omega_0*Omega_exp*pow(a,Omega_exp-1.); // Derivative w.r.t. the scale factor + Omega_pp = Omega_0*Omega_exp*(Omega_exp-1.)*pow(a,Omega_exp-2.); // Derivative w.r.t. the scale factor + g1 = g1_0*pow(a,g1_exp); + g2 = g2_0*pow(a,g2_exp); + g3 = g3_0*pow(a,g3_exp); + g3_p = g3_0*g3_exp*pow(a,g3_exp-1.); // Derivative w.r.t. the scale factor + } + else { //(pba->gravity_model_smg == eft_gammas_exponential) + Omega = exp(Omega_0*pow(a,Omega_exp))-1.; + Omega_p = Omega_0*Omega_exp*pow(a,Omega_exp-1.)*exp(Omega_0*pow(a,Omega_exp)); // Derivative w.r.t. the scale factor + Omega_pp = Omega_0*Omega_exp*pow(a,Omega_exp-2.)*exp(Omega_0*pow(a,Omega_exp))*(Omega_exp-1.+Omega_0*Omega_exp*pow(a,Omega_exp)); // Derivative w.r.t. the scale factor + g1 = exp(g1_0*pow(a,g1_exp))-1.; + g2 = exp(g2_0*pow(a,g2_exp))-1.; + g3 = exp(g3_0*pow(a,g3_exp))-1.; + g3_p = g3_0*g3_exp*pow(a,g3_exp-1.)*exp(g3_0*pow(a,g3_exp)); // Derivative w.r.t. the scale factor + } + - double alpha_M0 = pba->parameters_2_smg[0]; - double beta = pba->parameters_2_smg[1]; - double Omega = exp(alpha_M0*pow(a, beta)/beta)-1; + double c_over_H2 = (-pow(a,2.)*Omega_pp + 3.*(Omega + a*Omega_p/2.)*(rho_tot + p_tot)/rho_tot + 3.*(pvecback[pba->index_bg_rho_smg]+pvecback[pba->index_bg_p_smg])/rho_tot)/2.; - pvecback[pba->index_bg_tensor_excess_smg] = 0; - pvecback[pba->index_bg_mpl_running_smg] = alpha_M0*pow(a, beta); - pvecback[pba->index_bg_braiding_smg] = -pvecback[pba->index_bg_mpl_running_smg]; - pvecback[pba->index_bg_kineticity_smg] = (1-beta-pvecback[pba->index_bg_mpl_running_smg])*pvecback[pba->index_bg_mpl_running_smg] + 3*(2+pvecback[pba->index_bg_mpl_running_smg])*(pvecback[pba->index_bg_rho_smg]+pvecback[pba->index_bg_p_smg])/2/rho_tot + 3*(pvecback[pba->index_bg_mpl_running_smg]+2*Omega/(1+Omega))*(pvecback[pba->index_bg_rho_tot_wo_smg]+pvecback[pba->index_bg_p_tot_wo_smg])/2/rho_tot; + pvecback[pba->index_bg_kineticity_smg] = 2.*(2.*g1*pow(pba->H0,2.)/rho_tot + c_over_H2)/(1. + Omega + g3); + pvecback[pba->index_bg_braiding_smg] = -(g2*pba->H0/sqrt(rho_tot) + a*Omega_p)/(1. + Omega + g3); + pvecback[pba->index_bg_tensor_excess_smg] = -g3/(1. + Omega + g3); + pvecback[pba->index_bg_mpl_running_smg] = a*(Omega_p + g3_p)/(1. + Omega + g3); pvecback[pba->index_bg_M2_smg] = M_pl; + } @@ -2645,10 +2683,23 @@ int background_gravity_parameters( pba->parameters_2_smg[4]); break; - case planck_linear: - printf("Modified gravity: planck_linear with parameters: \n"); - printf("-> Omega_0 = %g \n", - pba->parameters_2_smg[0]); + case eft_alphas_power_law: + printf("Modified gravity: eft_alphas_power_law with parameters: \n"); + printf("-> M_*^2_0 = %g, c_K = %g, c_B = %g, c_T = %g, M_*^2_exp = %g, c_K_exp = %g, c_B_exp = %g, c_T_exp = %g\n", + pba->parameters_2_smg[0],pba->parameters_2_smg[1],pba->parameters_2_smg[2],pba->parameters_2_smg[3], + pba->parameters_2_smg[4],pba->parameters_2_smg[5],pba->parameters_2_smg[6],pba->parameters_2_smg[7]); + break; + + case eft_gammas_power_law: + printf("Modified gravity: eft_gammas_power_law with parameters: \n"); + printf("-> Omega_0 = %g, gamma_1 = %g, gamma_2 = %g, gamma_3 = %g, Omega_0_exp = %g, gamma_1_exp = %g, gamma_2_exp = %g, gamma_3_exp = %g \n", + pba->parameters_2_smg[0],pba->parameters_2_smg[1],pba->parameters_2_smg[2],pba->parameters_2_smg[3],pba->parameters_2_smg[4],pba->parameters_2_smg[5],pba->parameters_2_smg[6],pba->parameters_2_smg[7]); + break; + + case eft_gammas_exponential: + printf("Modified gravity: eft_gammas_exponential with parameters: \n"); + printf("-> Omega_0 = %g, gamma_1 = %g, gamma_2 = %g, gamma_3 = %g, Omega_0_exp = %g, gamma_1_exp = %g, gamma_2_exp = %g, gamma_3_exp = %g \n", + pba->parameters_2_smg[0],pba->parameters_2_smg[1],pba->parameters_2_smg[2],pba->parameters_2_smg[3],pba->parameters_2_smg[4],pba->parameters_2_smg[5],pba->parameters_2_smg[6],pba->parameters_2_smg[7]); break; diff --git a/source/input.c b/source/input.c index 4a75ebce..935bf5c7 100755 --- a/source/input.c +++ b/source/input.c @@ -1041,28 +1041,37 @@ int input_read_parameters( class_read_list_of_doubles("parameters_smg",pba->parameters_2_smg,pba->parameters_2_size_smg); } - if (strcmp(string1,"planck_linear") == 0) { - pba->gravity_model_smg = planck_linear; + if (strcmp(string1,"eft_alphas_power_law") == 0) { + pba->gravity_model_smg = eft_alphas_power_law; pba->field_evolution_smg = _FALSE_; pba->M_pl_evolution_smg = _TRUE_; flag2=_TRUE_; - pba->parameters_2_size_smg = 1; + pba->parameters_2_size_smg = 8; class_read_list_of_doubles("parameters_smg",pba->parameters_2_smg,pba->parameters_2_size_smg); } - if (strcmp(string1,"planck_exponential") == 0) { - pba->gravity_model_smg = planck_exponential; + if (strcmp(string1,"eft_gammas_power_law") == 0) { + pba->gravity_model_smg = eft_gammas_power_law; pba->field_evolution_smg = _FALSE_; pba->M_pl_evolution_smg = _TRUE_; flag2=_TRUE_; - pba->parameters_2_size_smg = 2; + pba->parameters_2_size_smg = 8; class_read_list_of_doubles("parameters_smg",pba->parameters_2_smg,pba->parameters_2_size_smg); } - + + if (strcmp(string1,"eft_gammas_exponential") == 0) { + pba->gravity_model_smg = eft_gammas_exponential; + pba->field_evolution_smg = _FALSE_; + pba->M_pl_evolution_smg = _TRUE_; + flag2=_TRUE_; + pba->parameters_2_size_smg = 8; + class_read_list_of_doubles("parameters_smg",pba->parameters_2_smg,pba->parameters_2_size_smg); + } + class_test(flag2==_FALSE_, errmsg, - "could not identify gravity_theory value, check that it is one of 'propto_omega', 'propto_scale', 'planck_linear', 'planck_exponential' ..."); + "could not identify gravity_theory value, check that it is one of 'propto_omega', 'propto_scale', 'eft_alphas_power_law', 'eft_gammas_power_law', 'eft_gammas_exponential' ..."); }// end of loop over models