Stability and precision

Stability, ICs and precision

hi_class checks every model for ghost and gradient instabilities, sets consistent initial conditions for the scalar field perturbations, and can evolve them quasi-statically when they oscillate much faster than the Hubble rate. This page documents those mechanisms and the precision parameters that control them (sections 3–7 of hi_class.ini).

Stability tests

After solving the background, hi_class scans the whole evolution and rejects the model (with an error naming the failed condition and the scale factor where it happened) if any of the four healthy-perturbations conditions is violated:

Condition	Instability	Quantity tested
D > 0	scalar ghost	kinetic term D = α_K + (3/2)α_B²
c_s² > 0	scalar gradient instability	scalar sound speed squared
M*² > 0	tensor ghost	effective Planck mass squared
c_T² > 0	tensor gradient instability	tensor speed squared, c_T² = 1 + α_T

All four quantities are printed in background.dat at output_background_smg ≥ 1 — see Output files. The tests only run once the model has been tuned to the requested Omega_smg (so intermediate shooting iterations are allowed to wander through unstable regions), and D is allowed to reach exactly zero, which can happen by underflow in models whose kineticity becomes arbitrarily small (e.g. quintessence with tiny initial φ′).

Three sets of parameters modify the tests:

• skip_stability_tests_smg (default no): set to yes to disable all four tests. Use with care — the perturbation equations are still integrated, and an actual instability will usually show up as runaway perturbations or solver failures.
• Thresholds cs2_safe_smg, ct2_safe_smg, D_safe_smg, M2_safe_smg (all default 0.): each test is relaxed to quantity > −|threshold|, i.e. slightly negative values are tolerated. This is meant to absorb numerical noise around zero, not genuine instabilities. Note the code takes the absolute value, so the sign of the threshold does not matter.
• a_min_stability_test_smg (default 0): the scan ignores everything at a < a_min. Some shipped models set it because the offending values are numerical artifacts of the deep radiation era, where the scalar is dynamically irrelevant but cancellations make the stability quantities noisy: gravity_models/brans_dicke.ini and the cccg_exp*.ini files use 1e-6, the quintessence models use 1e-10.

Background evolution options

Section 4 of hi_class.ini. For covariant theories the Friedmann equation is generally cubic in H, with multiple branches. Instead of solving the constraint algebraically at every step, hi_class integrates H dynamically and uses the constraint as a friction term that damps violations:

• hubble_evolution (default y): integrate H from its equation of motion. The constraint is only solved once, to set the initial H. Switching this off makes the code use the constraint equation directly, which will not work for most Lagrangian-based theories (branch ambiguities).
• hubble_friction (default 3.): coefficient of the constraint-damping term added to H′. Larger values pull the trajectory back to the constraint surface more aggressively; the sign is chosen automatically so the term is always restoring.

Perturbation initial conditions

At early times the scalar field perturbation is attracted to a solution driven by the metric perturbations of the radiation era. pert_initial_conditions_smg selects how the scalar is initialized when a mode starts fully dynamical (section 5 of hi_class.ini). The accepted values are:

• ext_field_attr (default) — external field attractor: assumes Ω_smg and all α's are small at early times, so the scalar does not gravitate; it simply follows the superhorizon metric potential h ∝ τ² created by radiation. Before initialization the code checks for a tachyonic instability of the scalar in radiation domination and aborts if it finds one (the ICs would not be trustworthy).
• gravitating_attr — gravitating attractor: allows the scalar to gravitate at early times (early modified gravity), solving the coupled radiation–scalar system with the assumption of constant Ω_smg and α's during radiation domination. The superhorizon evolution of all other species is modified accordingly (η and the neutrino velocity/shear ICs are recomputed). The code verifies that the standard adiabatic mode h ∝ τ² still dominates; if a faster-growing mode exists, curvature is not conserved on super-horizon scales, the connection to the inflationary amplitude is lost, and the run aborts. Isocurvature initial conditions are not implemented for this choice.
• single_clock — single-clock initial conditions set with respect to the photons: x_smg = δφ/φ̇ equals the common time shift δτ, implemented as x = −δ_γ/(4H) with x′ = 0.
• kin_only — solves the field equation keeping only the kineticity, x″ = 2(aH)²x, giving x ∝ (kτ)².
• zero — start with x = 0 and x′ = 0. Useful to check how much the results depend on the scalar's ICs (for attractor solutions the late-time observables should barely change).

Note that zero and kin_only appear only for completeness/debugging; the physical choices are the two attractors, described in detail in arXiv:1909.01828.

The associated tolerance parameters (defaults in parentheses; all can be set in the .ini):

Parameter	Default	Role
pert_ic_ini_z_ref_smg	1e10	redshift at which the initial-conditions stability tests above are performed
pert_ic_tolerance_smg	2e-2	tolerated deviation from the standard growth h ∝ τ² (n = 2) at the reference redshift; a negative value disables the test
pert_ic_regulator_smg	1e-15	minimum size of denominators in the IC expressions, regulated to prevent infinities; negative disables the regulator
pert_qs_ic_tolerance_test_smg	1. (the shipped hi_class.ini uses 10)	for modes that start quasi-static: maximum contribution of the scalar to superhorizon curvature non-conservation, as a fraction of the standard radiation contribution

