Decaying Cold Dark Matter

[kabirakesh]

I ran CLASS code to reproduce the fig 1 of Audren et al (http://arxiv.org/pdf/1407.2418.pdf).

I used following values in .ini files:

for 1) stable ................ omega_cdm = 0.12038 Gamma_dcdm = 0.0 omega_k = 0. r=1.0 ................

and, for 2) Decaying ........... omega_cdm = 0.12038 Gamma_dcdm = 20.0 omega_k = 0. r=1.0 .............

But the cl vs l plot is exactly the same in both cases.

How can I produce this plot? Or plz correct if I am missing something.

[ThomasTram]

Hi,

that is because you have specified the abundance of cdm instead of dcdm (decaying cold dark matter). There are two ways of doing this. 1) Specify the initial abundance of dcdm using omega_ini_dcdm or Omega_ini_dcdm 2) Specify the final abundance of dcdm and decay radiation using omega_dcdmdr or Omega_dcdmdr

If you want to make plots comparing the effect of Gamma, it is better to use (1) since you will then keep the redshift of matter-radiation equality fixed.

[kabirakesh]

Hi Tram, thanks for the quick answer. But still there is some problem.

Below is the excerpt from my notebook, showing input parameters and outputs numbers

17:stable 100*theta_s = 1.042143 Omega_ini_dcdm = 0.12038 Gamma_dcdm = 0.0 omega_k = 0. r=0.0 ln10^{10}A_s = 3.0980 output = tCl, pCl modes = s, t format = CLASS

Writing output files in output/DCDM07_...

18: decaying 100*theta_s = 1.042143 Omega_ini_dcdm = 0.12038 Gamma_dcdm = 20.0 omega_k = 0. r=0. ln10^{10}A_s = 3.0980 output = tCl, pCl modes = s, t format = CLASS

k_output_values =1,10

l_max_scalars= 2500 lensing= no

output/DCDM08...

.... Although the plot between decaying and stable cdm is more or less similar to that in Paper fig 1, problem starts with Tensors.

20: Tensors 100*theta_s = 1.042143 Omega_ini_dcdm = 0.12038 Gamma_dcdm = 0.0 omega_k = 0. r=0.2 ln10^{10}A_s = 3.0980 output = tCl, pCl modes = s, t format = CLASS

output/DCDM10_...

observation: omega_k is not having any effect. let us change that to Omega_k

................. 23:open 100*theta_s = 1.042143 Omega_ini_dcdm = 0.12038 Gamma_dcdm = 0.0 Omega_k = 0.02 r=0.0 ln10^{10}A_s = 3.0980 output = tCl, pCl modes = s, t format = CLASS

k_output_values =1,10

l_max_scalars= 2500 lensing= no

output/DCDM14_...

ya it is confirmed that omega_k with any value is giving flat universe, now for closed universe ....... 24:closed 100*theta_s = 1.042143 Omega_ini_dcdm = 0.12038 Gamma_dcdm = 0.0 Omega_k = -0.02 r=0.0 ln10^{10}A_s = 3.0980 output = tCl, pCl modes = s, t format = CLASS

k_output_values =1,10

l_max_scalars= 2500 lensing= no

output/DCDM15_...

But plot (of stable, decaying and Tensors) did not match with figure 1 in the paper, as the tensor curve has to be in the middle, but it went above Decaying CDM (see the attached fig below). Also "open" is below the stable, and closed above the stable ---opposite from that in fig 1 of the paper.

I would also appreciate if you will comment why any value for "omega_k" was giving flat universe results.

gnuplot script for the figure (plz comment if script is not right to get the zoom): set logscale x 10 set xlabel 'l' set ylabel 'l(l+1)C_l / 2{/Symbol p}' set xr [2:]

set yr [10^-10:1e-9]

set format y '%.0t*10^{%T}'

set key right set title 'CLASS output' plot [2:50] [-23.:-22.7]'output/DCDM07_cl.dat' u 1:(log($2)) w l title 'stable', \ 'output/DCDM08_cl.dat' u 1:(log($2)) w l title 'decaying', \ 'output/DCDM10_cl.dat' u 1:(log($2)) w l title 'tensors',\ 'output/DCDM14_cl.dat' u 1:(log($2)) w l title 'open', \ 'output/DCDM15_cl.dat' u 1:(log($2)) w l title 'closed'

[ThomasTram]

Hi again,

I can see two problems: omega_k is not implemented, only Omega_k. If you had been running from the Python wrapper you would get the message:

Error in Class: Class did not read input parameter(s): omega_k

From command line you can turn this warning on by using

write warnings = yes

as explained in explanatory.ini. The second problem is that you are using Omega_ini_dcdm, but the plot in the paper is generated using omega_ini_dcdm = 0.12038. This is much closer to the current best fit DM abundance from Planck.

[kabirakesh]

Hi again,

