Shooting

To fill the energy budget constraints, a shooting method is applied that fills the energy budget with dark energy. All implemented dark energy potentials have an overall scaling factor $c_1$. To tell CLASS to use this factor for shooting, the input parameter scf_tuning_index = 0 is to be passed (which is also the default value, in case nothing is passed). $c_1$ is the varied by the shooting algorithm within CLASS. Therefore, it is not a free parameter that needs to be estimated through the MCMC. What is passed by the user as the value of $c_1$ is only used as the initial guess for this variable. The shooting itself is not able to cover several orders of magnitude, if we have no idea about the value of $c_1$. We implemented two separate approaches to deal with this issue.

Analytical Guess

Logarithmic Shooting

The linear-space Newton solver fails if the step size is either too large (overshoots to c1<0) or too small (stalls). We can solve this by shooting $u = \log_{10}(c_1)$ instead of $c_1$ directly. The transform $u\rightarrow c_1 = 10^u$ is applied in input_try_unknown_parameters. The Jacobian $du/d\Omega = 1/(p\Omega\log(10))$ is scale-free ($c_1$ cancels), ensuring uniform convergence across decades.

Some potentials use $c_1$ to the power of something, such that $V\propto c_1^p$, with a potential-dependent exponent $p$:

• $p=1$: cosine (2), hyperbolic (3), exponential (6), Bean (8), double exponential (9)
• $p=4$: pseudo-Nambu–Goldstone (4)
• $p=4-c_2$: power-law (1), only when $c_3=0$
• $p=4+c_2$: inverse power-law (5)
• $p=(c_2+4)+c_1\phi^2$: squared exponential (7)

If $c_1$ is not the shooting parameter, $p=0$, or $c_1&lt;0$, the linear-space shooting algorithm is used as a fallback. This is implemented in the input_get_guess() function in input.c as follows:

case Omega_scf:
      /* *
       * KBL: Log-space shooting for the scalar field amplitude c1.
       *
       * When shooting on scf_tuning_index == 0 (the amplitude parameter c1),
       * we work in u = log10(c1) instead of c1 directly.  This handles the
       * multi-decade dynamic range of c1 that arises because Omega_scf ∝ c1^p,
       * where p depends on the potential:
       *
       *   potential   V(phi)                                        p
       *   ─────────   ──────                                        ──
       *   1 power-law c1^(4-c2)*phi^c2 + c3                        4-c2  (only if c3==0)
       *   2 cosine    c1*cos(c2*phi)                                1
       *   3 hyperbolic c1*[1-tanh(c2*phi)]                          1
       *   4 pNG       c1^4*[1+cos(phi/c2)]                          4
       *   5 iPL       c1^(4+c2)*phi^(-c2)                           4+c2
       *   6 exponential c1*exp(-c2*phi)                              1
       *   7 SqE       c1^(c2+4)*phi^(-c2)*exp(c1*phi^2)             (c2+4)+c1*phi^2
       *   8 Bean      c1*[(c4-phi)^2+c2]*exp(-c3*phi)               1
       *   9 DoubleExp c1*(exp(-c2*phi)+c3*exp(-c4*phi))              1
       *
       * The Jacobian du/dΩ = 1/(p * Ω_scf * ln10).
       *
       * For SqE, p_eff = d(ln V)/d(ln c1) = (c2+4) + c1*phi^2 is not
       * constant in c1, but evaluated at the initial guess it gives a
       * good enough Jacobian for the bracket search.  Monotonicity of
       * V(c1) guarantees a unique root.
       *
       * Conditions for log-space: scf_tuning_index == 0, c1 > 0, p ≠ 0,
       * and no additive constant breaking the proportionality V ∝ c1^p.
       * Falls back to linear-space (legacy) whenever these are not met.
       */
      if (ba.scf_tuning_index == 0)
      {
        double c1_guess = ba.scf_parameters[0];
        double c1_power = 0.0; /* exponent p in Omega_scf ∝ c1^p; 0 = cannot use log-space */

        switch (ba.scf_potential)
        {
        case 1: /* power-law: V ∝ c1^(4-c2) only if c3 == 0 */
          if (ba.scf_parameters[2] == 0.0)
            c1_power = 4.0 - ba.scf_parameters[1];
          break;
        case 2: /* cosine */
        case 3: /* hyperbolic */
        case 6: /* exponential */
        case 8: /* Bean */
        case 9: /* DoubleExp */
          c1_power = 1.0;
          break;
        case 4: /* pNG: V ∝ c1^4 */
          c1_power = 4.0;
          break;
        case 5: /* iPL: V ∝ c1^(4+c2) */
          c1_power = 4.0 + ba.scf_parameters[1];
          break;
        case 7: /* SqE: V = c1^(c2+4)*phi^(-c2)*exp(c1*phi^2)
                 * p_eff = d(ln V)/d(ln c1) = (c2+4) + c1*phi_ini^2
                 * Evaluated at the initial guess for c1. */
        {
          double phi_ini_sqe = ba.scf_parameters[ba.scf_parameters_size - 2];
          c1_power = (4.0 + ba.scf_parameters[1]) + c1_guess * phi_ini_sqe * phi_ini_sqe;
          /* Safety: if c2 ≈ -4 and c1*phi^2 ≈ 0, p_eff can be too small
             for a useful Jacobian.  Fall back to linear-space in that case. */
          if (c1_power < 0.01)
            c1_power = 0.0;
          break;
        }
        default:
          c1_power = 0.0;
          break;
        }

        if (c1_power != 0.0 && c1_guess > 0.0)
        {
          /* Log-space shooting: u = log10(c1), du/dΩ = 1/(p * Ω * ln10) */
          xguess[index_guess] = log10(c1_guess);
          dxdy[index_guess] = 1.0 / (c1_power * ba.Omega0_scf * _M_LN10_);
          pfzw->shooting_log_space[index_guess] = _TRUE_;
        }
        else
        {
          /* Fallback to linear-space (c1 ≤ 0, or p == 0, or unsupported potential) */
          xguess[index_guess] = c1_guess;
          dxdy[index_guess] = c1_guess / ba.Omega0_scf;
        }
      }
      else
      {
        /* Non-c1 tuning parameter: always linear-space */
        xguess[index_guess] = ba.scf_parameters[ba.scf_tuning_index];
        dxdy[index_guess] = ba.scf_parameters[ba.scf_tuning_index] / ba.Omega0_scf;
      }
      break;