Understanding decaying Dark matter in input.c and background.c

[bhaskardas13]

Hi,

I am working on a decaying cold dark matter model where the background equations are little different.

So, I need to understand the derivation of some equations in input.c and background.c.

I couldn’t understand the meaning of a_decay and derivation of the equation in input.c a_decay = pow(1+(gamma*gamma-1.)/Omega_M,-1./3.)

I am having trouble understanding this part of the code in input.c

case omega_ini_dcdm:
  Omega0_dcdmdr = 1./(ba.h*ba.h);
case Omega_ini_dcdm:
  /* This works since correspondence is Omega_ini_dcdm -> Omega_dcdmdr and
     omega_ini_dcdm -> omega_dcdmdr */
  Omega0_dcdmdr *=pfzw->target_value[index_guess];
  Omega_M = ba.Omega0_cdm+ba.Omega0_idm+Omega0_dcdmdr+ba.Omega0_b;
  gamma = ba.Gamma_dcdm/ba.H0;
  if (gamma < 1)
    a_decay = 1.0;
  else
    a_decay = pow(1+(gamma*gamma-1.)/Omega_M,-1./3.);
  xguess[index_guess] = pfzw->target_value[index_guess]*a_decay;
  dxdy[index_guess] = a_decay;
  if (gamma > 100)
    dxdy[index_guess] *= gamma/100;
  break;

In background.c the expression for the critical density fraction is given like this. I couldn’t understand how to derive this equation. I need to understand this so that I could derive it according to the model I am working on.

The part in background.c

if (pba->has_dr == _TRUE_) {
if (pba->has_dcdm == _TRUE_) {
  /**  - f is the critical density fraction of DR. The exact solution is:
   *
   * `f = -Omega_rad+pow(pow(Omega_rad,3./2.)+0.5*pow(a,6)*pvecback_integration[pba->index_bi_rho_dcdm]*pba->Gamma_dcdm/pow(pba->H0,3),2./3.);`
   *
   * but it is not numerically stable for very small f which is always the case.
   * Instead we use the Taylor expansion of this equation, which is equivalent to
   * ignoring f(a) in the Hubble rate.
   */
  f = 1./3.*pow(a,6)*pvecback_integration[pba->index_bi_rho_dcdm]*pba->Gamma_dcdm/pow(pba->H0,3)/sqrt(Omega_rad);
  pvecback_integration[pba->index_bi_rho_dr] = f*pba->H0*pba->H0/pow(a,4);

It would be very helpful if someone could explain or I could get any references where I can find the derivations.

Thank you very much

[ThomasTram]

Hi There is no reference, sorry! The first code blob is for finding a good starting guess for the shooting method. The idea is to find the decay-time, a_decay, and then assume instantaneous decay at this moment.

I think your second blob is the initial conditions for the decay radiation -- if you start early enough putting this to 0 should be OK, depending on your model. Otherwise you just assume radiation domination and derive it by solving your equations in the limit.

Cheers, Thomas

[bhaskardas13]

Hi Thank you for the response. I am sorry but I am not sure about the meaning of a_decay here. Does it mean the scale factor at the initial time when decay is starting? In my model the Gamma_dcdm is not constant in time so I need to understand the derivation of the expression for a_decay. Can you give a hint about how to get the expression for a_decay?

Similarly, can you give a hint about how do I get the expression for f-the critical density fraction of DR for my model... I mean by solving which equations.

Thank you very much Bhaskar Das

[ThomasTram]

a_decay is defined as the time when Gamma / H == 1, i.e. decay constant is similar to Hubble-rate. For the density of dr, you just solve the coupled system of background equations in radiation domination using e.g. Mathematica. If you cannot find an analytical solution, just use series expansions!