Skip to content

Shooting

Jump to bottom
Kay Lehnert edited this page Apr 1, 2026 · 8 revisions

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

$\Omega_\phi$ is chosen to fill the energy budget, i.e. we need the corresponding $\rho_\phi$. The distance conjecture informs us that $\Delta\phi\lesssim 10$, i.e. we can assume $\rho_\phi\approxV(\phi)/3$ as a slow-roll approximation. This means our target potential is $V\approx3\Omega_\phi H_0^2\approx10^{-7}$. From this, we can deduce an initial guess for $c_1$:

$$c_1\approx \left(\frac{V}{g\left(\phi_{ini}\right)\right)^{\frac{1}{p}}$$

for a potential of the form $V=c_1^pg\left(\phi\right)$. Initial tests showed that using this approximation reduced the number of shooting steps by 25% on average, across the different potentials. We implemented this initial guess not in CLASS, yet directly in the Cobaya YAML file, which we create for the different potentials with the following code snippet in create_cobaya_yaml.py:

if not is_attractor:
            # V_target = 3 * Omega_scf * (H0/c)^2
            # h = 0.6736, Omega_scf ≈ 0.685, c = 299792.458 km/s
            _V_T = 3.0 * 0.685 * (67.36 / 299792.458) ** 2  # ≈ 1.037e-7

            if potential == "cosine":
                # V = c1*cos(c2*phi) = c1*cos(xi_ini)  [xi = c2*phi]
                scf_c1 = {
                    "derived": f"lambda xi_ini: {_V_T} / np.cos(xi_ini)"
                    f" if np.cos(xi_ini) > 0.01 else 1e-7",
                    "drop": True,
                    "latex": "c_1",
                }
            elif potential == "hyperbolic":
                # V = c1*[1-tanh(c2*phi)] = c1*[1-tanh(chi_ini)]
                # Guard chi_ini < 15 to avoid 1-tanh(x) losing all digits.
                scf_c1 = {
                    "derived": f"lambda chi_ini: {_V_T} / (1 - np.tanh(chi_ini))"
                    f" if chi_ini < 15 else 1e-8",
                    "drop": True,
                    "latex": "c_1",
                }
            elif potential == "pNG":
                # V = c1^4*(1+cos(phi/c2)) = c1^4*(1+cos(xi_ini))
                scf_c1 = {
                    "derived": f"lambda xi_ini: ({_V_T} / (1 + np.cos(xi_ini)))**0.25"
                    f" if (1 + np.cos(xi_ini)) > 0.01 else 1e-1",
                    "drop": True,
                    "latex": "c_1",
                }
            elif potential == "exponential":
                # V = c1*exp(-c2*phi) → c1 = V_T*exp(c2*phi)
                # Rewrite as multiplication to avoid 1/exp → 0 division.
                # Clamp exponent to [-500, 500] against overflow.
                scf_c1 = {
                    "derived": f"lambda psi_ini, scf_c2, cdm_c: {_V_T}"
                    f" * np.exp(np.clip(scf_c2 * psi_ini / cdm_c, -500, 500))"
                    f" if abs(cdm_c) > 1e-6 else 1e-7",
                    "drop": True,
                    "latex": "c_1",
                }
            elif potential in ("Bean", "BeanSingleWell"):
                # V = c1*[(c4-phi)^2+c2]*exp(-c3*phi)
                # phi = psi_ini/c3, so exp(-c3*phi) = exp(-psi_ini)
                # Clamp result to [1e-15, 1e5] against tiny/huge denominators.
                scf_c1 = {
                    "derived": f"lambda psi_ini, scf_c2, scf_c3, scf_c4:"
                    f" np.clip({_V_T}"
                    f" / (max(((scf_c4 - psi_ini/scf_c3)**2 + scf_c2), 1e-30)"
                    f" * np.exp(np.clip(-psi_ini, -500, 500))), 1e-15, 1e5)",
                    "drop": True,
                    "latex": "c_1",
                }
            elif potential == "BeanAdS":
                # V = c1*[(c4-phi)^2+c2]*exp(-c3*phi), phi sampled directly
                # Clamp result to [1e-15, 1e5]; clamp exp argument.
                scf_c1 = {
                    "derived": f"lambda scf_phi_ini, scf_c2, scf_c3, scf_c4:"
                    f" np.clip({_V_T}"
                    f" / (max(((scf_c4 - scf_phi_ini)**2 + scf_c2), 1e-30)"
                    f" * np.exp(np.clip(-scf_c3 * scf_phi_ini, -500, 500))),"
                    f" 1e-15, 1e5)",
                    "drop": True,
                    "latex": "c_1",
                }
            elif potential == "DoubleExp":
                # V = c1*(exp(-c2*phi)+c3*exp(-c4*phi)), phi sampled directly
                # Clamp exp arguments and result against overflow.
                scf_c1 = {
                    "derived": f"lambda scf_phi_ini, scf_c2, scf_c3, scf_c4:"
                    f" np.clip({_V_T}"
                    f" / (np.exp(np.clip(-scf_c2 * scf_phi_ini, -500, 500))"
                    f" + scf_c3 * np.exp(np.clip(-scf_c4 * scf_phi_ini, -500, 500))),"
                    f" 1e-15, 1e5)",
                    "drop": True,
                    "latex": "c_1",
                }
            elif potential == "power-law":
                # V = c1^(4-c2)*phi^c2  (c3=0)
                # c1 = (V_t / phi^c2)^(1/(4-c2))
                scf_c1 = {
                    "derived": f"lambda scf_phi_ini, scf_c2:"
                    f" ({_V_T} / scf_phi_ini**scf_c2)**(1.0/(4.0 - scf_c2))"
                    f" if abs(4 - scf_c2) > 0.01 and scf_phi_ini > 0 else 1e-2",
                    "drop": True,
                    "latex": "c_1",
                }
            elif potential == "SqE":
                # V ≈ c1^(c2+4)*phi^(-c2)  (exp(c1*phi^2) ≈ 1 for small c1)
                # c1 = (V_t * phi^c2)^(1/(c2+4))
                scf_c1 = {
                    "derived": f"lambda scf_phi_ini, scf_c2:"
                    f" ({_V_T} * scf_phi_ini**scf_c2)**(1.0/(scf_c2 + 4.0))"
                    f" if (scf_c2 + 4) > 0.01 and scf_phi_ini > 0 else 1e-2",
                    "drop": True,
                    "latex": "c_1",
                }

Logarithmic Shooting

The linear-space Newton solver fails if the step size is either too large (overshoots to $c1&lt;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;

Clone this wiki locally