Coupling

The general function for a coupling between dark matter and dark energy is

$$\Xi=\frac{3H}{\left(\rho_\textnormal{DM}+\rho_\textnormal{DE}\right)^{\left(l-1\right)}} \left(\mathfrak{a}\rho_\textnormal{DE}^{l-k}\cdot\rho_\textnormal{DM}^k+\mathfrak{b}\rho_\textnormal{DM}^{l-k}\cdot\rho_\textnormal{DE}^k\right) +\dot{\phi}\left(\frac{\mathfrak{c}\rho_\textnormal{DM}}{M}+\mathfrak{d}\rho_\textnormal{DM}\prime\right).$$

This can be added to the scalar field potential to get an effective potential $V_\textnormal{eff}\prime=V\prime+\Xi/\dot{\phi}$. To achieve this, we implemented the general coupling function into background.c as

// KBL: General coupling function Xi/phi' that can be added to V' to get V'_eff including a general coupling between DE and DM
double coupling_scf(
    struct background *pba,
    double rho_cdm_prime, // derivative of dark matter density wrt scalar field respectively H
    double *pvecback)
{
  /*
   * The scalar field parameters are given as a list of parameters with the following format:
   * [0]     = c_1
   * [1]     = c_2
   * [2]     = c_3
   * [3]     = c_4
   * [4]     = q_1 (coupling)
   * [5]     = q_2
   * [6]     = q_3
   * [7]     = q_4
   * [8]     = coupling_exponent_1
   * [9]     = coupling_exponent_2
   * [10]    = phi_ini
   * [11]    = phi_prime_ini
   *
   * This function returns Xi/phi' right away, and can be added to V' without dividing by phi' again.
   */
  double q1 = pba->scf_parameters[pba->scf_parameters_size - 8];
  double q2 = pba->scf_parameters[pba->scf_parameters_size - 7];
  double q3 = pba->scf_parameters[pba->scf_parameters_size - 6];
  double q4 = pba->scf_parameters[pba->scf_parameters_size - 5];

  // Early return if no coupling is present to avoid numerical issues
  if (q1 == 0.0 && q2 == 0.0 && q3 == 0.0 && q4 == 0.0)
  {
    return 0.0;
  }

  // Cache frequently accessed array elements in locals
  double rho_cdm = pvecback[pba->index_bg_rho_cdm];
  double rho_scf = pvecback[pba->index_bg_rho_scf];

  // Direct coupling terms (always cheap, no transcendentals)
  double direct = q3 * rho_cdm + q4 * rho_cdm_prime;

  // If only direct coupling, skip all expensive pow/exp/log computations
  if (q1 == 0.0 && q2 == 0.0)
  {
    return direct;
  }

  // Guard against division by zero or near-zero phi_prime
  // Physical justification: if phi' is near zero, phi is not changing, and so its energy density
  // should not change by its own dynamics, only by the coupling to DM density evolution.
  double phi_prime = pvecback[pba->index_bg_phi_prime_scf];
  if (fabs(phi_prime) < 1e-20)
  {
    if (pba->background_verbose > 2)
    {
      printf("Warning: phi' is very small (%e). Using only the direct coupling terms q3=%e and q4=%e.\n",
             phi_prime, q3 * rho_cdm, q4 * rho_cdm_prime);
    }
    return direct;
  }

  double exp1 = pba->scf_parameters[pba->scf_parameters_size - 4];
  double exp2 = pba->scf_parameters[pba->scf_parameters_size - 3];

  // Precompute shared exponent differences
  double e1m1 = exp1 - 1.0;
  double e1me2 = exp1 - exp2;

  // Replace 5 independent pow() calls with 3 log() + 5 exp().
  // Each pow(base, y) is internally exp(y*log(base)), so sharing the 3 logarithms
  // across 5 bases saves 2 transcendental evaluations.
  double ln_sum = log(rho_cdm + rho_scf);
  double ln_scf = log(rho_scf);
  double ln_cdm = log(rho_cdm);

  double inv_sum_e1m1 = exp(-e1m1 * ln_sum); // 1 / (rho_sum)^(exp1-1)
  double scf_e1me2 = exp(e1me2 * ln_scf);    // rho_scf^(exp1-exp2)
  double cdm_e2 = exp(exp2 * ln_cdm);        // rho_cdm^exp2
  double cdm_e1me2 = exp(e1me2 * ln_cdm);    // rho_cdm^(exp1-exp2)
  double scf_e2 = exp(exp2 * ln_scf);        // rho_scf^exp2

  // Common prefactor: 3 * H / rho_sum^(exp1-1) / phi'
  double common = 3.0 * pvecback[pba->index_bg_H] * inv_sum_e1m1 / phi_prime;

  double term1 = q1 * scf_e1me2 * cdm_e2;
  double term2 = q2 * cdm_e1me2 * scf_e2;
  double indirect = common * (term1 + term2);
  double result = indirect + direct;

  if (pba->background_verbose > 6)
  {
    printf("Coupling term Xi/phi': %e\nComponents: %e, %e, %e, %e\n",
           result,
           common * term1,
           common * term2,
           q3 * rho_cdm,
           q4 * rho_cdm_prime);
  }

  return result;
}

