search

lesgourg / class_public

Public repository of the Cosmic Linear Anisotropy Solving System (master for the most recent version of the standard code; GW_CLASS to include Cosmic Gravitational Wave Background anisotropies; classnet branch for acceleration with neutral networks; ExoCLASS branch for exotic energy injection; class_matter branch for FFTlog)
223 stars 292 forks source link

Scale independent growth factor in DCDM scenario #436

Open ShazAlvi opened 3 years ago

ShazAlvi commented 3 years ago

In the process of working with growth factor in the decaying cold dark matter scenario, I discovered something that did not settle with me.

I am using the following LCDM parameters as the baseline case. 

params = {'output': 'tCl mPk dCl lCl','omega_b': 0.02237,'omega_cdm': 0.1201,'ln10^{10}A_s': 3.054, 'n_s': 0.9651,
                '100*theta_s': 1.04090, 'tau_reio': 0.0590, 'z_max_pk': 2.5, 'non linear': 'linear', 'lensing': 'yes','P_k_max_1/Mpc': 15.47}

The code used to compute the growth factor, in this case, is the following,

z = np.linspace(0, 2.5, nzmax)
MyClass = Class()
params = {'output': 'tCl mPk dCl lCl','omega_b': 0.02237,'omega_cdm': 0.1201,'ln10^{10}A_s': 3.054, 'n_s': 0.9651,
                '100*theta_s': 1.04090, 'tau_reio': 0.0590, 'z_max_pk': 2.5, 'non linear': 'linear', 'lensing': 'yes','P_k_max_1/Mpc': 15.47}
MyClass.set(params)
MyClass.compute()
for i, z_i in enumerate(z):
    D_0[i] = MyClass.scale_independent_growth_factor(z_i)

Then I implement the DCDM scenario in the following way,

MyClass.empty()
MyClass.struct_cleanup()
params={}
MyClass = Class()
params = {'output': 'tCl mPk dCl lCl','omega_b': 0.02237,'omega_cdm': 0.0,'ln10^{10}A_s': 3.054, 'n_s': 0.9651,
                '100*theta_s': 1.04090, 'tau_reio': 0.0590, 'z_max_pk': 2.5, 'non linear': 'linear', 'lensing': 'yes','P_k_max_1/Mpc': 15.47}
Y = 5.3687
Gamma_18 = 0.0
params['omega_ini_dcdm'] = Y*params['omega_b']
params['Gamma_dcdm'] = 30.85*Gamma_18
MyClass.set(params)
MyClass.compute()
for i, z_i in enumerate(z):
    D[i] = MyClass.scale_independent_growth_factor(z_i)

The way that this is implemented, the two cases should be equivalent to one another because Gamma_18 parameter is set to zero hence Gamma_dcdm is zero. Furthermore, the omega_ini_dcdm parameter is set so that it gives the correct omega_cdm density as used in the first LCDM case. The two growth factors from the two implementations are shown in the plot below, image I looked through the code and the difference seem to come from the factor $\rho_M$ in the following equation, image The difference occurs because $\rho_M$ factor does not include the dcdm density in the following lines in background.c in my version of CLASS,

rho_M = pvecback[pba->index_bg_rho_b];
  if (pba->has_cdm)
    rho_M += pvecback[pba->index_bg_rho_cdm];
  if (pba->has_idm_dr)
    rho_M += pvecback[pba->index_bg_rho_idm_dr];
  dy[pba->index_bi_D] = y[pba->index_bi_D_prime];
  dy[pba->index_bi_D_prime] = -a*H*y[pba->index_bi_D_prime] + 1.5*a*a*rho_M*y[pba->index_bi_D];

When the dcdm density is added in the $\rho_M$ factor, the difference between the two cases disappears which I think is very reasonable. I implemented the dcdm density in the above code in the following way,

rho_M = pvecback[pba->index_bg_rho_b];
  if (pba->has_cdm)
    rho_M += pvecback[pba->index_bg_rho_cdm];
  if (pba->has_idm_dr)
    rho_M += pvecback[pba->index_bg_rho_idm_dr];
  if (pba->has_dcdm == _TRUE_)
   rho_M += pvecback[pba->index_bg_rho_dcdm];
  dy[pba->index_bi_D] = y[pba->index_bi_D_prime];
  dy[pba->index_bi_D_prime] = -a*H*y[pba->index_bi_D_prime] + 1.5*a*a*rho_M*y[pba->index_bi_D];

The resultant plots after adding the dcdm density are the following, image image

Newer versions of CLASS also seem to miss the dcdm factor in the equation for growth factor so I want to ask if there is a reason to not include the dcdm density or is it a bug in CLASS?