Perturbation-level predictions of power-law f(Q) gravity from simulation-based inference

We present a comprehensive analysis of power-law symmetric teleparallel gravity, [Formula: see text], combining simulation-based parameter inference with perturbation-level predictions and multi-messenger forecasts. The model generalizes [Formula: see text]CDM through a single index [Formula: see text], with [Formula: see text] recovering standard cosmology exactly. By imposing the boundary condition [Formula: see text], we analytically eliminate the coupling constant [Formula: see text], reducing the model to three free parameters [Formula: see text] governed by an implicit Friedmann equation solved numerically via Brent’s method. We employ Conditional Likelihood-Free Inference (CoLFI), which trains neural density estimators on [Formula: see text] forward simulations per iteration, circumventing the intractable analytic likelihood inherent to this implicit equation. Using Pantheon+SH0ES supernovae, cosmic chronometers, Union 2.1, and DES-Y5 data, individually and in six combinations, the joint Pantheon+OHD analysis yields the tightest constraint: [Formula: see text], representing a [Formula: see text] departure from general relativity. The sign of the best-fit [Formula: see text] depends on the dataset combination, reflecting a residual [Formula: see text]–[Formula: see text] degeneracy that only multi-probe analyses can break. From the constrained posterior, we derive four perturbation-level predictions without requiring a perturbation-level likelihood: (i) the effective gravitational coupling [Formula: see text], which is enhanced at low redshift and drives excess structure growth; (ii) the gravitational slip [Formula: see text], departing from unity most strongly at [Formula: see text]; (iii) the growth rate [Formula: see text], consistent with BOSS/eBOSS data but distinguishable from [Formula: see text]CDM by Euclid; and (iv) the gravitational-wave luminosity-distance ratio [Formula: see text] at [Formula: see text], a 20% enhancement testable by LISA and the Einstein Telescope. Mapping the posterior onto the CPL parameterization yields [Formula: see text], placing the model in the phantom regime and geometrically disjoint from the DESI DR1 best-fit [Formula: see text]. Model comparison via AIC, BIC, and the Savage–Dickey Bayes factor indicates that current data accommodate a departure from GR but do not yet compel it.