-
Notifications
You must be signed in to change notification settings - Fork 0
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 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). 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
for a potential of the form 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",
}The linear-space Newton solver fails if the step size is either too large (overshoots to input_try_unknown_parameters. The Jacobian
Some potentials use
-
$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 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;If the initial guess is far away from the true value, overshooting can cause CLASS to crash. Algebraically, a value is then reached that is not only numerically unstable but also unphysical. To avoid this situation, we cap the maximum step size in input_find_root() to 1 decade, i.e. a factor of 10:
{
double max_dx = MAX(fabs(x1), 1.0);
if (fabs(dx) > max_dx)
dx = copysign(max_dx, dx);
}Since the rest of CLASS does not know about the log-shooting, the value is transformed back into the physical value in input_try_unknown_parameters():
pfzw = (struct fzerofun_workspace *)voidpfzw;
/** Read input parameters */
// This needs to be done with enough accuracy. A standard double has a relative
// precision of around 1e-16, so 1e-20 should be good enough for the shooting
for (i = 0; i < unknown_parameters_size; i++)
{
double val = unknown_parameter[i];
/* KBL: If this unknown parameter is in log-space (e.g. u = log10(c1)),
exponentiate to get the physical value before CLASS sees it */
if (pfzw->shooting_log_space != NULL &&
pfzw->shooting_log_space[i] == _TRUE_)
{
val = pow(10.0, val);
}
class_sprintf(pfzw->fc.value[pfzw->unknown_parameters_index[i]], "%.20e", val);
}The necessary flag (if it is a logarithmic parameter or not) is added to the fzerofun_workspace in input.h:
struct fzerofun_workspace
{
int *unknown_parameters_index;
struct file_content fc;
enum target_names *target_name;
double *target_value;
int target_size;
enum computation_stage required_computation_stage;
int *shooting_log_space; /**< KBL: array of size target_size; _TRUE_ if the
corresponding unknown parameter should be
exponentiated (10^u) before CLASS evaluation.
Used for log-space c1 shooting with the
hyperbolic SCF potential. */
};