Quasi-static approximation

When the scalar's effective mass is large, its perturbations oscillate on timescales far shorter than the Hubble time; integrating those oscillations is expensive and pointless for observables. hi_class can then replace the scalar's equation of motion with its quasi-static (QS) algebraic solution. method_qs_smg (section 6 of hi_class.ini) selects the strategy:

Value	Behavior
fully_dynamic (default; abbreviation fd)	always integrate the full scalar dynamics
quasi_static (qs)	force the QS solution at all times
automatic (a)	let the code decide, mode by mode and time by time (recommended for stiff models)
fully_dynamic_debug (fdd)	integrate fully, but also evaluate the automatic decision machinery (for comparison/debugging)
quasi_static_debug (qsd)	force QS, but also evaluate the decision machinery

What automatic does. For each k-mode the code samples, over the whole integration interval, an effective mass² of the scalar oscillations (in units of aH, including the c_s²k² pressure support), its time derivative, a radiation trigger, a friction coefficient and the decay rate ("slope") of the oscillation envelope; these are the mass2_qs, mass2_qs_p, rad2_qs, friction_qs, slope_qs columns of the perturbations_k*.dat files (Output files). A time is proposed as quasi-static when both

• mass² > trigger_mass_qs_smg² — the oscillations are fast compared to the expansion, and
• rad² > trigger_rad_qs_smg² — radiation perturbations no longer drive the scalar strongly,

and only at redshifts above z_fd_qs_smg (below it the fully dynamic evolution is always used). The proposed scheme is then compacted to the times where the status actually switches, and each dynamic→QS transition is delayed until the oscillation amplitude, estimated from the slope, has decayed by the factor eps_s_qs_smg, so that the mode lands on the QS solution rather than jumping onto it mid-oscillation. Finally, the scheme is fitted to the transition sequences the integrator implements. The code also checks that at a very early reference time (a_ini_test_qs_smg = 1e-14) all k-modes share the same status, and anticipates the start of integration for a mode if needed.

The knobs (defaults from include/precisions.h, all settable in the .ini):

Parameter	Default	Role
z_fd_qs_smg	10.	minimum redshift for QS: below it perturbations are always evolved fully dynamically
trigger_mass_qs_smg	1e3	mass threshold (in units of aH) above which QS may switch on
trigger_rad_qs_smg	1e3	radiation trigger threshold
eps_s_qs_smg	0.01	required decay of the oscillating solution before switching to QS
n_min_qs_smg / n_max_qs_smg	1e2 / 1e4	minimum/maximum number of sampling points used to scan the triggers in time (sampling adapts to how fast the mass varies)

Raising the triggers or lowering eps_s_qs_smg makes the scheme more conservative (more fully-dynamic evolution); the defaults were validated in arXiv:1909.01828.

Precision parameters

Section 7 of hi_class.ini, plus the smg-related adjustments hi_class makes to CLASS defaults.

• kineticity_safe_smg (default 0.): adds a constant to α_K. Historically used to damp the impact of rapid early-time oscillations that slow down the integrator; it should not affect observable scales, but use at your own risk — since v2 the recommended cure is the quasi-static approximation (method_qs_smg = automatic) instead.
• tol_initial_Omega_r (default 1e-4): standard CLASS parameter setting how close to total radiation domination the background integration must start. Early-dark-energy-like models where the scalar is non-negligible at early times can fail this check; loosening the tolerance lets the run proceed.
• tol_background_integration (default 1e-10): relative tolerance of the background ODE integration. hi_class tightens this from the CLASS default of 1e-2 — the shooting method and the stability scan need an accurately integrated scalar field. If a covariant model's background integration is unusually slow, this (together with background_Nloga) is the knob to inspect.
• background_Nloga: number of stored background table points. Upstream CLASS uses 40000; for smg runs hi_class automatically resets it to 3000 unless you set it explicitly, because in hi_class the background is recomputed at every shooting iteration and the larger table costs roughly an order of magnitude in run time for no measurable accuracy gain. Pure ΛCDM runs (Ω_smg = 0) keep the upstream default and remain directly comparable to CLASS.
• perturbations_sampling_stepsize: automatically capped at 0.05 for smg runs (CLASS default 0.1), otherwise the ISW effect is not sampled accurately in modified-gravity models.
• In addition, the shipped hi_class.ini sets a few CLASS parameters to safer-than-default values for modified gravity — start_small_k_at_tau_c_over_tau_h = 1e-4 (CLASS default 0.0015), start_large_k_at_tau_h_over_tau_k = 1e-4 (default 0.07), so that modes start being integrated earlier, and a denser multipole sampling l_logstep = 1.045, l_linstep = 50 (defaults 1.12, 40). These are set in the example input file, not in the code: if you write your own .ini from scratch and need highly accurate C_ℓ's, copy them over.

For the accuracy of the ~k > 10 Mpc⁻¹ modes, see also get_h_from_trace (discussed in hi_class.ini and printed as a warning at high input_verbose): hi_class can obtain the metric perturbation h′ either from the Einstein 00 constraint (default) or from the trace of the ij equation, which avoids a gauge-dependent singularity when a model crosses α_B = 2.