And thanks for the answer. But I am getting error on using "omega_ini_dcdm = 0.12038"

ini file:

All the parameters value are set by default to the Planck +WP best fit

Only few parameters will be varied.

100*theta_s = 1.042143 omega_ini_dcdm = 0.12038 gamma_dcdm = 0.0 Omega_k = 0.0 r=0.0 ln10^{10}A_s = 3.0980 output = tCl, pCl modes = s, t format = CLASS

k_output_values =1,10

l_max_scalars= 2500 lensing= no

# background_verbose = 1 thermodynamics_verbose = 1 perturbations_verbose = 1 transfer_verbose = 1 primordial_verbose = 1 spectra_verbose = 1 nonlinear_verbose = 1 lensing_verbose = 1 output_verbose = 1

terminal output: Error running input_init_from_arguments =>input_init_from_arguments(L:179) :error in input_init(&fc, ppr, pba, pth, ppt, ptr, ppm, psp, pnl, ple, pop, errmsg); =>input_init(L:415) :error in fzero_Newton(input_try_unknown_parameters, x_inout, dxdF, unknown_parameters_size, 1e-4, 1e-6, &fzw, &fevals, errmsg); =>fzero_Newton(L:1105) :error in func(x_inout, x_size, param, F0, error_message); =>input_try_unknown_parameters(L:3345) :error in background_init(&pr,&ba); =>background_init(L:533) :condition ((pba->H0 < _H0SMALL)||(pba->H0 > _H0BIG)) is true; H0=8.38694e-05 out of bounds (0.00010007<H0<0.000333567)

[ThomasTram]

Hi again,

the problem is that you are not setting the CDM abundance to 0, so you are effectively considering a cosmology with twice as much dark matter than what is observed. You must add either omega_cdm = 0.001 or

omega_cdm = 0.0
gauge = Newtonian

(Synchronous gauge requires a small amount of cold dark matter.)

Regarding the error message, (in case someone finds this issue by Google) the problem is that combining decaying DM and theta_s transforms the initial value problem into a boundary value problem (BVP) with two values fixed today. Then the code uses the shooting method to transform the BVP into a system of nonlinear algebraic equations which are then solved by Newton's method. The last step is highly sensitive to the initial guess provided by the user in input_get_guess(), and this guess may be bad once the code is far from standard values. However, I do have a patch that enhances the robustness for BVPs, and the patched version could deal with your input file as it were. This patch will be part of the next public release.

Cheers. Thomas

[kabirakesh]

Hi again,

Regarding error message, I understand. But now stable and dcdm are giving same plot: 28:stable 100*theta_s = 1.042143 omega_cdm = 0.001 omega_ini_dcdm = 0.12038 gamma_dcdm = 0.0 Omegak = 0.0 r=0.0 output/DCDM16...

no error...

.... 29:decaying 100*theta_s = 1.042143 omega_cdm = 0.001 omega_ini_dcdm = 0.12038 gamma_dcdm = 20.0 Omega_k = 0.0 r=0.0

"open" was giving previous error, so "theta_s" was commented out for "open" run:

[ThomasTram]

You really need to start using write warnings = yes You have now changed your "Gamma_dcdm" into "gamma_dcdm", so "Gamma_dcdm" takes the default value of 0.0.

[AndreV84]

@ThomasTram Maybe you could add on the following too or shall a new issue be created? Thanks from colleague: I have the following parameters input :

params = {
'Omega_Lambda' : '0.0',
'Omega_fld' : '0.0',
A_planck = 1.0

params = {
'Omega_Lambda' : '0.0',
'Omega_fld' : '0.0',
'Omega_smg' : '-1.0',
'gravity_model' : 'brans dicke',
'lensing': 'yes',
'l_max_scalars': str(lmax),
'a_min_stability_test_smg' : '3e-13',
'output_background_smg': '10',
'write parameters': 'no',
'write background': 'no',
'write thermodynamics': 'no',
'input_verbose' : '0',
'background_verbose' : '0',
'output_verbose' : '0',
'thermodynamics_verbose' : '0',
'perturbations_verbose' : '0',
'spectra_verbose' : '0',
'omega_b' : '0.022032',
#'omega_cdm' : '1e-3',
'omega_cdm' : '1e-6',
'gauge' : 'Synchronous',
'YHe': '0.245',
'H0' : '67.4',
#'pert_ic_tolerance_smg': '1e-6', # Exemple d'ajustement
}

but I get the follwing error after few iterations :

hi_classy.CosmoComputationError:

Error in Class: background_init(L:746) :error in background_solve(ppr,pba);
=>background_solve(L:2282) :condition (pba->min_cs2_smg < -fabs(pba->cs2_safe_smg)) is true; Gradient instability for scalar field perturbations with minimum c_s^2=-301.804

new issue created: https://github.com/lesgourg/class_public/issues/576