The values of $H$, as well as the DM and scalar field densities $\rho$ are internally computed. The coupling parameters $\mathfrak{a}$, $\mathfrak{b}$, $\mathfrak{c}$, and $\mathfrak{d}$ are user-defined numbers that are read from the *.ini-file respectively set through the python wrapper. $^1$ If all coupling factors are set to zero, the coupling function returns $0$ right away, to avoid numerical instabilities. If the coupling is present, an issue might arise from very small values of $\phi^\prime$, as this is in the denominator and the coupling contribution would numerically explode. However, since in this case, the scalar field itself does hardly change, its dynamics should not be strongly influenced by its almost zero change. We circumvent this issue by neglecting the scalar field dynamics for the coupling if those are numerically tiny, and only consider the contribution to the coupling stemming from dark matter itself.

If the DM mass depends on the scalar field, there is an additional coupling term depending on the derivative of the DM density with respect to the scalar field, $\rho\prime$. For our distance conjecture–inspired DM mass function

$$m\left(\phi\right)=\frac{m_0}{2}\left[1-\tanh\left(\mathfrak{c}\phi\right)\right],$$

the derivative is given by

$$m\prime=\frac{-m_0\mathfrak{c}}{1+\cosh\left(2\mathfrak{c}\phi\right)}.$$

Since $m$ is not part of the CLASS code, we express the derivative in terms of the DM density by using $\rho=mn$ with $n$ the number density that will cancel out, and expressing $\rho\prime$ as

$$\frac{\rho\prime}{\rho}\rho=-c\left(1+\tanh\left(c\phi\right)\right)\rho.$$

This is then implemented in background.c:

// Derivative of rho_cdm wrt scalar field phi
double rho_cdm_prime(
    struct background *pba,
    double phi,
    double *pvecback)
{
  if (pba->model_cdm == 2) // Interacting DM model
  {
    return -pba->cdm_c * (1 + tanh(pba->cdm_c * phi)) * pvecback[pba->index_bg_rho_cdm];
  }
  // else if (pba->model_cdm == 1) // Hubbleian DM model: rho_cdm_prime is wrt the scalar field, since it's for the coupling. d rho / d H is not coupled to the scalar field.
  // {
  //   return 0.0; // pvecback[pba->index_bg_rho_cdm] / pvecback[pba->index_bg_H] * pvecback[pba->index_bg_H_prime];
  // }
  else
    return 0.;
}

For performance reasons, we deliberately choose the tanh() function, which is numerically faster than an expansion into exponentials.

---

$^1$: In the code we use q_1–q_4 for the coupling constants, as we use $q_1$ to $q_4$ in a Fraktur typeface in our thesis. However, the MathJax rendering for subscripts is hard to read, and the Fraktur typeface is not supported either. Therefore, we use $a-d$ on this page for the coupling constants.