<figure><img src="image_57.png" /><figcaption>Fig. 50. A flowchart summarizing the main differences between different samplers.</figcaption></figure>

### 3.1.1. Markov chain Monte Carlo sampling algorithms

MCMC algorithms explore the parameter space through a proposal-acceptance mechanism, building a *chain* of states, or points in parameter space, from which a (decorrelated) fair sample from the posterior can be extracted. Robust convergence tests (usually based on the stability of a single Markov chain or similarity between parallel ones), are crucial to avoid biased inferences and inaccurate error estimates, guaranteeing that the MCMC samples reliably represent the target posterior distribution.

The family of MCMC algorithms encompasses a large number of different approaches. They can be classified under different criteria: the strategy followed to propose new points, the acceptance-rejection test performed on the proposals, the distribution from which they effectively sample (not always the target distribution), whether they pre-condition the parameter space, and how.

In the simplest approach, at each step $\theta$, a new state $\theta'$ is drawn from a proposal density $q(\theta'|\theta)$ and accepted with probability

$$ \alpha = \min \left( 1, \frac{\phi(\theta')q(\theta|\theta')}{\phi(\theta)q(\theta'|\theta)} \right), \qquad (3.2) $$

where $\phi(\theta)$ is the target distribution, i.e., the posterior in a Bayesian framework. This test is called the *Metropolis-Hastings* (MH) criterion. Rejected proposals increase the weight of the current state within the chain. A typical choice for the proposal distribution is a multivariate normal $\theta' \sim N(\theta, \Sigma_T)$, with a variable dispersion covariance $\Sigma_T$ tailored to parameter values at each step.

While unsophisticated, this method can produce accurate inference for simple distributions in large dimensions, but its efficiency depends critically on the choice of an appropriate proposal distribution, and otherwise deteriorates quickly for highly correlated or complex distributions.

Efficiency can be increased by using the current chain (or parallel chains) to compute an affine transformation that would decorrelate the posterior locally (or globally if sufficiently Gaussian). This affine transformation can be chosen so that proposals involve the re-calculation of only parts of the likelihood pipeline [1915] (as implemented in MontePython [1916,1917] and Cobaya [1918]). Parallelization, together with these strategies, can make MCMC well-suited for computationally expensive likelihood evaluations.

For more complicated posterior structures, such as curving degeneracies that would need a local transformation to be decorrelated, ensemble methods such as emcee<fnref n="17" /> [1919] can use the information on the position of multiple walkers to automatically adapt proposals. More recently, the authors of pocoMC<fnref n="18" /> [1920] have proposed the use of Sequential Monte Carlo (i.e., sequentially sampling from a family of distributions that interpolate between the prior and the posterior) while learning a generalized transformation that decorrelates arbitrary parameter spaces, using Normalizing Flows [1921]. pocoMC can compute model evidences along MC samples.

*Slice sampling* [1922] is an alternative adaptive approach to propose new Markov chain states, that automatically adjusts step sizes to the local shape of the distribution. For each parameter update, it samples uniformly from a randomly-oriented “slice” under the PDF curve, $S_y = \{\theta : y < \phi(\theta)\}$, where y is uniformly drawn from $[0, \phi(\theta_{current})]$. A proposal $\theta'$ is then uniformly sampled from this slice through a stepping-out and shrinking procedure. This adaptivity makes slice sampling particularly robust, requiring minimal tuning while maintaining

<footnote n="17">¹⁷ https://github.com/dfm/emcee</footnote>

<footnote n="18">¹⁸ https://github.com/minaskar/pocome</footnote>
good mixing properties. The algorithm naturally accommodates different scales in different regions of the parameter space, though its efficiency can decrease in higher dimensions. A modern implementation combining slice sampling with an ensemble approach, with applications in cosmology and astronomy, is zeus<fnref n="19" /> [1923]. Slice sampling has also been used to increase the efficiency and robustness of nested samplers (see below).

### 3.1.2. Hamiltonian Monte Carlo

Hamiltonian Monte Carlo (HMC), also known as hybrid Monte Carlo, is an MCMC approach that incorporates gradient information through an augmented parameter space with momentum variables. The system evolves according to Hamiltonian dynamics

$$H(\theta, p) = -\ln \phi(\theta) + \frac{1}{2} p^T M^{-1} p, \qquad (3.3)$$

where M is a mass matrix. In this approach, proposals for accepted/rejected new states are obtained by letting the current state evolve according to an energy-conserving trajectory for a set amount of time — in practice some number of steps of corresponding of small but finite duration. By following continuous trajectories determined by local dynamics, rather than diffusive random steps, HMC significantly reduces random-walk behavior, improving sampling efficiency for posterior with complicated geometries.

A key challenge in HMC is selecting an appropriate integration time for the dynamic evolution of the steps, i.e., the number and duration of finite steps. The No-U-Turn Sampler (NUTS) [1924] refines this approach by adaptively determining the trajectory length, ensuring efficient exploration without manual tuning. By dynamically terminating paths upon noticing signs of retracing, NUTS avoids redundant computations while preserving detailed balance, thus enhancing the efficiency, particularly in high-dimensional spaces.

In order to exploit the advantages of HMC, cosmological codes and likelihoods need to be able to compute derivatives with respect to input parameters and intermediate observables, respectively. The most common approach is to build these codes using automatically-differentiable numerical frameworks that are able to produce gradients for arbitrary inputs. The most capable and popular one at the time of writing this paper is the Python library JAX [1925]. There is a growing body of JAX-based cosmological and astronomical libraries (JAX-COSMO<fnref n="20" /> [1926], DISCO-DJ<fnref n="21" /> [1927]) and likelihoods (candl<fnref n="22" /> [1928]), as well as machine-learning cosmological emulators (CosmoPower-JAX<fnref n="23" /> [1929]). Another popular programming framework that allows for automatic-differentiation is julia,<fnref n="24" /> that has been employed for example in LimberJack.jl<fnref n="25" /> [1930], Capse.jl<fnref n="26" /> [1931] and Effort.jl<fnref n="27" /> [1932].

An implementation of NUTS that can be interfaced with JAX-based codes is numpyro<fnref n="28" /> [1933]. Using numpyro and diverse JAX-based cosmological libraries and emulators, recent works [1934,1935] have reported efficiency gains of up to orders of magnitude, especially for parameter spaces of large dimensionality that are expected when performing inference on data from next-generation surveys.

### 3.1.3. Nested sampling

Nested sampling [1936] transforms the evidence calculation into a one-dimensional integral over prior mass X accumulated up to a given

likelihood value

$$Z = \int_{0}^{1} \mathcal{L}(X) dX \quad \text{with} \quad X(\lambda) = \int_{\mathcal{L}(\theta) > \lambda} \pi(\theta) d\theta, \qquad (3.4)$$

where the likelihood is evaluated such that $\mathcal{L}(X(\lambda)) = \lambda$ for some likelihood-cutoff value $\lambda$. Numerically, this integral is computed as a weighted sum $Z \sim \sum_{i=1}^{m} w_i L_i$, with $w_i \sim \Delta X$. The algorithm maintains a set of live points updated at every step: the one with the lowest likelihood is dropped (and referred to as a dead point), and a new one is sampled restricted to having a likelihood larger than that of the recently-dead point. The prior volume shrinkage $\Delta X$ is computed probabilistically and assigned the likelihood value of the dead point. In addition to the evidence calculation, the pool of dead points may be used to construct a fair sample from the posterior.

The procedure to sample a point from the prior subject to a minimum likelihood, i.e., within an iso-likelihood contour, is the main hurdle of the algorithm. Some implementations, such as Multi-Nest<fnref n="29" /> [293], construct a minimal ellipsoid containing the set of live points, and sample uniformly from an enlarged version of it. Ultra-Nest<fnref n="30" /> [1937], implements a different sampling approach based on ellipsoids centered around the live points themselves, and includes more robust uncertainty estimation for the integrated evidence. Others, like PolyChord<fnref n="31" /> [1938], propose new points by running a short affine-invariant MCMC chain (in particular using slice-sampling) from one of the current live points. MCMC-based approaches have better dimensionality scaling, but tend to be slower for simpler problems. Other modern implementations like dynasty<fnref n="32" /> [1939] dynamically vary the number of live points to improve efficiency. All these approaches allow for efficient sampling of multi-modal distributions by using clustering algorithms to separate and evolve different subsets of live points in parallel.

In general, nested samplers are preferred when the geometry of the posterior distribution is complicated or exhibits multiple peaks. Due to its ability to compute the Bayesian evidence, nested sampling is also particularly valuable for model selection problems. Nested sampling can however be computationally very expensive for highly dimensional parameter spaces, but it parallelizes very efficiently up to the number of live points.

There have been in recent years various attempts at accelerating nested sampling using ML, focusing on efficiently sampling within iso-likelihood contours. nautilus<fnref n="33" /> [1940] and neuralike<fnref n="34" /> [1941] train a neural network on the spatial dependence of the likelihood in the vicinity of the current live set (only on the live set for neuralike) to construct an estimate of the iso-likelihood contour, and reject points predicted to be outside it before the true likelihood is evaluated, nessai<fnref n="35" /> [1942] takes a similar approach, training a normalizing flow on the set of live points to find a transformation of the iso-likelihood contour into a much simpler hyper-sphere.

### 3.1.4. Practical considerations for choosing an MC sampler

The sampling methods discussed above differ in both their theoretical foundations and practical applications, each offering distinct advantages and limitations in cosmological parameter inference.

Traditional Metropolis-Hastings MCMC uses local random-walk proposals, leading to $O(\sqrt{d})$ scaling in dimensionality $d$. While simple and reliable, these proposals often require many likelihood evaluations, making them computationally expensive, especially for high-dimensional parameter spaces or when the likelihood poorly constrains

<footnote n="19">19 https://github.com/minaskar/zeus</footnote>

<footnote n="20">20 https://github.com/DifferentiableUniverseInitiative/jax_cosmo</footnote>

<footnote n="21">21 https://github.com/ohahn/DISCO-EB</footnote>

<footnote n="22">22 https://github.com/Lbalkenhol/candl</footnote>

<footnote n="23">23 https://github.com/dpiras/cosmopower-jax</footnote>

<footnote n="24">24 https://github.com/JuliaLang/julia</footnote>

<footnote n="25">25 https://github.com/jaimerzp/LimberJack.jl</footnote>

<footnote n="26">26 https://github.com/CosmologicalEmulators/Capse.jl</footnote>

<footnote n="27">27 https://github.com/CosmologicalEmulators/Effort.jl</footnote>

<footnote n="28">28 https://github.com/pyro-ppl/numpyro</footnote>

<footnote n="29">29 https://github.com/farhanferoz/MultiNest</footnote>

<footnote n="30">30 https://github.com/JohannesBuchner/UltraNest</footnote>

<footnote n="31">31 https://github.com/PolyChord/PolyChordLite</footnote>

<footnote n="32">32 https://github.com/joshspeagle/dynesty</footnote>

<footnote n="33">33 https://github.com/johannesulf/nautilus</footnote>

<footnote n="34">34 https://github.com/igomezv/neuralike</footnote>

<footnote n="35">35 https://github.com/mj-will/nessai</footnote>
the parameters. Non-Gaussian posteriors, particularly multi-modal dis-
tributions, can further challenge efficient parameter space exploration.
A good dimensionality scaling can be retained for simple distributions
using decorrelating affine-invariant transforms. For more complicated,
highly non-Gaussian, and multi-modal distributions, ensemble and
sequential MC approaches can help, especially when in combination
with generalized decorrelating transforms such as NFs.

Slice sampling takes a geometric approach by sampling from level sets of the PDF, automatically adjusting to the local structure of the posterior. This self-adaptive behavior results in $O(d)$ scaling without requiring manual tuning, making it effective for problems with widely varying parameter scales, such as fitting galaxy LFs where parameters span multiple orders of magnitude. While computationally inexpensive in low dimensions, its efficiency can degrade in high-dimensional spaces.

HMC improves upon this by using gradient information to make physics-inspired proposals, achieving better $O(d^{1/4})$ scaling. This makes HMC particularly effective for problems with strong parameter correlations, such as hierarchical models in galaxy clustering. However, it requires differentiable likelihoods and careful tuning of integration parameters.

Nested sampling transforms the problem into a one-dimensional integration over nested likelihood contours. Despite its $O(d^3)$ scaling due to constrained sampling requirements when using simple acceptance/rejection sampling, it offers important advantages: direct computation of the evidence $\mathcal{Z}$ (essential for model comparison in DE studies) and natural handling of multi-modal distributions. This makes it particularly valuable in cosmology, GW physics, and gravitational lensing analyses. A better dimensionality scaling can be achieved by using slice-sampling for constrained sampling, or performing acceptance/rejection sampling on a surrogate model.

The choice of sampler often depends on specific requirements like parameter space dimensionality, likelihood characteristics, and the need for evidence calculations [1943–1946]. For problems in which the likelihood is extremely fast, $O(< 10^{-2} \text{ s})$, and the dimensionality is of orders or a few 10 s, simple approaches should suffice, such as affine-invariant MCMC, or a non-boosted nested sampler if a model evidence estimate is needed. For harder problems (very high dimensionality, strong non-Gaussianity, multi-modality), the use of the more advanced methods discussed above, such as preconditioned Sequential Monte Carlo and ML-boosted nested samplers, can prove optimal. For extremely slow likelihoods, $O(> 10 \text{ s})$, the approaches discussed in Section 3.2 such as *surrogate posterior methods* or *simulation-based inference* may be the only viable option, at the cost of some statistical robustness.

An illustration of the problem of dimensionality and the complexity of the distribution can be seen Fig. 51. On it we show how the runtime and the accuracy scale with dimensions for a multi-dimensional Gaussian distribution and the Rosenbrock distribution. For the accuracy test, we use the Wasserstein distance metric [1947], which measures the true distributional similarity between sampler outputs and reference distributions. One can see that all samplers show excellent performance on Gaussian problems across dimensions, but exhibit a dramatic accuracy drop-off for Rosenbrock problems beyond 6–8 dimensions, reflecting the inherent challenge of sampling correctly from complex, highly correlated distributions in higher dimensions.

3.1.5. *Cosmological inference frameworks*

There exist a number of public codes facilitating the integration
of MC samplers such as the ones mentioned above with cosmological
inference pipelines: numerical codes to calculate cosmological observ-
ables (Boltzmann codes such as CAMB [1948], CLASS [1949], and
PyCosmo [1950] and their extensions, as well as emulators) and
experimental likelihoods for cosmological and astrophysical surveys
and data.

The most used ones in the last years are CosmoSIS<fnref n="36" /> [1951] Mon-
tePython<fnref n="37" /> [1916,1917] and Cobaya<fnref n="38" /> [1918] (which uses Get-
Dist<fnref n="39" /> [1952] for analyzing MC results). The choice for a particular
one will be determined by the scientific problem at hand, since each
provides a different combination of features and integrates different
theoretical codes and experimental likelihoods.

Their use is thoroughly recommended since they facilitate not only the interfacing between the different parts of the inference pipeline, but also help in specifying the models and priors as well as the configuration of samplers and theoretical codes, and provide sane defaults for each. They usually either include in their source or automate the installation of cosmological codes and likelihoods, as well as MC samplers.

3.1.6. *Modeling problems: Prior dependence and projection effects*

Analyses of present-day data require a significant modeling effort in order to achieve accurate theoretical predictions that can be compared to data without biasing the results. This includes the modeling of complicated physics and systematic or non-linear effects, which require introducing additional parameters (nuisance parameters) that need to be kept free during the sampling of the corresponding posterior distribution. If the data is not sensitive to all parameters in such a high-dimensional spaces, this can lead to a flat likelihood surface $\mathcal{L}(D|\theta)$ in Eq. (3.1), which in turn can result in a sensitivity to the choice of prior $\pi(\theta)$. If the prior distribution is not informed by past experiments or theory, the dependence of the results on $\pi(\theta)$ may be unwanted and needs to be explored in a sensitivity analysis. A second — more subtle — dependence of the one-dimensional posteriors on the prior is introduced at marginalization. Since the marginalization procedure takes into account the posterior volume, for very flat or very non-Gaussian likelihood surfaces, there may be an (unwanted) up-weighting of regions with large (prior) volume. These so-called projection or prior volume effects can lead to shifts of the marginalized posterior away from the maximum-likelihood estimate or “best fit” point (see also Section 3.6), and may bias the constraints on the physics of the problem if the choice of prior was not sufficiently informed.

An example of this in the context of the $\sigma_8$ tension is the Effective Field Theory of Large Scale Structure (EFTofLSS) (e.g., see Refs. [1177, 1178,1953–1955]), a perturbation-theory based technique to describe the power spectrum of biased tracers up to mildly non-linear scales, that is applied to LSS measurements. In such an approach, several new parameters are introduced, as they enter the perturbative expansion of the galaxy density field. These parameters will be free parameters in the context of an MCMC analysis and, therefore, require a choice of prior. However, the arbitrary nature of the choice of such priors can become a problem for the analysis, since the inferred cosmological parameters can depend on the particular choice of prior. It has been shown how different prior choices, e.g., extending the range of allowed values to extreme intervals, can significantly affect the outcome of the analysis, with projection effects shifting some parameters, such as the primordial amplitude $A_s$ and the amplitude of matter fluctuations, $\sigma_8$, with respect to their true values [1956,1957]. It was further demonstrated that different — albeit theoretically equivalent — choices of the EFT parameterization lead to discrepant inferred cosmological parameters [1957,1958]. Ref. [1959] confirms a significant impact of priors on the nuisance parameters on the inferred cosmological parameter using frequentist profile likelihoods. Approaches to choose well motivated priors on the EFT parameters include the use of a Jeffreys prior [719,1960] or simulation-inferred priors [1961–1963]. Such an effect hinders the robustness of the analysis results, as strong

<footnote n="36">36 https://github.com/joezuntz/cosmosis?tab=readme-ov-file</footnote>

<footnote n="37">³⁷ https://github.com/brinckmann/montepyton_public</footnote>

<footnote n="38">³⁸ https://github.com/Cobaya Sampler/cobaya</footnote>

<footnote n="39">³⁹ https://github.com/cmbant/getdist</footnote>
<figure><img src="image_58.png" /><figcaption>Fig. 51. Illustration of how different sampling methods scale with dimensionality and target distribution complexity. The left panel shows runtime scaling with increasing dimensions, while the middle and right panels display sampler accuracy using two different metrics. Results are shown for both a simple multivariate Gaussian distribution (solid lines) and the challenging banana-shaped Rosenbrock distribution (dashed lines). The traditional accuracy metric combines deviations from mean and function evaluation errors. The Wasserstein accuracy metric uses optimal transport theory to measure distributional similarity between samples and theoretical distributions.</figcaption></figure>

assumptions, only partially motivated by theoretical considerations, need to be done on the choice of priors for the extra parameters.

A similar projection effect, shifting the values of parameters away from their true value, can be introduced also when considering priors on the baryon abundance $\Omega_{b,0} h^2$. Several nuisance effects modeled in LSS observables require the value of this parameter, which is however not strongly constrained by LSS observations. For such a reason the choice of prior for the abundance of baryons can be significant for the analysis, as switching from a uniform prior to a Gaussian constrain as provided by BBN [798] can change the outcome in a significant way [1957].

In the context of the Hubble tension, an impact of prior effects on the inferred value of $H_0$ has been reported and explored for several extended models, for example: EDE (see Section 4.1.1) [1964–1970], NEDE (see Section 4.1.2) [1971,1972], number of relativistic species (see Section 4.1.3) [1973–1975], Brans–Dicke model (see Section 4.3.1) [1976] decaying DM (see Section 4.4.3) [1977] and more.

### 3.1.7. Model comparison criteria

The increasing complexity of cosmological models and datasets necessitates robust statistical frameworks for model comparison and consistency checks. Information criteria and Bayesian methods provide quantitative tools for selecting between competing models while accounting for model complexity and data support.

There exist several information criteria that can be used to compare models [1978]. Common criteria that balance goodness-of-fit of a model against its complexity are the Akaike Information Criterion (AIC) and the Bayesian Information Criterion

$$ \text{AIC} = -2 \ln(\mathcal{L}_{\text{max}}) + 2k, \quad \text{BIC} = -2 \ln(\mathcal{L}_{\text{max}}) + k \ln(N), \quad (3.5) $$

where $\mathcal{L}_{\text{max}}$ is the maximum likelihood, $k$ is the number of free parameters, and $N$ is the effective number of data points on which the likelihood is defined. When comparing a test model versus a baseline one (e.g., $\Lambda$CDM) using ICs (AIC and BIC), we compute the difference in criterion values as

$$ \Delta\text{IC}_{\text{test}} = \text{IC}_{\text{baseline}} - \text{IC}_{\text{test}}. $$

The model with the lowest IC is preferred [1979], meaning, in this sign convention, $\Delta IC > 0$ favors the test model, while $\Delta IC < 0$ favors the baseline model. The value of $|\Delta IC|$ indicates the strength of the preference: $|\Delta IC| \ge 2$ (weak), $|\Delta IC| \ge 6$ (medium), $|\Delta IC| \ge 10$ (strong).

The Deviance Information Criterion (DIC) extends AIC to hierarchical models by accounting for parameter uncertainty:

$$ \text{DIC} = -2 \ln(\mathcal{L}(\hat{\theta})) + 2p_D, \quad \text{with} \quad p_D = 2[\ln(\mathcal{L}(\hat{\theta})) - \langle\ln(\mathcal{L}(\theta))\rangle], \quad (3.6) $$

where $\hat{\theta}$ is the posterior mean of the free parameters, $p_D$ the effective number of them, and $\langle\ln(\mathcal{L}(\theta))\rangle$ the mean log-likelihood over the posterior samples. Unlike AIC and BIC, DIC uses the entire posterior distribution rather than just the maximum likelihood estimate, making it particularly suitable for Bayesian hierarchical models where the effective number of parameters may be less than the actual number due to prior constraints [1980].

The Bayes Factor (BF) provides a fully Bayesian approach to model comparison

$$ \text{BF}_{12} = \frac{\frac{\bar{z}_1}{\bar{z}_2}}{\frac{\int \mathcal{L}_1(\theta_1)\pi_1(\theta_1) d\theta_1}{\int \mathcal{L}_2(\theta_2)\pi_2(\theta_2) d\theta_2}}. \quad (3.7) $$

Jeffreys’ scale is commonly used to interpret $\ln(\text{BF}_{12})$ values: 0–1 indicates inconclusive evidence, 1–3 suggests positive evidence, 3–5 shows strong evidence, and values larger than 5 represent very strong evidence in favor of model 1 [1979,1981]. While AIC and BIC are easily computed from MCMC chains, Bayes factors require accurate evidence estimation, typically obtained through nested sampling or thermodynamic integration methods. As a fully Bayesian approach, the value of the Bayes Factor depends strongly on the choice of prior density, unlike for information criteria, so a prior sensitivity analysis is advisable.

### 3.1.8. Tension metrics

Tension metrics are statistical tools that are used to quantify the agreement (or lack thereof) between the estimation of cosmological parameters obtained with different data sets, beyond the widely applied “rule of thumb”. According to this rule the distance between two posterior 1-d distributions can be quantified in 1-dimensional standard deviations by $N_\sigma = (\mu_A - \mu_B)/\sqrt{\sigma_A^2 + \sigma_B^2}$ where $\mu_{A/B}$ and $\sigma_{A/B}$ refer to the means and variances obtained with each data set. Most studied tension metrics can be grouped into two families: (i) those based on Bayesian evidence [1982,1983] and (ii) those based on the posterior distribution [1984–1986]. While group (i) answers the question about what hypothesis is preferred by the data under the assumed model, group (ii) intends to establish the statistical significance between the posteriors of data sets A and B within the parameter space of both
experiments. For each metric, an estimator is computed and a corresponding probability $P$ is defined that quantifies the agreement or disagreement between data sets in terms of its corresponding equivalent 1-dimensional $N_σ$ defined as: $P = Erf(N_σ/\sqrt{2})$. In most cases the probability $P$ is related to the probability-to-exceed (PTE) of an estimator $Q$ that follows a $χ_d^2$ distribution, defined as

$$ PTE = \int_{Q}^{\infty} \chi_{d}^{2}(x) dx. \qquad (3.8) $$

(i) **Metrics defined on Bayesian evidence ($Z$).** The Bayesian ratio statistic $R$ corresponding to two datasets A and B, and the combination AB of both of them, can be written as [1982]

$$ R = \frac{\hat{Z}_{AB}}{\hat{Z}_{A}\hat{Z}_{B}}. \qquad (3.9) $$

$R \gg 1$ is interpreted as both datasets being consistent, while $R \ll 1$ means that the datasets are inconsistent. However, it has been shown that the $R$ estimator depends on the choice of the prior and can therefore underestimate inconsistencies. This estimator has two primary contributions, one from the unlikeliness of two datasets being in agreement, the information ratio $I$, and another from their disagreement, the suspiciousness $S$ [1983]

$$ \ln I = D_A + D_B - D_{AB}, \qquad \ln S = \ln R - \ln I. \qquad (3.10) $$

Here $D$ refers to the Kullback-Leibler divergence between prior and posterior, quantifying the information gain/compression produced by the given data set. The suspiciousness estimator $S$ remains unaltered under a change of the prior widths as long as this change does not significantly alter the posterior. In addition, if the posterior is a $d$-dimensional Gaussian, the quantity $d - 2 \ln S$ follows a $χ_d^2$ distribution. Under that approximation, the probability of two datasets being discordant by chance can be computed as the PTE of $d - 2 \ln S$.

(ii) **Metrics based on the posterior distribution.** First, we discuss some metrics based on quadratic estimators that require Gaussian posteriors. For these, the probability $P$ of the data sets agreeing is defined as $P = 1 - PTE$ where $PTE$ is the probability-to-exceed defined above. Let us call $\hat{\theta}_A$ the mean of the parameters inferred with dataset A and $\hat{\theta}_B$ the one obtained with dataset B, while $\hat{\Sigma}_A$ and $\hat{\Sigma}_B$ refer to their respective covariance matrices. The method called parameter differences in standard form relies on the quadratic estimator [1984]

$$ Q_{DM} = (\hat{\theta}_B - \hat{\theta}_A)^T (\hat{\Sigma}_B + \hat{\Sigma}_A)^{-1} (\hat{\theta}_B - \hat{\theta}_A), \qquad (3.11) $$

This estimator follows a $χ_v^2$ distribution with $v = Rank[\hat{\Sigma}_B + \hat{\Sigma}_A]$ and can be regarded as a generalization of the rule of thumb. The method called parameter differences in updated form measures how one data set updates the other and is based on [1984]

$$ Q_{UDM} = (\hat{\theta}_{AB} - \hat{\theta}_{A})^{T} (\hat{\Sigma}_{AB} - \hat{\Sigma}_{A})^{-1} (\hat{\theta}_{AB} - \hat{\theta}_{A}), \qquad (3.12) $$

Here $\hat{\theta}_{AB}$ are the inferred parameters obtained with the joint datasets A and B, while $\hat{\Sigma}_{AB}$ is the respective covariance matrix. This estimator follows a $χ_v^2$ distribution with $v = rank[(\hat{\Sigma}_{AB} - \hat{\Sigma}_A)]$ degrees of freedom. Next, we describe the Goodness-of-fit loss which evaluates the likelihood function at the maximum of the posterior of the joint distribution and of each distribution separately [1984]

$$ Q_{DMAP} = 2 \ln \mathcal{L}_A(\theta_{p,A}) + 2 \ln \mathcal{L}_B(\theta_{p,B}) - 2 \ln \mathcal{L}_{AB}(\theta_{p,AB}), \qquad (3.13) $$

Here $\theta_{p,A/B}$ are the Maximum a posteriori (MAP) parameters considering the dataset A/B. The estimator $Q_{DMAP}$ follows a $χ^2$ distribution with $\Delta v = v^A + v^B - v^{AB}$ degrees of freedom, where for a given data set/combination $v = N - \text{tr}[\Sigma_\pi^{-1} \Sigma_p]$, with $N$ being the effective number of data points defining the likelihood, and $\Sigma_\pi$, $\Sigma_p$ the covariance matrices of the prior and the posterior, respectively. This metric is effective at evaluating how well the theoretical prediction fits the data, but does not evaluate directly the disagreement between inferred parameters.

The next method, called exact parameter shift (EPS) does not require Gaussian posteriors and is based on the parameter difference probability [1984,1986]

$$ P(\Delta\theta) = \int P_A(\theta) P_B(\theta - \Delta\theta) d\theta, \qquad (3.14) $$

where $\Delta\theta = \hat{\theta}_A - \hat{\theta}_B$ is the difference between the means of the posterior parameters that correspond to datasets A/B. The probability to identify the tension is calculated by summing over all the values of $P(\Delta\theta)$ over the iso-contour corresponding to no difference $\Delta\theta = 0$

$$ \Delta = \int_{P(\Delta\theta) > P(0)} P(\Delta\theta) d\Delta\theta. \qquad (3.15) $$

Usually, this integral is computed numerically from the posterior samples. In cases of strong disagreement, this estimator may be difficult to compute.<fnref n="42" />

Many authors use the tools presented here to discuss tensions in the inferred value of cosmological parameters with recent data sets [699, 1987,1988]. It has been discussed [1987] that the different metrics answer different questions. For example, $Q_{UDM}$ quantifies how a given data set updates the values of the inferred parameters of another data set, while $Q_{DMAP}$ quantifies the difference of the goodness of fit between the individual and joint data sets to the underlying theoretical model. On the other hand, to quantify the tension between the value of inferred parameters, the appropriate metrics are $Q_{DM}$ and EPS, emphasizing that the former requires posteriors to be sufficiently Gaussian, while the latter does not. Fig. 52 shows two examples of 2-d contours and posteriors that correspond to different data sets and the corresponding $N_σ$ obtained with each metric, together with the ones computed with the 1-d rule of thumb for each parameter. If EPS can be computed, usually the equivalent $N_σ$ that corresponds to DM and EPS indicate similar values, while there is a difference with the ones obtained using the rule of thumb [1987]. Note that the rule of thumb only quantifies the tension in a given parameter, while the metrics presented above quantify the disagreement in the whole parameter space.

## 3.2. Machine learning based inference techniques

**Coordinator:** Jurgen Mifsud

**Contributors:** Alba Domi, Anto Idicherian Lonappan, Benjamin L'Huilier, Celia Escamilla-Rivera, Clecio Roque De bom, Daniela Grandón, Daria Dobrycheva, David Valls-Gabaud, Denitsa Staicova, Elena Giusarma, Filippo Bouchè, Germano Nardini, Ippocratis Saltas, Iryna Vavilova, J Alberto Vazquez, Jacobo Asorey, Jenny G. Sorce, Judit Prat, Konstantinos Dialektopoulos, Leandros Perivolaropoulos, Luca Izzo, Luis A. Escamilla, Matteo Martinelli, Purba Mukherjee, Rahul Shah, and Ruth Lazkoz

In traditional approaches, analytical likelihood functions are employed to characterize the PDFs of observed data, enabling parameter estimation through posterior inference. However, certain datasets exhibit errors that cannot be adequately described by simple analytical distributions, such as multivariate Gaussian distributions. Consequently, the true likelihood in such cases is often complex and lacks an explicit analytical form, making traditional parameter inference methods challenging to implement. To address this issue, various likelihood-free techniques have been developed to circumvent the direct computation of likelihood functions.

<footnote n="41">⁴¹ Reprinted from Physics Letters B, Volume 855, id.138844 Leizerovich, M., Landau, S.J., Scóccola, C.G. Tensions in cosmology: A discussion of statistical tools to determine inconsistencies Copyright (2024), with permission from Elsevier.</footnote>

<footnote n="42">⁴² We point out that in the case of strong disagreement, care must be taken with the sampling of the tails of the distributions.</footnote>
<figure><img src="image_59.png" /><figcaption>**Fig. 52.** Tensions metrics applied to (a) Planck 2018 vs Pantheon Plus + SH0ES (PPS), (b) CC vs Pantheon Plus + SH0ES. In all cases, blue contours correspond to PPS and gray contours to the joint analysis. $N_σ$ refers to the 1-d standard deviations that correspond to a 1-d PDF associated with each metric estimator that quantifies the disagreement between different data sets (see the beginning of this subsection). $Q_{UDM}$ quantifies how the parameters inferred from one dataset are updated when another is incorporated into the analysis. $Q_{DMAP}$ is related to the difference in goodness of fit to the different data of the underlying model. $Q_{DM}$ and EPS quantify the difference in the value of the inferred parameters. $N_{Ω_m}$ and $N_{H_0}$ quantify the tension in each inferred parameter applying the 1-d rule of thumb.⁴¹
Source: Adapted from Ref. [1987].</figcaption></figure>

In standard Bayesian inference, posterior distributions are typically explored using methods such as MCMC sampling, variational inference, or other Bayesian computational techniques. These conventional approaches generally rely on evaluating the likelihood function for the models and parameters under consideration. However, for highly complex and computationally intensive models, the likelihood function may become intractable, and simulations can demand significant time and resources, rendering parameter inference impractical. Therefore, parameter estimation techniques that mitigate or overcome these challenges are highly beneficial, particularly in cosmological studies. Likelihood-free inference has emerged as a promising paradigm for Bayesian inference in scenarios involving complex generative models, leveraging only forward simulations to perform analysis.

With the current generation of large-scale structure observational programs, such as the European Space Agency's Euclid mission [437] and LSST [897], sub-percent-level precision measurements are anticipated. The upcoming decade is thus expected to witness an unprecedented increase in the quantity, diversity, and quality of multi-wavelength astronomical observations of the large-scale structure. This surge in data will necessitate the development of highly sophisticated computational tools. At this stage, the limitations may shift from data quality or availability to the capabilities of statistical and data-driven methodologies. ML techniques have demonstrated substantial potential in addressing the computational challenges associated with traditional statistical methods, positioning them as valuable tools for advancing cosmological analyses.

For instance, current ML techniques were adopted to the scatter in cluster mass estimates [1989] and to distinguish between standard and MG theories from statistically similar WL maps [1990]. Such techniques have also been found to be useful for the next generation CMB experiments [1991], N-body simulations [1992], cosmological parameters inference [1993], DE model comparison [1994], supernova classification [1995] and strong lensing probes [1996].

Cosmology is experiencing an unprecedented growth in the volume and complexity of astronomical datasets, driven by large-scale surveys like Euclid, LSST, and DESI. These datasets are used to construct data vectors for cosmological parameter inference, with which we can shed

light on current parameter tensions within ΛCDM. To fully leverage the wealth of cosmological information encoded in these datasets, ML methods have emerged as a powerful tool.

In this section, we will discuss a number of ML based inference techniques, with an emphasis on those techniques which were adapted to address cosmological tensions. Neural networks are becoming ever more widely employed in physics, including in the field of cosmology. We refer the interested reader to a review of the core concepts surrounding the use of neural networks in Refs. [1997,1998]. In this section we discuss the cosmological applications of artificial neural networks (Section 3.2.1), convolutional neural networks (Section 3.2.2), Bayesian neural networks (Section 3.2.3), and deep learning (Section 3.2.4).

### 3.2.1. Artificial neural networks

The method by which ANN architectures [1661] are used to learn to mimic the iterative process by which an MCMC analysis takes place is described in this section together with some details of the adopted methodology to optimize the ANN structure for comparative analyses. ANNs are non-parametric techniques in the sense that they do not contain the cosmological parameters themselves, unlike the likelihood functions. This may extend cosmological analyses to otherwise overly complex systems, or reduce the computational requirements of regular cosmological model analyses making the approach more accessible. Inspired by biological neural networks, an ANN is a setup composed of a collection of neurons that are organized into layers [1999–2001]. This family of algorithms offers a powerful base on which to apply complex numerical tasks. ANNs are optimized to work in computer systems with a large number of threads which has become a feasible prospect in recent years with the development of very powerful graphical processing units (GPUs), thus making the process of training ANNs competitive with traditional techniques. In recent years, ANNs have become very useful in meeting the accuracy and computational efficiency required for analyses in cosmology. ANNs are constructed with input and output layers to represent the MCMC parameter set that is being sampled and output likelihood respectively. These layers are connected to a series of internal, or hidden, layers that are structured in a way to optimize
<figure><img src="image_60.png" /><figcaption>**Fig. 53.** The network structure used for parameter estimation in Ref. [2002]. The structure of the left panel is for one dataset, while the multi-branch network of the right panel is for multiple datasets {$D_1, D_2, D_3, ..., D_n$} that are from different experiments.</figcaption></figure>

how closely the ANN can imitate the MCMC iterations. Each layer in this network is composed of neurons, each of which is connected to the neurons of the layers preceding and after it, as illustrated in the left panel of Fig. 53 for a single data set.

For the task of inferring cosmological parameter constraints, the observational data is fed to the input layer, then the information of observational data passes through each hidden layer, and finally the cosmological parameters are outputted from the output layer. Specifically, each layer accepts a vector, the elements of which are called neurons, from the former layer as input, then applies a linear transformation and a nonlinear activation on the input, and finally propagates the current result to the next layer.

ANNs have been widely implemented as emulators of cosmological observables, such as the lensing power spectrum [2003,2004], CMB source functions [2005,2006], the CMB temperature anisotropies power spectrum and the matter power spectrum [2004,2007]. The main purpose of these applications is to accelerate the statistical analysis, when expensive Boltzmann codes and big volumes of data prevent us from efficiently sampling the posterior distributions. In particular, CosmicNet I aims to remove bottlenecks from Einstein-Boltzmann solvers (such as CLASS or CAMB) by training neural network emulators that learn the mapping from four ΛCDM parameters to source functions of CMB anisotropies. The resulting trained ANN model is then injected into the CLASS code for public use. This work was later extended to include extensions of a flat ΛCDM, such as DDE, spatial curvature. Emulators are crucial for summary statistics of cosmological probes that lack an analytical prescription that models their dependence on cosmology and other parameters. This is the case of higher-order statistics of cosmic fields, such as the peak counts (a distribution of local maxima in a field) and the one-point PDF. As emulators are only approximations of the target observables, recent works have explored how to propagate the uncertainty derived from such approximation, in order to safeguard parameter inference against biases [2008].

In recent years, likelihood-free analyses relied on ensembles of neural networks as density estimators of the sampling distribution of the data. In Refs. [2002,2009], a cosmological likelihood-free inference technique has been developed. The left panel of Fig. 54 depicts the schematic diagram of the likelihood-free architecture adopted in Ref. [2002] containing the main processes of training and parameter estimation. First, a class object for the cosmological model that contains the simulation method of the measurements is developed, where the simulation method here is used to generate the training set. Then, the initial parameters are set without any constraints in their selection. Subsequently, the training set will be simulated automatically, where two sampling methods can be considered to generate cosmological parameters: sampling uniformly in a hypercube or a hyper-ellipsoid.

It should be noted that the posterior distribution is unknown before the first estimation using the ANN model. Thus, the cosmological parameters here cannot be generated in a hyper-ellipsoid. Therefore, the cosmological parameters are uniformly generated in a hypercube using the set of initial parameters. Consequently, an ANN model will be constructed automatically based on the respective training set. At the same time, the training set is then preprocessed before training the network, while noise will be automatically generated based on the observation errors and added to the training set. The training set will then be normalized, and the training set will be fed to the ANN model, consisting of thousands of epochs.

Cosmological parameters are then estimated using the well-trained ANN model. In simple terms, the ANN model actually learns a mapping between the data space of the measurements and the parameter space of cosmological parameters, as depicted in the right panel of Fig. 54 via a mapping between the data space of measurements and the cosmological parameter space. Therefore, in order to infer the respective posterior distribution, the corresponding distribution of the measurements is fed to the ANN model. Indeed, a large number of data-like samples of the measurements using the observational data are generated and fed to the ANN model to obtain the corresponding ANN chain, similar to the concept of an MCMC chain.

The ANN method in Refs. [2002,2009] was shown to be very successful in estimating cosmological parameters with one data set, which was then extended to multiple data sets by utilizing a multi-branch network. The latter multi-branch network is illustrated in the right panel of Fig. 53, in which it was successfully tested with CMB, SNIa, and BAO datasets in Refs. [2002,2009], as illustrated in the derived results of Fig. 55 with the Planck 2015 data set. Thus, the ANN method was shown to be capable of estimating parameters using data sets of multiple experiments and is a direct competitor of the customary MCMC technique.

### 3.2.2. Convolutional neural networks

With the advent of precision cosmology, we are reaching an era in which it is extremely important to avoid any potential biases in the cosmological analysis when extracting cosmological parameters. For example, most of the time large-scale structure analyses are done using summary statistical observables like 2-point correlation functions, after assuming a cosmology to go from angular coordinates to Cartesian ones. This type of analysis can bias the cosmological parameters and a tomographic analysis can avoid this [2010].

One alternative to the use of summary statistics is the use of CNN architectures. They can be used to make direct inferences on cosmological parameters from astrophysical images. For instance, Ref. [2011] used CNNs to constrain axions as DM from 21 cm line intensity mapping
<figure><img src="image_61.png" /><figcaption>Fig. 54. The left panel depicts the schematic diagram of the likelihood-free architecture adopted in Ref. [2002]. In the right panel, a visualization of a mapping between the data space of measurements and the cosmological parameter space is illustrated.</figcaption></figure>

<figure><img src="image_62.png" /><figcaption>Fig. 55. One-dimensional and two-dimensional marginalized distributions [2002] constrained from Planck 2015 power spectrum. The blue solid lines are the results of the ANN method, the red dashed lines represent those of the MCMC method, and the gray circles are the fiducial values of the cosmological parameters.</figcaption></figure>

probes. In Fig. 56 we show a schematic of the adopted CNN model architecture. It took as input a 64 × 64 × 128 lightcone data cube where each voxel has dimension 4.68 cMpc×4.68 cMpc×1.48 Mhz (cMpc means co-moving Mpc).

If DM were described by axion-like particles (ALPs), the small scale structure would change the value of $\Omega_{m,0}$ leading to a change in $H_0$ and $\sigma_8$, potentially solving the potential tensions. To make a direct inference from CNNs, we need a training set of simulations of the given astrophysical image for different sets of parameters (e.g., axion mass or $H_0$). The input layer is basically the simulation with the angular maps and the frequency maps while the final output layer is the vector with the cosmological parameters inferred (in this case, the axion mass). The input layer is filtered and pooled down before being inputted as a vector to the connected neural network layers. The data is passed through the network a number of epochs updating the corresponding weights each time.

The inferred results, as depicted in Fig. 57, show that the CNN successfully recovered the true values of the axion mass in the testing data with a precision of ~ 20% across a broad range of masses. To evaluate the uncertainty in the recovered values, 100 additional independent simulations were generated for both an axion model and a standard CDM model. For an axion mass of $M_X = 10^{-21}$ eV, a mass uncertainty of $^{+5.94}_{-4.80} \times 10^{-21}$ eV was inferred. It was further concluded that this method could potentially enable the detection of axion DM using SKAO1-Low with 68% confidence if the axion mass is $M_X < 1.86 \times 10^{-20}$ eV. However, these findings depend on the Planck 2015 cosmological parameters and the specific design parameters of future SKAO1-Low configurations.

A CNN was also adopted to perform the typically computationally expensive task of estimating the parameters of GW events. The considered CNN in Ref. [2012] is able to produce posterior distributions that in all cases are compatible with the already published results. The schematic of the respective CNN is depicted in Fig. 58, whereas
<figure><img src="image_63.png" /><figcaption>Fig. 56. The architecture of the CNN adopted in Ref. [2011], with input data cube of 64 × 64 × 128 voxels where each voxel has dimension 4.68 cMpc × 4.68 cMpc × 1.48 Mhz (cMpc means co-moving Mpc).</figcaption></figure>

<figure><img src="image_64.png" /><figcaption>Fig. 57. The distribution of the CNN inferred predictions [2011] on 100 noiseless simulations fixed axion mass of 10⁻²¹ eV, along with the vertical dashed lines showing the 1-σ (68%) credible interval. We should remark that the CNN could predict the axion mass because the CNN is trained only with the data including axions.</figcaption></figure>

the inferred results for the event GW200224_222234 are illustrated in Fig. 59. In the case of event GW200224_222234, all the estimated parameters were found to be in accordance with those published in the GWTC-3 catalog. Having said that, the effective spin was the worst performing parameter. Furthermore, the uncertainty yielded by the CNN on the sky position was found to be too large to accurately pinpoint the event, but it still is unbiased and is expected to produce a first alert for the instruments that look for electromagnetic counterparts to start pointing their instruments to a given patch in the sky. The latter would be more feasible with the upcoming more accurate pipelines which are expected to improve the precision of the sky localization. Overall, the GW event results indicate a good response from CNNs. One of the most important advantages is its computational efficiency, in which the CNN model could provide several posterior samples in

a fraction of the time required by Bayesian inference methods, thus facilitating multi-messenger astronomy.

### 3.2.3. Bayesian neural networks

A BNN that outputs Gaussian PDFs was adopted in Ref. [2013]. The BNN has five input nodes for the five-band grizy photometry. A parameter grid search was performed in order to optimize for free parameters, including the number of epochs, number of layers, number of nodes per layer, learning rate, loss function, activation function, and optimizer. The BNN has a final output node that produces a mean and standard deviation assuming a Gaussian distribution for each photo-z prediction.

A key attribute of the BNN model is the production of photo-z uncertainties, which are needed for using photo-z results in cosmological analyses. It was also found that the BNN produces accurate uncertainties. In Ref. [2013], the HSC Public Data Release 2 [2014], which is designed to reach similar depths as LSST but over a smaller portion of the sky was adopted for the BNN. In total, the data consists of 286,401 galaxies with broadband grizy photometry and known spectroscopic redshifts. The considered galaxy sample covers a redshift of 0.01 < z < 4, however the majority of the sample lies between 0.01 < z < 2.5. In the BNN analysis of Ref. [2013], 80% of the galaxies were used for training, 10% for validation, and 10% for testing. The performance of the BNN technique is clearly illustrated in Fig. 60, in which the BNN was reported to have superior photo-z estimations with respect to other competing models.

The capability of BNNs was confronted in Refs. [2015,2016] on different avenues. The first is how effective BNNs can be in recognizing the distinct features in the power spectrum for a particular modification to the concordance model of cosmology, such as $f(R)$ or Dvali-Gabadadze-Porrati (DGP) [2017]. The second avenue that was exploited in Refs. [2015,2016] is their ability to detect a deviation from ΛCDM in the power spectrum irrespective of the particular modification. Two BNNs with the same architecture were trained in Refs. [2015,2016], the first for five labels divided between ΛCDM and the considered four extensions, whilst the second trained to distinguish between the two labels of ΛCDM and non-ΛCDM. Due to the fact that only four redshift bins were considered in these works, it was beneficial to treat the data as four separate time series and use one-dimensional convolutional layers.

The architecture of the network used to train both the five-label and two-label networks is depicted in the left panel of Fig. 61. Initially, the adopted structure consisted of three 1D convolutional flip out layers with 8, 16 and 32 filters, kernel sizes of 10, 5, and 2 with strides of 2, 2, and 1, respectively. Each of the first two 1D convolutional layers was followed by a max pooling layer with a pool size of 2 and a pooling stride of 2 for the first max pooling layer and a pooling stride of 1 for the second max pooling layer. After both of these max pooling layers there is a batch normalization layer. Following the final convolutional layer there is a global average pooling layer to reduce the filter size to
<figure><img src="image_65.png" /><figcaption>**Fig. 58.** Architecture of the CNN used in Ref. [2012]. The first set of layers are 1-dimensional convolutions and 1-dimensional maxpooling ones, followed by a set of fully connected dense layers. Between each of the dense layers, a dropout of 20% is applied to prevent overfitting during the training stage.</figcaption></figure>

<figure><img src="image_66.png" /><figcaption>**Fig. 59.** Full posterior distribution inferred [2012] from the CNN and MCMC approaches for the GW200224_222234 event.</figcaption></figure>

one in order to pass it to a dense layer with 32 nodes. Finally, after a further batch normalization there is a softmax layer consisting of five or two neurons for either the five- or two-label networks, respectively.

It was found that a five-label network trained to classify between $\Lambda$CDM, $f(R)$ gravity, DGP gravity, wCDM and a “random” class provided more reliable predictions than a two-label network trained to distinguish simply between $\Lambda$CDM and non-$\Lambda$CDM. While generally being less sensitive to variations in the noise distribution, it can also determine whether a power spectrum does not belong to any class included in the training set. Since the selection of the correct model is crucial when performing conventional statistical analyses such as with MCMCs, this ability could prove beneficial in indicating prospective models to consider. However, the network used in this work is currently limited to classification tasks while the notion of model selection

on firm statistical grounds in the context of BNNs remains an open problem.

It could be concluded that BNNs may provide a powerful new means to search for hints of new physics in cosmological datasets. In particular, it is anticipated that they will serve as a powerful “filter”, allowing us to narrow down the theory space before moving on to constrain model parameters with MCMCs while perhaps even signaling the presence of new physics that does not belong to any known model.

### 3.2.4. Deep learning

Neural networks have lately become ubiquitous in cosmology [2001,2009,2019–2032]. In this regard, the LADDER - Learning Algorithm for Deep Distance Estimation and Reconstruction - suite is a novel deep learning algorithm designed to learn the cosmic distance ladder in a model-independent, non-parametric manner [2018,
<figure><img src="image_67.png" /><figcaption>**Fig. 60.** Visualization of the neural network (top left) and BNN (top right) performance compared to the BNNN with outlier removal criteria $\sigma_z < 0.5$ (bottom left) and $\sigma_z < 0.3$ (bottom right) [2013].</figcaption></figure>

<figure><img src="image_68.png" /><figcaption>**Fig. 61.** The right panel illustrates the BNN workflow as adopted in Refs. [2015,2016], whilst the left panel depicts the corresponding BNN architecture employed for both the five-label and two-label classification tasks. The height of each block illustrates the dimension size for each layer, while the number of blocks per layer corresponds to the number of filters. The dense blocks embedded in the first three transparent layers indicate the kernels for the first three one-dimensional convolutional layers scaled by their respective size.</figcaption></figure>

2018]. The schematic overview of the training algorithm is depicted in Fig. 62, whilst the corresponding algorithm is illustrated in Fig. 63. Trained on late-time datasets like the Pantheon SNIa dataset [2033], LADDER incorporates associated errors and complete covariance information. It interpolates from the joint distribution of a randomly chosen subset of the dataset to estimate target variables and errors

simultaneously, effectively handling correlations and the sequential nature of the data. This leads to robust predictions resilient to input noise and outliers, even in data-sparse regions. Optimized by physically motivated metrics such as monotonicity and smoothness, LADDER leverages the Gaussian nature of the observations by employing the Kullback-Leibler (KL) divergence as the loss function during training.
<figure><img src="image_69.png" /><figcaption>Fig. 62. Schematic overview of the deep learning algorithm adopted in Ref. [2018].</figcaption></figure>

<figure><img src="image_70.png" /><figcaption>Fig. 63. The deep learning algorithm adopted in Ref. [2018].</figcaption></figure>

The LSTM network, found to perform best with the Pantheon dataset, is included in the LADDER suite, available publicly on GitHub (https://github.com/rahulshah1397/LADDER) under a MIT License and version 1.0 is archived in Zenodo [2018].

Having learned the cosmic distance ladder, LADDER is a powerful tool for cosmological applications. It can independently verify the consistency of SNIa datasets like Pantheon+ [32] and DES [1257], and serve as a pathology test for qualitatively different datasets, such as 2D vs. 3D BAO. Additionally, it acts as a model-independent calibrator for high-redshift datasets, including GRBs and QSOs. Moreover, previous ML approaches in cosmology struggled with stable error prediction at higher redshifts, rendering them unsuitable for precision cosmological tests. LADDER overcomes these challenges and shows promise in extrapolating beyond available data, useful for simulating intermediate-redshift data or augmenting current data to higher redshifts.

In Ref. [2034] LADDER is employed to recalibrate SDSS BAO and DESI BAO in a model-independent manner, which helps address the longstanding $H_0$ and $S_8$ tensions simultaneously. Traditionally, BAO distances are inferred using the sound horizon at the drag epoch ($r_d$), which is calibrated based on CMB data under the assumption of the $\Lambda$CDM model. In contrast, LADDER, trained on the Pantheon dataset, derives $r_d$ purely from observational data without relying on cosmological model assumptions, instead incorporating merely an astrophysical prior on the absolute magnitude of SNIa ($M_B$). When

combined with Planck 2018 data, this recalibrated BAO framework, along with Pantheon, yields joint constraints on $\Lambda$CDM that significantly alleviate the $H_0$ and $S_8$ tensions, due to an “in-plane shift”, despite their negative correlation. These results highlight the potential biases introduced by CMB-based BAO calibration and demonstrate the robustness of data-driven approaches, offering a promising new avenue for precision cosmology.

Beyond conventional deep learning architectures, more advanced generative models have been developed to tackle complex challenges in cosmological inference. These methods, such as normalizing flows [2035] and diffusion models [2036,2037], extend deep learning capabilities by efficiently learning high-dimensional PDFs and accelerating likelihood-free inference. Unlike traditional MCMC approaches, which often struggle with computational inefficiency, generative models can rapidly approximate complex posteriors while maintaining high accuracy.

Normalizing flows have shown promise in addressing cosmological tensions, such as discrepancies in the Hubble constant and matter fluctuation amplitude ($S_8$). By modeling joint PDFs across datasets, they enhance uncertainty quantification and provide a robust statistical framework for analyzing these tensions [2038,2039]. One example is the EMUFLOW framework [2040], which utilizes normalizing flows to model joint posterior distributions from multiple datasets, enabling efficient constraint combination without increasing parameter space dimensionality. Additionally, methods employing normalizing flows
have been developed for non-Gaussian tension estimation, improving the quantification of agreement levels between different cosmological experiments and offering insights into discrepancies such as the $H_0$ tension [2041]. Further developments in normalizing flow architectures have introduced approaches that incorporate symmetries of the universe, such as translation and rotation equivariance, optimizing parameter inference for cosmological fields [2042]. Other extensions [2043] leverage hierarchical wavelet-based decompositions to improve robustness across spatial scales, while iterative flow-based techniques facilitate the transformation of complex PDFs into tractable forms, improving efficiency in cosmological data analysis.

Diffusion models, leveraging stochastic iterative denoising processes, have also been applied successfully to cosmological parameter inference. These models excel in uncertainty quantification and high-fidelity generative sampling of cosmological fields, particularly in large-scale structure and WL studies. A recent application integrates diffusion models with Hamiltonian Monte Carlo methods [2044,2045], enabling accurate posterior sampling and accelerated convergence in extracting key cosmological parameters such as $\Omega_{m,0}$ and $S_8$. The ability of diffusion models to reconstruct high-fidelity statistical distributions while remaining robust to observational noise makes them valuable tools for precision cosmology. As ML techniques continue to evolve, integrating generative models into cosmological inference is expected to enhance the precision and reliability of parameter estimation, especially for upcoming surveys like LSST and Euclid. These advancements underscore the growing role of deep learning-driven methodologies in shaping the future of cosmological analysis.

### 3.2.5. Conclusion

We reviewed the application of ML methods for the inference of cosmological parameters, emphasizing their versatility and capacity to address the inherent complexity of astronomical datasets. By integrating ML techniques into cosmological inference, we enable efficient and accurate predictions, facilitating faster convergence of traditional algorithms such as MCMC and advancing implicit likelihood or likelihood-free inference using deep learning. These innovations are particularly valuable in handling the high-dimensional, noisy, and often incomplete nature of cosmological data.

Across the subsections of this review, it is evident that the ML techniques considered provide powerful new tools to uncover potential signals of new physics within cosmological datasets. These methods act as a "filter", allowing researchers to effectively narrow down the theory space and refine subsequent analyses with standard techniques like MCMC. Furthermore, ML based approaches may even reveal evidence of new physics that does not fit within the framework of existing theoretical models, opening new avenues for exploration.

While alternative methodologies exist, we emphasize the transformative potential of advanced learning algorithms in optimizing the extraction of information from cosmological data. By demonstrating their utility, we aim to showcase how these novel approaches can complement or enhance traditional inference techniques. The innovative integration of ML into the cosmological framework promises to address current challenges, such as improving parameter constraints and identifying hidden structures or patterns in data, thereby paving the way for more accurate reconstructions of cosmological phenomena.

We also highlight the possibility of yet unexplored ML algorithms that may surpass the performance of those outlined in this review. This underscores the importance of continued exploration and innovation within the ML and cosmology communities. By encouraging the adoption and further development of ML based inference techniques, we aim to inspire new efforts to tackle unresolved questions in cosmology from fresh perspectives. Such endeavors would significantly improve the precision and reliability of cosmological analyses, contributing to a deeper understanding of the Universe's fundamental nature and addressing ongoing tensions in cosmological observations.

In conclusion, the integration of ML methods into cosmological research represents a paradigm shift, offering not only computational efficiency but also novel insights into the underlying physics of the Universe. We urge the community to actively engage in the exploration of these advanced tools, as they hold the potential to revolutionize the study of cosmology in the years to come.

### 3.3. Reconstruction techniques

**Coordinator:** Luis Escamilla, Daniela Grandón

**Contributors:** Adrià Gómez-Valent, Anil Kumar Yadav, Anjan Ananda Sen, Anto Idicherian Lonappan, Ariadna Montiel, Arianna Favale, Benjamin L'Huillier, Biesiada Marek, Celia Escamilla-Rivera, David Benisty, David Valls-Gabaud, Eoin Ó Colgáin, Filippo Bouchè, Iryna Vavilova, Isidro Gómez-Vargas, J. Alberto Vázquez, Jenny G. Sorce, Jenny Wagner, Jurgen Mifsud, Konstantinos Dialektopoulos, Luis E. Padilla, Matteo Martinelli, Miguel A. Sabogal, Purba Mukherjee, Rafael C. Nunes, Rahul Shah, Rocco D'Agostino, Ruth Lazkoz, Víctor H. Cárdenas, and Wojciech Hellwing

Reconstruction techniques have become crucial for understanding the phenomenology of DE, the expansion history, and the growth of large-scale structures in our Universe. They aim to capture specific features or general trends in the data, test a predefined physical model against observations, or explore a flexible, theory-agnostic framework to elucidate the phenomenology of cosmological functions. In general, reconstructions can be categorized into parametric and non-parametric approaches. They mainly differ in their underlying assumptions about a cosmological model and the functional form of the reconstructed quantity. The former imposes an explicit form on the function being reconstructed, thereby reducing uncertainties in the result. However, the a priori assumption of a specific model may limit the flexibility of the reconstruction and introduce biases regarding the properties of the cosmological function of interest. To mitigate this, some studies instead perform a non-parametric ("free-form") reconstruction within a theory-agnostic framework. Non-parametric reconstructions are particularly useful for gaining insights into the dynamics of cosmological functions when the true underlying function is not directly observable but can be inferred from observed data.

Many works on the parametric and non-parametric approaches investigate the impact of $H_0$ priors, to see whether the $H_0$ tension is also imprinted in the derived cosmological reconstructions. Hence, these techniques can shed light on the cosmological tensions and advance our understanding of their effects on a broader context. Moreover, a joint analysis of parametric and non-parametric reconstructions using cosmological observations is valuable to robustly identify departures from the $\Lambda$CDM model that emerge despite the differences of the methods used (e.g., see Ref. [2046-2048] for examples). In summary, these methods serve as data-driven frameworks to decode the phenomenology of cosmological functions from complex datasets, and offer an alternative approach to distinguish between competing cosmological models.

In what follows we present some examples of reconstructions and how they can be classified according to their methodology.

#### 3.3.1. Parametric reconstructions

In a parametric reconstruction, the desired physical quantity is modeled using an analytical function with a set of free parameters that describe its behavior across cosmic time or spatial scales. The free parameters of the postulated function, and the cosmological parameters, are hence jointly inferred in light of the observational data. Hence, parameterizations are a crucial tool for modeling unknown physical phenomena in a mathematically tractable way. As an example, in the absence of a fundamental and well-defined theory of DE, several functions, i.e., $w(z)$, $\rho(z)$, have been parameterized in a broad number of different ways. These parametrizations track the dynamics of a
DE back in time simply extending the cosmological model with extra parameters. While parametric methods are computationally efficient and easy to interpret, they are inherently limited by the assumptions built into the chosen parameterization. A poorly chosen model may lead to biased results that fail to capture the true complexity of the underlying physics. As such, model selection criteria for various parametric candidates are critical in assessing the robustness of the results.

### 3.3.2. Phenomenological parameterizations

Phenomenological parameterizations focus on capturing specific behaviors, often without direct derivation from fundamental physics. Their primary goal is to test whether certain trends (such as evolving equations of state, deviations from General Relativity, or modified growth of structure) are compatible with observational data. These parameterizations provide a practical way to explore a broad range of possibilities without requiring a detailed theoretical foundation. Some well-known examples include:

* CPL parameterization [297,298] for the DE equation of state

$$w(z) = w_0 + w_a z/(1 + z), \qquad (3.16)$$

which is the general go-to parameterization used to test if any dynamic behavior can be present in DE (this parameterization is further explained in Section 2.3.4). The DESI collaboration recently employed this parameterization in their analysis of year-one BAO data [699] and found evidence suggesting a preference for a DDE model with a significance exceeding $2\sigma$.

* The Jassal–Bagla–Padmanabhan (JBP) parameterization [2049], also for the DE equation of state

$$w(z) = w_0 + w_a z/(1 + z)^2, \qquad (3.17)$$

which serves as a modification of the CPL designed to exhibit different behaviors at large redshifts.

* Interaction term between DE and DM [2050–2053]

$$Q = \xi H \rho_{DM}, \qquad (3.18)$$

with the parameter $\xi$ governing the strength of the interaction. It is proposed as a way to represent the interaction as a function of time ($H$) and energy density ($\rho_{DM}$). Given its phenomenological nature, this approach can be extended further, as demonstrated in Ref. [2054], where the interaction term was modified to depend not only on DM but also on DE. In this case, the interaction takes the form $Q = \xi H(\rho_{DM} - \alpha\rho_{DE})$ with the new parameter $\alpha$ regulating the strength of the dependence on DE. Moreover, the authors in Ref. [2052] discuss the interaction between DM and DE in the context of the second law of thermodynamics and its consequences on DE evolution. For a more complete discussion on DE-DM interactions please refer to Section 4.4.3.

* A parametric form for the DE density $X(z)$ was developed by Refs. [2055–2057] and further tested in Refs. [2046,2058]. The quadratic parameterization corresponds to

$$X(z) \equiv \frac{\rho_{DE}(z)}{\rho_{DE}(0)} = 1 + (4x_1 - x_2 - 3) \left(\frac{z}{z_m}\right) - 2(2x_1 - x_2 - 1) \left(\frac{z}{z_m}\right)^2, \qquad (3.19)$$

where $z_m$ is the maximum redshift in the data set, $x_1$ and $x_2$ are the free parameters of the model to be constrained in light of the data. When $x_1 = x_2 = 1$, the parameterization reduces to $\Lambda$CDM where the DE density is constant $X(z) = 1$. This approach aims to study evidence for DE evolution by allowing its energy density to change with redshift. Reconstructing $X(z)$, instead of $w(z)$, is advantageous since this quantity is more directly traced by the luminosity distance. A cubic parameterization was also explored in Ref. [2058] to allow for an extra degree of freedom and to

investigate whether further transitions with respect to $X(z) = 1$ are found in the reconstruction.

* The growth rate parameterization [2059]

$$f(z) = \Omega_{m,0}(z)^\gamma, \qquad (3.20)$$

where $\gamma$ is named the growth index, which for $\Lambda$CDM has a value of $\sim 0.55$. A recent application of this parameterization can be found in Ref. [2060], where an analysis using $f\sigma_8$ and Planck CMB revealed a tension of $4.2\sigma$ with the concordance model's expected value of $\gamma = 0.55$, suggesting a potential need for modifications to the way the standard model predicts structure growth.

* The gravitational slip [2061]

$$\eta = \frac{\Phi}{\dot{\psi}}, \qquad (3.21)$$

where $\Phi$ and $\dot{\psi}$ are the metric potentials in the perturbed Einstein equations. If $\eta$ deviates from 1 there may be an indication for MG. In Ref. [2062] four different models for $\eta$ were analyzed using forecasts from WL, galaxy clustering, and supernova data to evaluate how effectively future surveys could constrain the slip parameter. Furthermore, in Ref. [2063] a deviation from GR was found using GW data along with the slip parameter.

As we can see, phenomenological parameterizations allow for a broad exploration of deviations from the standard model without committing to a specific underlying theory.

### 3.3.3. Physically motivated parameterizations

Physically motivated parameterizations, in contrast, are derived from fundamental theories and ensure consistency with known physics. These models are often constructed from Lagrangian formulations, field equations, or extensions of Einstein's General Relativity, ensuring that theoretical principles like energy conditions and stability constraints are satisfied. Some of these parameterizations which can be found in the literature are:

* Scalar fields for DE. A well-known one being the Quintessence [295], where a scalar field $\phi$ with a potential $V(\phi)$ leads to an equation of state

$$w(z) = \frac{\dot{\phi}/2 - V}{\dot{\phi}/2 + V}, \qquad (3.22)$$

having the peculiarity of behaving as, in which the kinetic term has the wrong sign, $\dot{\phi}^2/2 \to -\dot{\phi}^2/2$, leading to $w(z) > -1$. Complementary to it we have the Phantom scalar field [2064] where $w(z) < -1$. The combination of them, resulting in the interacting-two-scalar-field model named Quintom [2065–2067] has also been studied due to its peculiarity of being able to cross the Phantom-divide line ($w = -1$) [2068,2069]. K-essence scalar field [2070], which can be seen as an extension of Quintessence where non-canonical kinetic terms are included in the Lagrangian.

* Horndeski theories of MG [665,2071–2073], which provide the most general scalar-tensor theories leading to second-order field equations.

* $f(R)$ theories [2074] which modify the Ricci scalar in the Einstein–Hilbert action

$$S = \int d^4x \sqrt{-g} \left( \frac{1}{16\pi G} f(R) + \mathcal{L}_m \right), \qquad (3.23)$$

where the $f(R)$ can take many forms. In Ref. [2075] two different parameterizations for $f(R)$ were studied, them being $f(R) = R - 2\Lambda(1 - e^{-R/\Lambda b})$ and $f(R) = R - 2\Lambda(1 - \frac{1}{1+(R/\Lambda b)^n})$, finding concordance with the standard model when using BAO, SN and CC data. Similarly, in Ref. [2076] a preference for the standard model is recovered when the two models $f(R) = R - \beta/R^n$ and $f(R) = R - \alpha \ln R - \beta$ are used in tandem with simulated GW data.
* For DM models, we encounter WDM [2077–2081], where a nonzero velocity dispersion modifies the matter power spectrum. Another example is Self-Interacting DM [2082–2085], in which a cross-section $\sigma/m$ accounts for self-interactions, affecting structure formation at small scales while preserving the CDM behavior on larger scales.

Both phenomenological and physically motivated parameterizations aim to describe the evolution of cosmic components, test deviations from $\Lambda$CDM, and assess the viability of new models against observations. Phenomenological parameterizations provide a flexible way to analyze trends, while the physically motivated ones ensure consistency with fundamental physics. The choice between the two depends on whether the focus is on empirical data fitting or theoretical consistency.

### 3.3.4. Model-independent parameterizations

While both phenomenological and physically motivated parameterizations help constrain the behavior of cosmological functions, they inherently require assumptions about the functional form. However, these assumptions can limit the flexibility of the analysis, as the chosen parameterization may not fully capture the true underlying features of the function being reconstructed. To overcome this limitation, one can employ more sophisticated parameterization methods that minimize the number of underlying assumptions. Some of these approaches, referred to as model-independent methods, remain parametric but offer significantly greater flexibility, allowing them to capture a wider range of possible features in the data.

#### 3.3.4.1. Basis representation.
Basis-representation parameterizations provide a flexible and systematic way to reconstruct cosmological functions by expanding them in terms of a chosen set of basis functions. Common examples found in the literature include:

* Given its simplicity and ease of use, the Taylor expansion has been applied to many different cosmological functions. As an example, in Ref. [2086] the DE equation of state, through its pressure, was expressed as a truncated Taylor series

$$p = p_0 + \sum_{n=1}^{N-1} \frac{1}{n!} \frac{d^n p}{d\rho^n} \bigg|_0 (\rho - \rho_0) + \mathcal{O}[(\rho - \rho_0)^N]. \quad (3.24)$$

This idea can be extended even further. By performing a Taylor expansion of the Hubble parameter around $z=0$, one arrives at the cosmographic approach [2087–2090], which will be explored in more detail in a later subsection.

* In a manner similar to the Taylor series expansion, one can instead adopt a different basis, such as a Fourier series. In Ref. [2091], this approach was applied to the DE equation of state, revealing a preference for an oscillatory behavior at late-times. The equation of state was found to cross the phantom divide multiple times, which has been found to be a stable solution for the Kinetic Gravity Braiding model proposed in Ref. [2092], a minimally coupled Horndeski model. This however does not hold for other single scalar field models, where gravitational instabilities appear when crossing the phantom divide.

* The Padé approximation could be seen as an extension to the Taylor series, being itself a rational of two polynomial series. It has the advantage of a faster convergence rate, but potentially using more new terms/parameters. Its uses are well varied in the literature, going from the luminosity distance [434,2093–2096], the Hubble expansion rate [2096], the DE EoS [2095] and even $f(R)$ [2097]. In Ref. [2093] the luminosity distance was expressed as $d_L = \frac{a_1 z}{1+b_1 z+b_2 z^2}$, which is usually called the (1, 2) Padé approximant for $d_L$, finding some agreement with Planck's CMB data at small redshifts but allowing DE to be dynamical.

* Wavelet expansion. Wavelets work as a basis and have the peculiarity of being localized in frequency and configuration space,

which is useful for capturing both global trends and local fluctuations. In Ref. [2098], a wavelet-basis was used for the DE EoS in the next manner

$$w(z_j) + 1 = \sum P_l \psi_l(z_j), \quad (3.25)$$

with $z_j$ being the redshift points at which $w$ is calculated, $P_l$ the coefficients of the expansion and $\psi_l$ the wavelet function. In this particular work the wavelets used were the Haar ones (where $\psi_l(x) = 1$ for $0 \le x \le 1/2$, $\psi_l(x) = -1$ for $1/2 < x \le 1$ and $\psi_l(x) = 0$ otherwise) and, using SNIa, WMAP and BAO data, a hint for dynamical DE was found.

* Reconstructions using an orthonormal basis set of functions (ONB) has been used to reconstruct the normalized cosmic expansion function, $H(z)/H_0$, from the Pantheon dataset [2099], without assuming a specific cosmology. In this approach, the luminosity distance function $D_L(a, c)$ was expanded into orthonormal basis functions $\phi_a(a)$ as

$$D_L(a, c) = \sum_{a=0}^{N_B-1} c_a \phi_a(a) = c \circ \Phi, \quad (3.26)$$

where $a$ is the scale factor. The $c_a$ elements of $c \in \mathbb{R}^{N_B}$ denote the weights of the $N_B$ basis functions. These terms are summed up in the short-hand notation of the right hand side. The choice of the ONB can be adapted to the problem at hand, as has been done in Ref. [2099], but other general choices like Chebycheff polynomials can be used. Then the values of $c_a$ are obtained by minimizing the $\chi^2$ function [2100]. Using a special ONB adapted to approximate a $\Lambda$CDM universe with its first few basis functions, [2099] showed that the reconstructed cosmic expansion function from the Pantheon set of SNIa still yields very broad confidence bounds, such that $\Lambda$CDM and other competitive models of DE are still statistically consistent with the dataset.

#### 3.3.4.2. Interpolation methods.
Another model-independent approach to parameterization, similar to the basis representation method, involves interpolation techniques. In this approach, a set of “positions” is defined and then connected in a specific manner, with the chosen connection method determining the type of interpolation used. These positions act as free parameters, whose optimal values are determined through Bayesian parameter inference algorithms, such as MCMC or Nested Sampling.

An example of such a class of methods consists in binning the function one wants to reconstruct and treat each one of its binned values as a free parameter of the analysis, e.g., one could aim at constraining the value of $w(z)$ at a given number of redshift $z_l$. This completely removes assumptions made on the trend of the function, but significantly increases the number of free parameters to take into account. When using this approach the bins are usually connected with hyperbolic tangents to preserve continuity. As a function it can be represented as

$$w(z) = w_1 + \sum_{l=1}^{N-1} \frac{w_{l+1} - w_l}{2} \left( 1 + \tanh \left( \frac{z - z_l}{\xi} \right) \right), \quad (3.27)$$

where $N$ is the number of bins, $w_l$ the amplitude of the bin, $z_l$ the position where the bin begins in the $z$ axis and $\xi$ is a hyperparameter which governs how “smooth” is the transition from one bin to the other (a higher value represents a higher smoothness). Being one of the most utilized model-independent interpolation methods it has a wide range of applications such as in reconstructing: the neutrino mass [2101]; the primordial power spectrum and inflaton potential [2102]; the DE energy density [2103,2104]; the DE EoS [2103,2105,2106]; the interaction kernel in an IDE model [2021,2107]; a DE perfect-fluid model [2108]. Another example of a binned analysis (or interpolation) is the reconstruction of the functions $\mu(z)$ and $\Sigma(z)$, which describe deviations from General Relativity (e.g., see Refs. [2109–2111]). In
Refs. [2112,2113], joint constraints on these functions and on the DE density parameter $\Omega_{DE}(z)$ were obtained combining CMB, BAO, and SN data. The results of this analysis show hints for deviations from General Relativity, mainly in the $\Sigma(z)$ function which encodes deviations from the gravitational lensing effect expected in the standard model, while also reducing both $H_0$ tension to $\sim 2\sigma$ (when comparing the inferred value of $H_0 = 69.44 \pm 1.3 \text{ km s}^{-1}\text{Mpc}^{-1}$ with the SH0ES one) and the $S_8$ tension (by inferring a value of $S_8 = 0.780 \pm 0.033$ which is in line with DES Y3's measured value of $S_8 = 0.769 \pm 0.016$).

Another choice of interpolation is commonly referred to as “nodal reconstruction” [2114]. Here the interpolations are done using linear, cubic, or higher order splines, to fill in the gaps between a certain number of nodes. The simplest example is the linear interpolation, that is, given two coordinates $(z_i, f_i)$ and $(z_{i+1}, f_{i+1})$, the function behaves as follows

$$f(z) = f_i + \frac{f_{i+1} - f_i}{z_{i+1} - z_i}(z - z_i), \quad z \in [z_i, z_{i+1}]. \qquad (3.28)$$

This method has been used to reconstruct several important cosmological quantities, i.e., the primordial power spectrum [2115–2118], the DE energy density and the DE EoS [2103,2119], the expansion rate of the Universe [2120] and the neutrino mass [2101]; just to mention a few.

An unusual type of interpolation makes use of GP. The mathematical formalism behind it was already presented in Section 3.1, and further details of its applications in reconstruction will be discussed in the following subsection. This interpolation works in a similar manner to the nodal reconstruction, but with the obvious distinction that the connection between the varying nodes is made through GP. This approach can mitigate some issues that the other two approaches present, such as the choice of the bins' width and positions [2121], while also presenting the advantage of being infinitely differentiable. Given its relatively new implementation it has not been used widely. In Ref. [2121] was first implemented with the equation of state parameter, and in Ref. [1994] was used to reconstruct the interaction kernel of an IDE model.

### 3.3.5. Cosmography

The degeneracy among various theories proposed in recent years to explain DE has driven the investigation of methods that enable the study of cosmic expansion without relying on predefined cosmological models. This is the case of the well-known cosmographic approach [2087–2090], which depends only on the cosmological principle and uses series expansions of the luminosity distance around the current time. The power of the cosmographic method arises from the fact that it involves observables that can be directly compared with data and guarantees independence from any assumed DE equation of state. Thus, the cosmographic method has been extensively utilized to distinguish between various theoretical models that appear similar when compared to observations [2122–2125].

However, the cosmographic method faces two main challenges that may limit its effectiveness as an accurate tool for describing cosmic expansion. The first challenge is the necessity for a substantial amount of data to distinguish between the cosmological constant and a DDE. The second issue concerns the use of high-redshift data to explore potential deviations from the standard cosmological model. The latter aspect conflicts with the core principle of the standard cosmographic technique, which relies on a Taylor expansion series centered around $z = 0$. As a result, large error propagations due to convergence issues often compromise the effectiveness of the method. Over time, various alternatives to the standard cosmographic method have been explored to address these limitations. One approach involves using auxiliary variables and developing expansion series for cosmological observables based on re-parametrizations of the redshift variable converging as $z \gg 1$ [434,1394,2126].

An alternative way to tackle the aforementioned challenges is to utilize rational polynomials that can heal convergence problems typical of the standard Taylor-based cosmography. A prominent example of this technique is offered by Padè polynomials, which can be calibrated to maximize the convergence radius, thereby enabling a more stable fitting process at high redshifts [2093–2095]. The benefits of Padè polynomials have recently been exploited in various theoretical contexts to explore potential deviations from Einstein's gravity [2096,2127,2128]. Furthermore, the cosmographic method based on Padè polynomials, as well as the auxiliary y-redshift parameterization, has been adopted to evaluate in a model-independent way the cosmological tensions in the $H_0$ and $\sigma_8$ measurements suggested by recent observations [2129]. Another strategy to extend the convergence radius of the cosmographic series and address the subjectivity in truncating the expansion, which may still be an issue with the Padé method, involves using Chebyshev polynomials. It has been shown that these polynomials can significantly reduce the uncertainties in higher-order terms of the cosmographic series, thus offering a precise description of the Universe's evolution at late-times [2130,2131].

The Weighted Function Regression Method, first employed in Ref. [422], addresses the subjectivity involved in choosing the truncation order of cosmographical expressions. It does so by considering all cosmographical orders when reconstructing cosmological functions, and it applies robust Bayesian weights to penalize the use of additional parameters. In practice, the contribution of higher orders can be neglected since their weights are significantly suppressed. This method has been used to reconstruct the Hubble and deceleration functions, as well as to estimate the Hubble and deceleration parameters in a model-independent way with low-redshift data, see Ref. [422,2132], respectively.

### 3.3.6. Non-parametric methods

Non-parametric reconstructions are data-driven methods that do not rely on a specific functional form of the reconstructed function and hence strive to minimize model assumptions. In other words, the aim is to describe the data without imposing strong a priori constraints. Among the most common applications of non-parametric methods we find GP regression, ANNs, the local regression smoothing method, and the iterative smoothing method. These methods vary in their flexibility to reconstruct highly non-linear relations between variables and how the uncertainties on the reconstruction are obtained. We proceed to describe them and their applications in cosmology.

#### 3.3.6.1. Gaussian processes.
Both powerful and versatile tools in ML and statistical analysis [2133,2134], GP are a Bayesian technique that generalizes distributions over functions, extending Gaussian distributions into function space [2135]. This approach allows for the reconstruction of a function $f(x)$ at every point $x$ using an observational data set $\{(x_i, f(x_i) + \sigma_i) | i = 1, \dots, N\}$, without needing to assume a predetermined specific functional parameterization. The reconstructed function and its derivatives are Gaussian random variables with a mean $\mu(x)$ and variance cov $[f(x), f(x)]$ at each data point $x$. The functions at different points $x$ and $\tilde{x}$ are related by a covariance function $K(x, \tilde{x})$, commonly referred to as a kernel, which depends on a small set of hyperparameters. Although there is a wide range of covariance functions available in the literature [2134,2136], the hyperparameters are generally constant, as their values characterize the function's smoothness rather than modeling such behavior. Most cosmological applications of GP use a zero mean function and optimize the hyperparameters [2137]. A strict Bayesian approach is to marginalize over the hyperparameters. As for the mean function, while stationary processes can be described by a constant mean function, non-stationary ones such as the cosmic distances or growth may not. The authors in Ref. [2138] showed how marginalizing over a reasonable family of mean functions and over the hyperparameters yields a robust and unbiased result. An application of this full marginalization approach
is to use GP as a forward model. One can generate *untrained* GP by
sampling the hyperparameters, and calculating the likelihood. This al-
lows forward modeling of quantities such as the growth rate $f\sigma_8(z)$ and
distances while combining several datasets [2138–2144]. In cosmology,
GP techniques have been successfully employed to reconstruct the
dynamics of cosmological functions with minimal physical assumptions
about the Universe’s geometry or the nature of its main components.
These include the background the expansion rate of the Universe
$H(z)$ [419,422,423,429,2047,2145–2149] as can be seen in Fig. 64,
the deceleration parameter $q(z)$ [1376,2150–2152], the distance duality
relation [2153–2155], and the cosmological jerk parameter [2156]. To
explore the dynamics of DE, examples include reconstructions of the
DE energy density [2046,2058,2142], the scalar field dynamics [2157],
the equation of state $w(z)$ [425,2158,2159], and interaction between
DE and DM [2160–2163]. These applications have served to study
MG models [2147,2147,2164–2176], cosmic curvature [2177–2184],
and GWs [672,2185–2187], inflaton speed of sound’s profile [2188],
and many other perspectives [424,430,2101,2189–2198]. In particular,
the authors in Ref. [2188] adopt a slightly different approach. They
demonstrated how mildly-informative physical priors can be imposed
on a GP reconstruction in a Bayesian way so that robust constraints
on physical parameters can be extracted along with a non-parametric
reconstruction. In summary, all these applications have proven par-
ticularly valuable for testing the robustness of ΛCDM and identifying
potential deviations from our current understanding of the Universe.

3.3.6.2. *Artificial neural networks.* ANNs have also seen a steep increase in their application to cosmology (for further information on this please refer to Section 3.2.1, Section 3.2.2, Section 3.2.3 and Section 3.2.4). The two main advantages of ANNs over other non-parametric techniques are that they can capture (highly) non-linear relations in complex datasets, and in principle do not require the data to follow a specific statistical distribution. For this reason, the cosmological community has developed various neural network reconstructions of cosmological functions. Some cosmological functions reconstructed with ANNs include the distance modulus $\mu(z)$ from the Pantheon and Pantheon+ compilation [1994,2018,2020,2023,2025,2034] and the Hubble function $H(z)$ [2001,2020,2026,2199,2200]. In particular, a novel approach was developed in Ref. [1994] where the authors train a recurrent neural network to learn the mapping from redshift $z$ to distance modulus $\mu(z)$ using the Pantheon SNIa sample. This model, combined with a BNN, is also able to propagate the uncertainties into the predicted function $\hat{\mu}(z)$ and thereby into cosmological constraints. This approach was also implemented to calibrate GRBs at high redshift [2021]. Other applications include ANNs reconstructions using LSST simulated data [2201], CNNs to reconstruct the BAO signal [2202], ANNs for growth data $f\sigma_8$ [2001,2019], and velocities of rotation curves [2203]. For visual reference, in Fig. 64 an ANN was used to reconstruct $H(z)$. Other approaches generate ANN models for small datasets and their errors using feedforward neural networks with Monte Carlo Dropout and hyperparameter grid optimization [2001] or GAs [2023].

As with other non-parametric reconstruction methods, the goal is
to generate an ANN model based on observational data and then
compare the resulting reconstruction with theoretical predictions from
parametric methods to allow for robust conclusions. In this direction,
different cosmological models have been analyzed using neural net-
work reconstructions. Examples include ΛCDM [2001,2020,2184], the
CPL model [1994,2001,2020,2201], Chaplygin gas models [1994], and
Horndeski gravity [2026], among others.

3.3.6.3. *Iterative smoothing method.* The authors in Ref. [2204] used this non-parametric method to smooth supernova data using a Gaussian smoothing function. This method aims to reconstruct cosmological quantities such as the expansion rate, $H(z)$, and the equation of state of DE, $w(z)$, in a model-independent manner. The only assumptions

made by this method are the smoothing scale and a guess background
model for the quantity under study. However, it is necessary to use
a bootstrapping method to determine the optimal guess model. In
Ref. [2204], the authors used an iterative method to estimate the guess
model. They started with a simple cosmological model, such as ΛCDM,
as the initial guess model and then the results were used as the next step
in the iteration. With each iteration, it was expected the guess model to
become more accurate, thus giving a result that is less and less biased
toward the initial guess model used. Indeed, it was noted that using
different models for the initial guess does not affect the final result as
long as the process is iterated several times.

The method requires careful consideration of the smoothing scale. A very small smoothing scale gives an accurate but noisy guess model, therefore after a few iterations, the result will become too noisy to be of any use. It is better to use a larger smoothing scale for smoother results. On the other hand, the bias of the final result decreases with each iteration, since with each iteration one gets closer to the true model. The bias decreases non-linearly with the number of iterations. In Ref. [2204], points out that after about 10 iterations, for moderate values of the smoothing scale, the bias is acceptably small. Also, beyond this, the bias still decreases with the number of iterations but the decrease is negligible while the process takes more time and results in larger errors on the parameters. In Refs. [2205–2207] the method was improved to reconstruct the distance modulus in a model-independent way and then employed in Refs. [2208,2209], to name just a few examples. In particular, this method was applied to reconstruct the cosmic expansion history and growth to test curvature, DE, and GR [2207,2210–2212].

3.3.6.4. *Local regression smoothing and simulation extrapolation (Loess+Simex).* The LOcally wEighted Scatterplot Smoothing (Loess), also known as local polynomial regression, is a non-parametric method for data analysis that does not require specifying a predefined relationship between dependent and independent variables. It generalizes standard least-squares methods and is widely used for non-parametric simple regression in various disciplines. Loess aims to depict the global trend of a dataset by fitting low-degree polynomials to subsets of data around each observation, giving more weight to points closer to the target observation using a Kernel-based weighting system. Typically, first or second-order polynomials are used since higher orders do not significantly enhance results and increase computational complexity. The process is applied iteratively to cover the entire data range, resulting in a comprehensive trend depiction. Important aspects of the method include selecting the number of data points for each fit through a smoothing parameter called span, $s$, determining the polynomial degree, and choosing the weight function form. Additionally, confidence intervals are constructed by calculating the variance of the fitted values, assuming Gaussian distributed errors. Despite some bias in the estimation, the cross-validation technique provides accurate confidence intervals around the Loess curve. To account for the observational errors on real data, additional methods are needed. To this end, the authors in Ref. [2213] proposed to combine the Loess method with the simulation-extrapolation method, also called Simex, which is a simulation-based method designed to minimize bias resulting from the inclusion of covariates that are prone to errors. Estimates are derived by introducing additional measurement errors, using a form of resampling. This resampling helps identify the pattern of measurement error. After estimating the pattern, final estimates are obtained by extrapolating back to the scenario where there is no measurement error.

Loess and Simex were first used together in cosmology in Ref. [2213] to reconstruct the expansion history $H(z)$. Later, this technique was also used to reconstruct galaxy rotation curves [2214], the distance modulus [2215], and other cosmological quantities [2048,2216,2217], with also accurate results.
<figure><img src="image_71.png" /><figcaption>**Fig. 64.** Non-parametric reconstructions of the Hubble function $H(z)$ using CCs. Two different approaches were used: GP regression (left) and ANNs (right). The ANN reconstruction was made following the approach in Ref. [2001]. The analysis models observational uncertainties using Monte Carlo Dropout, as also employed in Refs. [1994,2184,2201]. The cosmic chronometer data from [410] and theoretical predictions from the ΛCDM model are included for reference.</figcaption></figure>

## 3.4. Bio-inspired algorithms in model selection

**Coordinator:** Reginald Christian Bernardo

**Contributors:** Anto Idicherian Lonappan, Antonio da Silva, Arrianne Crystal Velasc, Celia Escamilla-Rivera, David Valls-Gabaud, Dinko Miklovic, Erika Antonette Enriquez, Filippo Bouche, Isidro Gómez-Vargas, J. Alberto Vázquez, Jenny G. Sorce, John K. Webb, Jurgen Mifsud, and Renier Mendoza

We briefly discuss GAs (Sections 3.4.1–3.4.2), their applications to cosmology (Section 3.4.3) and potential relevance to the understanding of cosmological tensions (Section 3.4.4).

### 3.4.1. Genetic algorithm

A GA is a biology-inspired optimization strategy that mainly takes elements of natural evolution in order to single out one solution that is the fittest from a pool of similarly naturally selected individual solutions. GAs are powerful optimization methods, called meta-heuristic because they do not use derivatives to find the optimum, and they guarantee to find the best solution under certain conditions, despite the challenges posed by local optimality [2218]. This has been used in a wide variety of scientific problems such as in high energy physics [2219] and GW astronomy [2220], and is known to be able to find global optimum and resolve tiny differences between seemingly degenerate solutions to a problem. GA is particularly well-suited to bypass issues related to complex, high-dimensional parameter spaces and multimodal functions. In cosmology, it was introduced as a means to overcome the bias in selecting a cosmological model in order to infer the properties of DE [2221]. This was then followed by Ref. [2222] and Ref. [2223], which has further marketed GA as an alternative tool for cosmological analysis through the estimation of uncertainty. An excellent recent introduction to GA for cosmological parameter estimation is given in Ref. [2224].

GAs operate with a population of individuals (possible solutions), where each individual is characterized by a chromosome, which in turn is described by a set of genes. These individuals make up the population that evolves over successive generations. The key ingredients of GA are described in Table 4, while the evolutionary process is illustrated in Fig. 65.

In essence, GAs sample from a population of individuals ranked by a fitness scale. The fittest individuals are more likely to survive and reproduce through crossover and mutation, creating offspring with improved traits. The best or elite individuals are stored to guarantee that the population is improving through the generations. This iterative process continues until the optimal solution, or the fittest individual,

<figure><img src="image_72.png" /><figcaption>**Fig. 65.** Flowchart illustrating the key steps in GAs [2225].</figcaption></figure>

emerges. In cosmology, the gene can be thought of as each cosmological parameter, such as $w_0$ and $w_a$ in CPL cosmology, and the chromosome as the string of genes, e.g., $(H_0, \Omega_{m,0}, w_0, w_a)$ [2224]. A GA can, of course, be applied in a variety of ways; notably, in the pioneering works [2221–2223] the GA has been used in the context of grammatical evolution (GE), where the elements and key ingredients of GAs are utilized to search for solutions within function spaces, alleviating the arbitrariness of cosmological parameterization. More recently, GAs have been used to eliminate the arbitrary choice of a kernel function in GP regression [2226] and to optimize neural networks [2023]. It has also found applications in spectroscopic modeling [2227–2229], and in conjunction with information criteria for model selection [2230].

### 3.4.2. Machine learning and GA variants

In this section, we dig a little further into the details on GAs in order to get a better grip of its intricacies, and use it as a template to introduce some of its more familiar variants that have also been considered in the field. The basic steps of GAs are summarized in the pseudocode shown in Algorithm 1, based on the notation in Ref. [2223] and Ref. [2225].

Algorithm 1 highlights that GAs are fundamentally an optimization strategy, inspired by the principles of natural selection and genetics. In cosmology, however, particularly through the lens of grammatical evolution, GAs transcend traditional optimization by serving as an ML method. The resulting combined method utilizes the adaptive search capabilities of GA and GE to incrementally evolve a grammar that approximates cosmological functions, such as cosmic expansion and growth rates, from input data. By discovering functional forms that capture the underlying dynamics of cosmological phenomena,
<table><caption>**Table 4**
Key ingredients of GAs [2225].</caption><tr><td>Fitness-function</td><td>Determines the fitness of individuals in a population. This function can be tailored to the specific problem, such as minimizing the Euclidean distance or optimizing the likelihood in cosmological applications. The fitness function is a predefined metric used to rank solutions.</td></tr><tr><td>Selection</td><td>Defines the portion of the population that will advance to the next generation. A common method is the ' roulette wheel' selection, where fitter individuals have higher chances of being selected for reproduction.</td></tr><tr><td>Elitism</td><td>Ensures that the individuals that passed the selection process can be directly added to the next generation. Elitism guarantees that the best fitness value will not get worse after every generation.</td></tr><tr><td>Crossover</td><td>Refers to the process of combining the genetic information of two parents to produce offspring. This mechanism allows for the exchange of genes between individuals, promoting the inheritance of favorable traits while introducing genetic diversity into the population.</td></tr><tr><td>Mutation</td><td>Involves the random alteration of an individual's genes. This step is crucial for introducing new genetic variations, ensuring that the population does not stagnate and continues to evolve toward better solutions. This allows GAs to escape local solutions and obtain global solutions.</td></tr></table>

**Algorithm 1 GA pseudocode.**

**Notation:** $f$: fitness function, $[b_1, b_2]$: population range, $D$: population dimension, $r$: mutation rate, $N$: population size,

$G$: max. number of generations, $m$: individual index, $t$: generation index, $M(t)$: population at generation $t$, $m^*$: best solution

**Input:** $f$, $[b_1, b_2]$, $D$, $r$, $N$, $G$

**Output:** $m^*$, $f(m^*)$

1. Define the maximum number $G$ of generations and the mutation rate $r$

2. Set the generation counter $t \leftarrow 0$

3. Generate an initial random population $M(0)$ consisting of $N$ individuals having dimension $D$ and range $[b_1, b_2]$

4. **while** $t < G$ **do**

5. Compute and save the fitness $f(m)$ of each individual $m$ in the current population $M(t)$

6. Sort $M(t)$ with respect to the fitness of each $m$

7. Define selection probabilities $p(m)$ for each $m$ in $M(t)$

8. Choose a proportion of $M(t)$ based on the selection probabilities $p(m)$

9. Apply crossover on parent individuals to produce offspring

10. Apply mutation on offspring with a probability based on the mutation rate $r$

11. Store the new generation $M(t + 1)$ of individuals

12. Select the elite individuals in $M(t + 1)$ to be preserved for the next generations

13. $t \leftarrow t + 1$

14. **end while**

15. Return the best individual $m^*$ of the final generation and its fitness $f(m^*)$

GAs in this context demonstrate their dual role in both optimizing and learning, revealing intricate patterns that may be governing our Universe. This versatility is notably owed to the pioneering works of Refs. [2221–2223], which fully fleshed out GA’s flexibility for extracting insights on the elusive DE given late-time data. This opens up GAs to the same wealth of applications ML has been applied to in cosmology, such as in the calibration of very high redshift observables [2021], non-parametric cosmological reconstructions [2215], Bayesian deep learning for DE [1994,2231,2232], and in testing the validity of routines such as the cosmographic approach [2233].

Other than ML, another broad scientific discipline where GAs are prominent is optimization, specifically under the categories of metaheuristic, nature-inspired, and evolutionary algorithms. Within this realm, GAs have various variants, such as memetic algorithms, which are designed to prevent premature convergence by incorporating local search strategies [2234]. Given their inherently distributed nature, GAs stand to benefit significantly from the advent of quantum computers and their parallelism. Quantum GAs exploit quantum mechanical phenomena such as superposition and entanglement to implement

quantum evolutionary concepts and operators [2235]: quantum chromosomes, entangled crossovers, quantum elitism, etc. The quantum nature of these algorithms introduces a non-zero probability of adding new genetic material to the population, thereby mitigating premature convergence [2236].

In recent years, several quantum GAs have been developed [2236–2239], featuring either a fully quantum architecture or the integration of quantum operators within a classical GA framework. Implementing these innovative GAs on actual quantum devices has yielded enhanced performances compared to their classical counterparts, despite the limitations of today’s noisy quantum hardware [2236,2238]. The first application of quantum GAs in cosmological analysis is currently underway, aiming at minimizing the chi-squared function of different cosmological probes: SNIa, baryon acoustic oscillations, CMB [2240]. This development marks a significant step forward, combining the strengths of quantum computing and GAs to tackle complex problems in cosmology.

One of the more familiar variants of GAs in the astrophysics and cosmology community is particle swarm optimization (PSO). PSO is inspired by the social behavior of birds flocking or fish schooling, where each particle (or candidate solution) adjusts its position in the search space based on its own experience and the experience of neighboring particles. The algorithm involves iteratively updating the velocities and positions of the particles, guiding them toward the best solutions found so far. This method has been used in Ref. [2241] to infer cosmological parameters given CMB data and in Ref. [2242] to calibrate semi-analytic galaxy formation models in a fixed cosmological background. PSO has been employed even earlier in astrophysics, such as in Ref. [2243] to search for periodic orbits in three-dimensional galactic potentials. Recently, the method has found application in GW analysis as seen in Refs. [2244].

It is now becoming clearer that GA and PSO approaches are only the tip of the iceberg, or perhaps an opening to Pandora’s box, of mathematical optimization techniques available for astrophysical and cosmological research. The recent review by Refs. [2244] discusses several optimization methods, including particle swarm optimization, in high energy and astrophysics applications. In the following sections, we shall introduce a Bayesian tool inspired by PSO, the Approximate Bayesian Computation-Sequential Monte Carlo (ABC-SMC) [2245–2247], which has been applied recently in cosmological model selection and parameter estimation [2248,2249]. Work on a comparative study of nature-inspired optimization methods applied to cosmological parameter estimation, including MCMC, GA, Improved Multi-Operator Differential Evolution [2250], and the Philippine Eagle Optimization Algorithm [2251], is in progress, showing optimistic results for parameter estimation and cosmological reconstruction [2252].

### 3.4.3. GA for cosmology

In cosmology, GAs have mostly been used in the context of grammatical evolution [2221–2223]. This was motivated by the fact that the Universe’s most dominant component at late-times is dark, and a non-parametric reconstruction via GA-GE seems particularly well
suited for such a problem. In grammatical evolution, one deals with a set of grammar functions, often polynomials are considered as well as trigonometric, logarithm, and exponential functions, that are exposed to the GA optimization scheme. The output of GA-GE is thus a functional approximation of the “true” function that has spawned the observed data; this can then be characterized with a fundamental model, such as ACDM or $w_0w_a$ CDM CPL, if the data were cosmological in nature. Uncertainty estimation in the context of GA-GE was also studied in Refs. [2221–2223], suggesting various converging estimates by a Fisher matrix formalism, bootstrapping a la Monte Carlo, and a path integral approach.

Recently, GA-GE and other ML strategies were used in Ref. [2253] to draw insights from late-time expansion data. As with any ML approach, however, the output of the reconstruction via GA-GE can generally be quite dependent on the choice of the hyperparameters, for GA-GE, this would be the grammar functions. This can be put to good use if there is prior knowledge of the functional behavior of the underlying fundamental model, such as in Ref. [2254] where trigonometric grammar functions have been utilized to search for oscillating features in the primordial power spectrum. This turned out to be a promising avenue to address simultaneously the tensions in the Hubble constant and the matter power spectrum [2255,2256]. Furthermore, the GA-GE approach has been applied to SN, baryon acoustic oscillations, CC, RSD, and normalized Hubble rate measurements to symbolically reconstruct the cosmological functions in Ref. [2257], and was also used to forecast constraints on the ACDM model in Refs. [2258–2260].

Another way to use GAs in cosmology that has gained some traction recently is to directly use it to estimate the cosmological parameters, e.g., phenomenological DE parameters $w_0$ and $w_a$, in a given model [2224]. In this way, the uncertainty can similarly be estimated using a Fisher matrix approach or via bootstrapping, resulting in an approximation to the posterior of the cosmological parameters, that can be used with MCMC to strengthen the robustness of the analysis. We illustrate this very briefly with a quick GA computation (Fig. 66), to show the constraints obtained by GAs and MCMC in spatially-curved ACDM and flat CPL given cosmic chronometer measurements [2261] and SN observations [32,33,2262]. For the GA parameters, we consider the likelihood $L$ as a fitness function, a “tournament” type selection with a rate of 30%, an adaptive mutation rate of (80%, 20%), and a “scattered” crossover type with a probability of 50%. In our notation for adaptive mutation, the first number in the tuple (a, b) denotes the fraction a of the genes in a chromosome that are going to be mutated for low quality solutions (ranked accordingly by fitness), and the second number corresponds to the fraction b of chromosomes that are to be mutated for high quality ones.

In non-flat ACDM, it can be seen that the GA and MCMC results agree with each other on the inferred cosmological parameter space to a high degree, including correlations [2263]. The GA output population can furthermore be used to chip in to the results; as shown in this case it is also consistent with MCMC. However, in the more complex CPL model, the analogous constraints reveal a slight hint of disassociation between the results of MCMC and GA, particularly with the DE parameters. The disagreement is not to be worried about regardless since the deviations are well within the inferred confidence regions of either method. This quick exercise nonetheless shows GA in action in cosmological parameter estimation, and reinforces one of our messages that GAs can act as a supporting tool to MCMC for cosmological analysis. This is especially true for models that are more complex than ACDM, despite the different nature of MCMC methods that sample the posterior distribution, while GAs maximize the likelihood.

GAs and their variants have been considered in applications that assess improve the output of other methods, such as in optimizing Neural Networks [2023,2029] and speeding up Bayesian inference [1941]. Similarly ABC-SMC and GA have been applied to tackle the issue of

kernel selection of GP [2264] when it is applied for cosmological reconstruction [2226,2265]; leading toward an ML hybrid that no longer depends on any one arbitrary choice of a kernel function. Furthermore, ABC-SMC has shown to be quite useful for parametric reconstruction and model selection, as shown in Ref. [2248], in the end catering to its users with one model and its parameters that best suit the data. Another notable application of GAs is given in Ref. [2266] to test the viability of an emergent universe scenario using Planck data and BBN.

### 3.4.4. GA for cosmological tensions

We emphasize that GAs for cosmological parameter estimation are not meant as a replacement for MCMC. GAs are rather a support tool for the traditional MCMC, offering unique advantages given its ability to navigate complex, non-linear, and high-dimensional parameter space, for example, to have faster parameter estimations and more information to define the priors of an MCMC. Subsequently, GAs can become very useful in tackling cosmological tensions and systematics, especially given the fact that biases might not be traceable by other methods. This has been teased out in recent results to be discussed.

Ref. [2254] used GAs to search for local features in the primordial power spectrum using Planck data. GAs were applied as a reconstruction tool by using its main ingredients in a grammar space of functions, such that the overall implementation looks for the most suitable parameterization of the power spectrum all while constraining the parameters of each parameterization. It resulted in significant improvements to the fit, compared with prior Planck likelihoods [2254]. Such local features were shown to be able to reduce simultaneously both the Hubble and the $\sigma_8$ tensions, giving a possible unified fundamental solution to the cosmic tensions traceable to inflation [2255,2256]. Notably, GA's capacity to break degeneracy between otherwise acceptable solutions has proven to be an invaluable lens for identifying local features in the primordial power spectrum.

On the other hand, Ref. [2267] used GAs to constrain the comoving sound horizon at the baryon drag epoch. They also checked it against traditional recombination codes or numerical fitting via the Eisenstein-Hu approach. GAs were again applied in an independent cosmological model fashion by working on grammar functions. Subsequently, GA constraints to the sound horizon at baryon drag prevented usual cosmological parameter biases that affect traditional methods. Because the Hubble constant from BAO measurements relies on constraints to the sound horizon at baryon drag epoch, GAs may well be able to help in understanding the cosmic tensions and the role of systematics. The observable is also directly influenced by the matter density, giving this a route for GA to chip into the $\sigma_8$ tension.

Ref. [2268] furthermore used GAs to test cosmology at the perturbative level. They used GAs to constrain growth rate observations, based on a synthetic ACDM cosmological model. GAs led easily to the reconstruction of the underlying cosmology and ruled out other often considered phenomenological and MG models including wCDM, designer f(R), and the Hu-Sawicki model. This impressive application highlights another aspect of GAs — the ability to see a signal through noise, even when the noise is as huge as in growth rate observations. [2269] used GAs together with cosmological background and perturbations data, leading to a support for phantom DE behavior in the dark sector in order to alleviate the cosmic tensions.

Ref. [2249] used a GA variant ABC to pit together Hubble constant priors that are representative of the Hubble tension. This has seen the solution consistently evolve toward the direction of the Planck values, independent of the shape of the priors and the data sets (CCs, SN, and RSD) used in combination. However, this is heavily influenced by the use of CCs which were considered as a baseline data set in the analysis, since CCs prefer low values of the Hubble constant. Nonetheless, this serves as another example of when GAs can be used to choose which priors are more consistent for given data sets. An earlier work has used ABC in the same vein [2248], but in a more extensive way by using the GA variant to select between ACDM and phenomenological
<figure><img src="image_73.png" /><figcaption>**Fig. 66.** Constraints on a spatially-curved ΛCDM and a flat CPL ($w_0w_a$CDM) models with CCs [2261] and SN [32,33,2262]; blue corresponds to MCMC results; black to GA-Fisher (GA Fisher matrix uncertainty estimation hybrid); orange are the localized samples (trimmed outliers) in the GA final population containing the best solution. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

parameterizations of DE, sometimes referred to as XCDM. The results of this work have been quite intriguing, all evolved toward XCDM with a preference for Hubble constant values consistent with the Planck constraint, regardless of the choice of priors, including Hubble constant priors, and data sets considered. However, this has similarly relied on CCs as an anchor of the data, and the resulting low values of the Hubble constant may be traced to this.

Fig. 67 presents a comparison of $H_0$ and $S_8 = 0.8 \sqrt{\Omega_{m,0}/0.3}$ [6] constraints from MCMC, GA, and ABC-SMC using a combination of late-time data. The results highlight the ability of GAs, discussed extensively in this review, to address cosmological tensions by corroborating MCMC outcomes and even providing its own distinctive results. All cases examined show remarkable consistency between GAs and MCMC, with GAs successfully capturing correlations within the parameter space. Notably, all independent methods considered generally support a low $S_8$ and a high $H_0$, although this conclusion is contingent on the datasets used. Additionally, ABC, which compares ΛCDM and wCDM while constraining their parameters, demonstrated narrower posterior distributions compared to MCMC, particularly when the algorithm extends over numerous generations before selecting a preferred model [2248]. Nonetheless, ABC constraints remain consistent due to the method's capacity to mitigate prior dependence [2248,2249].

To sum up, GAs offer several significant advantages as they are particularly effective for solving complex optimization problems where traditional methods falter. GAs are robust and flexible, capable of finding optimal or near-optimal solutions in large, multi-dimensional search spaces. GAs are also inherently parallel, allowing to explore multiple solutions simultaneously and increase the likelihood of discovering high quality solutions efficiently. Understandably, GAs also have some disadvantages, such as computational efficiency, requiring significant processing power and time, especially for large-scale problems. Its performance also relies heavily on the choice of its hyperparameters such as population size, mutation rate, and crossover rate. Regardless, this is certainly one method that the community may find worth further investing into, for its flexibility and applicability to an incredible range

of scenarios, embodied by the various ways GAs and their variants have been applied in astrophysics and cosmology.

When looking for a needle in a haystack, one may turn to MCMC to identify the most probable region where the needle could be. Reversely, one may turn to GAs to identify this needle, systematically sifting through the haystack and evolving toward the solution with precision and persistence. This needle may well be the solution to the cosmic tensions, if not a huge hint to this cosmological conundrum.

### 3.5. Inference from cosmological simulations

**Coordinator:** Lei Zu

**Contributors:** Alan Heavens, Andrew Liddle, Anto Idicherian, Lonappan, Antonio da Silva, Benjamin L'Huillier, Chi Zhang, Daniela Grandón, David Benisty, Elena Giusarma, Filippo Bouchè, Houzun Chen, Jenny G. Sorce, Jenny Wagner, Jurgen Mifsud, Luz Ángela García, Marika Asgari, Nikolaos E. Mavromatos, Oleksii Sokoliuk, Ruth Lazkoz, Vasiliki A. Mitsou, and Wojciech Hellwing

In this section, we explore cosmological simulation and their connection to cosmological tensions. We begin by introducing the role of cosmological simulations in modeling nonlinear scale physics and the relevant tools. Next, we examine how cosmological simulations serve as a powerful tool for investigating physics beyond the ΛCDM. Finally, we discuss key observations and associated tensions at small scales, showing how cosmological simulations bridge theoretical predictions and observational data.

#### 3.5.1. Cosmological simulation

When studying the large-scale structure evolution of the Universe, the density perturbation $\delta = \frac{\delta\rho_m}{\rho_m}$ becomes much greater than unity ($\delta \gg 1$) at small scales like $k \gtrsim 1 \hbar/\text{Mpc}$. In this nonlinear regime, linear perturbation theory fails to accurately describe the evolution of structures. To overcome this limitation, N-body simulations, a particle based method, are employed to model the behavior of phase space particles—representing a group of DM particles within a given cosmological
<figure><img src="image_74.png" /><figcaption>**Fig. 67.** $H_0$ and $S_8 = \sigma_8\sqrt{\Omega_{m0}/0.3}$ constraints on a spatially-flat ΛCDM and u:CDM models using MCMC, GA, and ABC. Data sets used are CCs [2261], standardized distances (Pantheon+ /SH0ES [32,33,2262]), growth rate measurements (compiled in Ref. [2270]), and baryon acoustic oscillations (DESI year-1 [699]). Red vertical bands are from Planck [192], blue from KiDS-450 [697,850,858,903], and gray from Ref. [34]. An asterisk in the superscript in the label denotes that the Pantheon+ /SH0ES data was not considered. In the ABC label, the bold faced model is the preferred one. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

framework. Additionally, some cosmological simulations incorporate
hydrodynamics to account for the effects of baryons. These simulations
produce detailed predictions of structure formation, enabling direct
comparisons with observational data. While computationally intensive,
cosmological simulations are essential for understanding the Universe
at scales dominated by nonlinear effects. They serve as a critical link be-
tween theoretical cosmological models and observations. By providing
insights into nonlinear physics, cosmological simulations play a pivotal
role in addressing tensions between the standard ΛCDM model and
observational data.

**3.5.2. Nonlinear effects**

N-body simulations describe structure formation at scales domi-
nated by the nonlinear effects of gravity. In these simulations, “parti-
cles” represent DM, characterized by properties such as mass, position
and velocity. The initial conditions for these simulations are derived
from theoretical models of the early Universe, typically set at high
redshifts (z ~ 100) when the Universe was nearly homogeneous.
On small scales like k ≥ 10 h/Mpc, hydrodynamic effects become
significant. The inclusion of hydrodynamics is necessary in simulations
to accurately capture the influence of baryonic processes. GADGET
(GAlaxies with DM and Gas intEracT) is one of the most widely used N-
body and hydrodynamics simulation codes [2271–2273]. It combines
N-body methods for gravitational interactions with Smoothed Particle
Hydrodynamics (SPH) for modeling gas dynamics, enabling the study
of both collisionless systems (e.g., DM) and hydrodynamic phenomena.
Building on GADGET, other advanced codes have been developed, such
as AREPO [2274], which uses a moving mesh approach instead of
SPH, and GIZMO [2275], which replaces modern SPH with a mesh-
free hydrodynamics solver for improved accuracy. SWIFT is also a
fully open-source highly-parallel, versatile, and modular coupled hy-
drodynamics, gravity, cosmology, and galaxy-formation code [2276].
In addition, Adaptive Mesh Refinement based codes such as RAM-
SES (Refined Adaptive Mesh with Static Expansion) [2277] are also
extensively used for cosmological and astrophysical simulations. Due

<figure><img src="image_75.png" /><figcaption>**Fig. 68.** The matter power spectrum with the same ΛCDM cosmological parameters for $k \sim 0.1–10$ h/Mpc at $z=0$. The green solid line shows the linear results. The black solid line represents the *N*-body simulation results by using the publicly available code **GIZMO** while the orange dotted (blue dash-dotted) line is calculated by the nonlinear analytical code **Halofit** (HMCode). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

to the computational intensity of high-resolution simulations, several
groups have made their results publicly available. Examples include the
**Millennium Simulation** [2278], **Horizon Run Simulations** [2279],
**Illustris and IllustrisTNG** [2280,2281], **Bolshoi Simulation** [2282],
**EAGLE** [2283,2284], **Quijote Simulations** [2285]. Researchers can
utilize these simulated results for their own science targets.

Using some N-body simulation results as benchmarks [1021,2286],
several semi-analytic models are developed to predict the nonlinear
matter power spectrum, such as Halofit [889,890] and HMCode [891,
2287]. These semi-analytic models significantly reduce computational
<figure><img src="image_76.png" /><figcaption>**Fig. 69.** The ratio of the nonlinear matter power spectrum at $z=0$ in DM-proton scattering scenario [2290]. The red solid lines denote the emulator developed in Ref. [2290], the black dashed lines are the N-body simulation results from GIZMO, the blue dashed-dotted lines are the correction from HMCode and the orange dotted lines are the correction from HaloFit.</figcaption></figure>

time while remaining applicable across a wide range of cosmological parameters. Halofit employs simulation results to construct fitting formulas, while HMCode, a refined version of the halo model, introduces modifications that enhance its consistency with simulation outcomes. Additionally, ML techniques have been increasingly applied as emulator to predict nonlinear structure formation. For example, the **Quijote Simulations** have been widely employed in ML applications for cosmology, providing large-scale training datasets for deep learning models aimed at improving parameter inference and nonlinear structure formation predictions [2288,2289].

**Fig. 68** shows the matter power spectrum in the standard ΛCDM model for $k \sim 0.1–10$ h/Mpc at $z = 0$. The linear matter power spectrum (green solid line) significantly underestimates the power at nonlinear scales $k > 0.1$ h/Mpc. Nonlinear results from N-body simulations (black solid line, GIZMO) align closely with analytical models such as Halofit (orange dotted line) and HMCode (blue dash-dotted line). Within the ΛCDM framework, both Halofit and HMCode achieve accuracy within 10% compared to the N-body simulation. HMCode also exhibits superior performance when incorporating extensions beyond ΛCDM, such as baryonic feedback, WDM, massive neutrinos, DDE, and MG [891, 2287]. However, if the initial matter power spectrum or the late-time dynamical evolution significantly deviates from the ΛCDM scenario, the reliability of HMCode predictions is doubtful. As shown in Fig. 69, when addressing new physics beyond ΛCDM, at nonlinear scales, both HMCode and Halofit fail to provide accurate nonlinear corrections, see Ref. [2290] for details. In such cases, conducting cosmological simulations correctly becomes essential to study the structure evolution in the specific model.

### 3.5.3. Cosmological simulations beyond ΛCDM

To explore new physics beyond the ΛCDM framework at nonlinear scales, various groups have developed cosmological simulations in different cosmological scenarios. One of the most studied extensions is the WDM model. The free-streaming motion of WDM, due to its thermal velocity, suppresses structure formation on small scales. Numerous cosmological simulations of WDM [1722,2291–2293] have investigated the resulting mass distributions and halo properties. These studies reveal suppressed mass distributions below a scale determined by the

mass of WDM particle, and shallower density slopes systematically compared to CDM.

Another critical area of interest is the total mass of neutrinos, which influences the evolution of the total energy density and introduces free-streaming effects. Several studies have incorporated massive neutrinos into cosmological simulation codes using additional particle sets or relativistic approaches [2294–2296]. A higher total neutrino mass results in a stronger suppression of the matter power spectrum at small scales. For instance, the suppression of the total matter power spectrum at $k \sim 1$ h/Mpc is approximately 5% for $m_v = 0.06$ eV, and increases to 25% for $m_v = 0.3$ eV, see Ref. [2295] for details.

In addition to neutrino effects, cosmological simulations have also been extended to study alternative gravity models and DE effects beyond CDM. Various computational tools have been developed to simulate structure formation under MG theories, including ECOSMOG [2297], MG-GADGET [2298], ME-GADGET [2299], and MG-AREPO [2300], MG-PCOLA [2301], MG-QUIJOTE [2302]. ECOSMOG incorporates MG models, including $f(R)$ gravity, chameleon models, and symmetron theories, using an AMR framework to achieve high-resolution nonlinear simulations. ME-GADGET allows for an arbitrary phenomenological (with varying $H(a)$ and $G_N(a)$) and IDM (with varying $m_{DE}(a)$) models to be used in the gravity-only simulation setting. MG-PCOLA is a fast N-body solver extending the COLA method to MG, balancing speed and accuracy for large-scale structure studies. Recent MG-PCOLA simulations incorporating massive neutrinos reveal that their inclusion leads to additional suppression of structure growth. In particular, in scalar-tensor models, this suppression enhances the damping of small-scale power beyond what is expected from standard neutrino free-streaming, reinforcing the interplay between neutrinos and MG. MG-QUIJOTE is a large suite of MG simulations tailored for ML applications and parameter inference, covering various MG scenarios and offering a benchmark dataset for cosmological studies.

Furthermore, if DM has additional non-gravitational interactions — for instance, with baryons/electrons [2290,2303,2304], neutrinos [2305–2307], or dark radiation [2308–2310] — this will introduce additional pressure that suppresses small-scale structure formation. These interactions modify the Boltzmann equation in the early Universe by adding energy and momentum exchange terms, creating a pressure that counteracts gravitational collapse. This suppression occurs primarily during the early Universe, when matter densities are high enough to trigger collisions. As the Universe expands, the interaction rate falls below the Hubble rate, i.e., DM decouples. After that, the interactions become negligible on galactic scales. Thus the primary deviations between IDM and CDM arise during the early Universe, influencing the initial conditions for structure formation. A common approach involves calculating the linear evolution of the matter power spectrum first, including DM interactions, during the early phase when density perturbations remain small and linear theory is applicable. These results are then used as initial conditions (at $z \sim 100$) for later cosmological simulations when DM interactions can be neglected. Notably, interactions often imprint unique oscillation features in the initial power spectrum, similar to those observed in baryons. However, these oscillations diminish during nonlinear evolution as the Universe evolves. As a result, these distinctive oscillation structures are expected to be observable primarily at high redshifts [2290,2310,2311].

### 3.5.4. Cosmological simulations in nonlinear observation and cosmic tension

Cosmological simulations model particle behavior and predict DM properties at nonlinear scales, providing a crucial tool for comparing theoretical predictions with observations such as WL, Lyman-$a$, halo density profile, upcoming cosmic 21 cm, and so on. Various discrepancies between ΛCDM model and observations have emerged at different scales, including the $S_8$ tension, core-cusp problem, too-big-to-fail problem, the diversity of dwarf galaxy rotation curves, missing satellites problem, and so on. Cosmological simulations serve as a
powerful bridge between theory and observations, helping to address discrepancies and exploring the physics beyond standard ΛCDM. In this section, we discuss key nonlinear observations and the related tensions between ΛCDM and the measurements.

### 3.5.4.1. Weak lensing.
WL provides a direct method for mapping the late-time large-scale structure of the Universe by statistically analyzing the shape distortions of numerous background galaxies induced by foreground matter fields. By calculating the two-point correlation functions of the shapes of galaxy pairs, WL measurements have widely been used for exploring the matter distributions in our local Universe. WL is sensitive to the scale down to $k \gtrsim 1$ h/Mpc at low redshift (typically $z < 2$), where nonlinear gravitational evolution significantly affects the matter power spectrum [855,856,860,862,869]. On smaller scales ($k \gtrsim 10$ h/Mpc), baryonic processes, such as gas cooling, feedback from star formation, and AGN, leading significantly systematic uncertainties to WL measurement. To mitigate uncertainties in baryonic physics, many analyses mask or marginalize over these small-scale regions [862,869, 989]. Several other approaches have been developed to address this uncertainty. For example, one method that preserves small-scale measurements is the Bayesian Model Averaging (BMA) framework presented in Ref. [849]. In this approach, BMA combines the individual posterior distributions of multiple competing baryonic feedback models derived from hydrodynamical simulations to provide robust cosmological constraints. Furthermore, the nonlinear evolution of cosmic structures at late-times introduces strong non-Gaussianities in the matter density field, requiring estimators beyond two-point functions to capture this additional information [829]. These estimators, often referred to as higher-order or non-Gaussian statistics, are useful for breaking parameter degeneracies [2312,2313], self-calibrating systematics [983], and overall improving cosmological constraints compared to two-point functions [837,842,844,847,2314,2315]. The theoretical modeling of these estimators relies heavily on cosmological simulations, both to characterize their dependence on cosmological parameters and to build the likelihoods, where the covariance matrix is typically estimated from hundreds of simulations at a fixed cosmology [883]. An alternative to adding estimators of higher-order statistics is to use the full information from the field using Bayesian hierarchical modeling, e.g., see Ref. [2316–2320], or Bayesian simulation-based inference based on neural compression of the field, e.g., see Ref. [2321,2322], which avoid the need to estimate covariance matrices. Both these approaches rely heavily on simulations. We refer the reader to Section 2.1 for a detailed discussion on WL.

In ΛCDM framework, WL surveys have reported an $S_8$ tension (see Section 2.1 for details): the $S_8$ parameter measured from WL survey is lower than the value inferred from Planck CMB data, with a significance of $2\cdot3\sigma$ [859,869,2323]. Despite some uncertainties in observational measurements, various beyond ΛCDM models have been proposed to address this issue, such as suppressed small-scale matter power due to interactions between DM and dark radiation [2310,2324–2326], DM-baryon interactions [2327], or DM-neutrino scattering [2328,2329]. Cosmological simulations play a crucial role in addressing these issues. Upcoming more precise large scale surveys like Vera C. Rubin Observatory [897], Euclid [437] and China Space Station Telescope (CSST) [2330] are expected to provide crucial insights and help resolve these tensions.

### 3.5.4.2. Lyman-α forest.
Lyman-α forest refers to a series of absorption lines observed in the spectra of distant QSOs. These lines are caused by the interaction of light from the QSO with intervening clouds of neutral hydrogen gas along the line of sight. As a result, the Lyman-α forest is widely used as a tracer of neutral hydrogen gas distribution and, by extension, the underlying matter distribution, assuming a correlation between mass and neutral gas. The Lyman-α forest probes the matter power spectrum in the mildly nonlinear regime over a broad range of redshifts, down to scales of $1–80$ h/Mpc (at $z = 2–6$ in ground-based

observation) [2331–2333]. Compared to WL surveys, the Lyman-α forest is sensitive to smaller scale structures at higher redshifts. Since the Lyman-α forest mainly traces the IGM in the cosmic web, directly correlated with the neutral hydrogen gas distribution, its signal is highly influenced by the baryonic physics [2334]. High-precision hydrodynamical simulations that accurately account for baryonic physics are essential for studying the Lyman-α forest. Using these simulation results, researchers can investigate physics beyond the ΛCDM, including WDM [1550,1553,2333,2335], fuzzy DM [1551,1563,2336] and neutrino mass [2337,2338]. Recent simulations suggest that small tensions may exist between ΛCDM cosmological simulation and Lyman-α observation (see Section 2.3.8 for the anomalies in Lyman-α measurements), which can be explained by mechanisms such as the dark photon heating or DM-neutrino interaction [2339,2340]. An improved understanding of baryonic physics is crucial to addressing these issues and refining our interpretation of Lyman-α observation.

### 3.5.4.3. Cosmic 21 cm.
The cosmic 21 cm signal refers to the hyperfine transition of neutral hydrogen atoms, which emits or absorbs radiation at a wavelength of 21 cm (frequency of $1.42$ GHz). This signal is a powerful tool for probing the evolution of the Universe at high redshift, across a wide range of cosmic epochs, from the Dark Ages (prior to star formation) to the EoR and beyond [2341,2342]. The 21 cm signal traces the distribution of matter and can map the evolution of large-scale structures over cosmic time. Although the angular resolution of the 21 cm signal is limited to tracing ionized bubbles on scales of several Mpc, which corresponds to very large structures, it remains a powerful tool for exploring the nonlinear scale structure. The small scale structure significantly impacts the halo collapsing. In general, any model that suppresses or modifies the amplitude of DM fluctuations on small scales could affect the 21 cm cosmic dawn signal. Several previous works have discussed using the future cosmic 21 cm data to probe the physics beyond ΛCDM [2343–2348]. One key advantage of the cosmic 21 cm signal is its ability to trace the mass distribution at very high redshifts, where the original information is preserved due to limited impact from structure formation [2311,2349]. By modeling the dynamics of DM and the feedback of baryons, cosmological simulations play a pivotal role in understanding and interpreting the cosmic 21 cm signal [2311,2348,2350]. Upcoming cosmic 21 cm observations, such as the Low-Frequency Array (LOFAR) [2351], the Hydrogen Epoch of Reionization Array (HERA) [2352], and the SKAO [2353], are expected to open a new window into the high redshift Universe and provide novel insights into the nature of DM.

### 3.5.4.4. Dwarf galaxy scale puzzles.
At small scales, such as those of dwarf galaxies, significant tensions exist between predictions from N-body simulations and observational data [2354]. Several key discrepancies are summarized below:

**Core-cusp Problem** Gravity only N-body simulations of DM in a CDM Universe predict that DM halos should have a “cusp” — a sharp rise in density toward their center, known as Navarro-Frenk-White profile, or NFW [2355]. However, observations of low-mass galaxies of different luminosity and morphology [358] indicate that their DM halos instead have a core, a central region with nearly constant density [2356,2357].

**Diversity Problem for Rotation Curves** Simulations predict that galaxies with the same maximum circular velocity should have similar halo density profiles. However, observations reveal a significant variation in the central densities of these galaxies [2358, 2359].

**Missing Satellites Problem** Simulations predict a large number of small DM halos, but observations show a much smaller number of dwarf satellite galaxies around larger galaxies, such as the MW [2360,2361].
**Too-Big-to-Fail Problem** The number of massive DM subhalos predicted by simulations is inconsistent with the number of observed bright satellite galaxies in systems like the MW [2362, 2363].

These small-scale issues highlight a disconnect between N-body simulations and observations, particularly in our local group. They suggest that gravity-only CDM N-body simulations are insufficient to fully describe the Universe on these scales. Several hypotheses have been proposed to resolve these discrepancies, involving both astrophysical processes and modifications to the standard ΛCDM model. Stellar and supernova feedback can redistribute baryonic matter, influencing the gravitational potential and altering the density profile of the DM and the halo mass function [2364,2365]. The feedback process can also eject gas from massive subhalos, suppressing star formation to explain the too big to fail problem [2366–2370]. Another way to solve this small-scale puzzle requires modifications to ΛCDM. Several models have been proposed such as Self-Interacting DM (SIDM) [2082] and fuzzy DM [2371]. SIDM introduces DM-DM interactions that redistribute energy, forming cores in halos, and leading to more diverse halo density profiles and a different subhalo mass function. N-body simulations incorporating SIDM have been used to test these predictions and compare them with the observations. Many SIDM simulations were performed and used to discuss the DM halo distribution after the first proposal [2372–2377]. Early studies constrained the SIDM cross section through observations, such as strong lensing [2378]. The results show that the parameter space for the SIDM model to solve the small scale problem have already been constrained. However, subsequent simulations with higher resolution and improved algorithms revealed that earlier constraints may have been overestimated [2379–2382]. For fuzzy DM, the introduction of ultra-light bosons with kiloparsec-scale wavelength naturally suppresses structure formation at these scales. Unlike classical particles, fuzzy DM behaves more like a wave, requiring the solution of quantum-mechanical wave equations in N-body simulations. Studies have shown that fuzzy DM naturally forms solitonic cores at the centers of halos due to quantum pressure while simultaneously suppressing small-scale structure formation as a consequence of the uncertainty principle [2383–2385]. This provides a compelling explanation for the core-cusp problem and the missing satellites problem. For now, resolving small-scale issues remains an ongoing challenge, as both baryonic feedback processes and alternative DM properties beyond CDM have their own strengths and limitations. High-resolution cosmological simulations, combined with advanced algorithms to incorporate baryonic physics and additional DM interactions, provide a crucial tool for understanding these small-scale puzzles. These simulations become a bridge between the theoretical predictions from proposed models and the observational data. They help us better understand the existing cosmological tensions.

## 3.6. Profile likelihoods in cosmology

**Coordinator:** Adrià Gómez-Valent

**Contributors:** Elisa G. M. Ferreira, Emil Brinch Holm, Giacomo Galoni, Laura Herold, Paolo Campeti, Sophie Henrot-Versillé, Vivian Poulin, and William Giarè

### 3.6.1. Motivation

In cosmology, it is very common to encounter high-dimensional parameter spaces. Efficient exploration of these spaces through the sampling of their corresponding multivariate distributions is commonly achieved using Monte Carlo techniques. In these analyses, statistical information is contained in the resulting MCMC, from which one can infer the posterior distribution and other derived Bayesian products. In parameter spaces of dimension greater than two one cannot visualize the posterior distribution in the original parameter space. One needs

to apply the so-called marginalization technique to project the results into spaces of lower dimension and compute the marginalized posterior distributions of the individual parameters and their constraints at the desired confidence level. This is also done to obtain the contour plots in all the two-dimensional planes of interest. Marginalization is extremely efficient and is a very practical tool that eases the visualization and interpretation of the results in Bayesian analyses. However, when the original (non-marginalized) PDF has non-Gaussian features, the marginalized posterior distributions obtained do not exclusively reflect how the parameters are distributed based on their ability to explain the data but also according to their integrated probability weight. This introduces volume (or marginalization) effects, which can potentially bias conclusions in a significant way. In other words, if deviations from Gaussianity are large, the marginalized posterior can hide points in parameter space that fit the data well. Another potential drawback of the marginalization method is the impact of the prior, which, even if assumed to be flat, is still informative in all cases. For instance, a flat prior on a variable x or on its logarithm, log x, might induce important changes in the constraints of a model, even if the boundaries of the priors are perfectly consistent. This is due to volume effects, again. Hence, volume effects can significantly impact the interpretation of results extracted from MCMCs. Detecting these biases is of utmost importance. In the context of cosmology, accurately quantifying these biases is essential, e.g., for a correct assessment of cosmological tensions in the ΛCDM model and beyond, as well as in data compression performed by galaxy and WL surveys, see Section 3.6.4 for details.

A mismatch between the point in parameter space that maximizes the original posterior distribution and the point constructed with the values of the parameters that maximize the individual one-dimensional posteriors already hints at the non-negligible impact of marginalization effects. One can investigate the impact of marginalization effects using the profile likelihood (PL). The PL is a method from frequentist statistics to construct parameter intervals, and therefore, by construction inherently independent of a prior. It ensures, in particular, that the results (on the parameter inference) do not depend on the choice of the parameterization of the problem. This part of the work aims to explain the basics of PLs and is organized as follows. In Section 3.6.2 we describe how to build confidence intervals using PLs. In Section 3.6.3 we describe several methods to compute them, both approximately from the MCMCs or, more precisely, using advanced numerical optimization techniques with the help of concrete codes. Finally, in Section 3.6.4 we give some examples of applications of the PLs in cosmology, providing a rich list of references.

### 3.6.2. Confidence intervals from profile likelihoods

While Bayesian credible intervals are statements about the degree of belief of the value of the true parameter, frequentist confidence intervals are based on the coverage: frequentist confidence intervals are said to cover the true value of the parameter with a confidence level α if — after repetition of the experiment — a fraction α of the experiments contain the true value of the parameter (e.g., see Ref. [1467] for a review, or Ref. [2386] in the context of cosmology). Confidence intervals with correct coverage can be constructed using the Neyman construction [2387]. This construction requires knowledge of the full PDF P(μ|μ<sub>true</sub>), where μ are the maximum likelihood estimates (MLE) of the model parameters μ estimated from observations and μ<sub>true</sub> the (hypothetical) true values of these parameters. Since this probability is typically not known a priori, it needs to be estimated by generating mock realizations of the data. This makes the Neyman construction computationally expensive. Fortunately, a theorem by Wilks [2388] helps in many practical applications, which states: In the limit of a large data set, the log-likelihood ratio

$$ \Delta\chi^2 = -2 \log R(x, \mu), \qquad R(x, \mu) = \frac{\mathcal{L}(x|\mu)}{\mathcal{L}(x|\hat{\mu})}, \qquad (3.29) $$
follows a $\chi^2$-distribution. Here $\mathcal{L}$ is the likelihood function and $x$ is the data vector. In this case, the Neyman construction can be replaced by a much simpler graphical construction based on the profile likelihood. Wilks theorem holds trivially if the p.d.f. is Gaussian (for every possible value of the true parameters). We note, however, that Wilks' theorem holds for the full-dimensional likelihood (for any possible value of the true parameters) and not the profile likelihood. Although this effect vanishes in the large data limit, it is often very difficult to assess whether this limit is reached in practice [2389], and the confidence intervals are, therefore, only approximate (see Ref. [2386] for checks of the asymptotic assumptions in cosmology).

The profile likelihood for a parameter of interest $\mu$ (in one dimension) can be obtained by evaluating Eq. (3.29) for different values of $\mu$. In the presence of other cosmological parameters and nuisance parameters, $\nu$, the profile likelihood is calculated by setting $\mu$ to various values within a range of interest and minimizing $\chi^2(\mu, \nu) = -2 \log \mathcal{L}(\mu, \nu)$ with respect to all parameters $\nu$. The global MLE, or “best-fit”, corresponds to the minimum value $\chi^2_{\min}$ by design. A graphical frequentist confidence interval for $\mu$ can now be constructed using the difference $\Delta\chi^2(\mu) = \chi^2(\mu) - \chi^2_{\min}$. When $\mu$ is far from its physical boundary, a confidence interval at confidence level $\alpha$ is obtained by applying a fixed threshold $\Delta\chi^2_{\text{th}}$ to $\Delta\chi^2(\mu)$. This threshold is chosen such that the cumulative distribution function of the $\chi^2$ distribution with one degree of freedom equals $\alpha$. Notable values are $\Delta\chi^2_{\text{th}} = 1$ for 68% CL and $\Delta\chi^2_{\text{th}} = 3.84$ for 95% CL [2390]. This method is applicable to both parabolic (i.e., for a Gaussian-distributed parameter) and non-parabolic $\Delta\chi^2(\mu)$ because of the invariance of the MLE under reparameterization.

When the parameter estimate is near its physical boundary, the conventional graphical profile likelihood construction for frequentist confidence intervals can be inadequate. It may result in empty intervals and fail to maintain the frequentist coverage property if the choice between an upper limit or a two-sided interval is made based on the data (the so-called “flip-flopping”). An adapted Neyman construction, also known as Feldman-Cousins (FC) prescription [2391] in the particle-physics literature, addresses these issues. To describe this method, we restrict the discussion to one-dimensional data $x$ and a single model parameter $\mu$ for clarity. For each value of $\mu$ (with an unknown true value) and each observable $x$, we compute the likelihood ratio $R(x, \mu)$, with $\mu$ (and its MLE $\hat{\mu} \equiv \mu_{\text{best}}$) restricted to its physically allowed values. The confidence belt at the desired CL $\alpha$ is constructed by selecting an acceptance interval $[x_1, x_2]$ for each $\mu$ satisfying

$$ \begin{cases} R(x_1, \mu) = R(x_2, \mu), \\ \int_{x_1}^{x_2} P(x|\mu) dx = \alpha, \end{cases} \qquad (3.30) $$

where $P(x|\mu)$ is the p.d.f. of $x$ given $\mu$, and values of $x$ are added to the acceptance interval in order of decreasing likelihood ratio. The confidence belt is formed by the union of all acceptance intervals $[x_1(\mu), x_2(\mu)]$. Intersecting the confidence belt with a line at $x = x_0$, where $x_0$ is the value of $x$ minimizing $\chi^2$ (i.e., the value observed in an experiment), yields the confidence interval $[\mu_1, \mu_2]$ for $\mu$. The adapted Neyman or FC prescription provides a method to determine the endpoints of confidence intervals, which smoothly transitions between an upper limit and a two-sided interval, ensuring exact frequentist coverage even for parameters with non-Gaussian p.d.f.'s (and if Wilks' theorem does not hold). This is in contrast to the conservatism (over-coverage) inherent in Bayesian limits in the same context [2391,2392]. While overcoverage might not be as problematic as undercoverage, it reduces the ability to reject false hypotheses. This highlights the value of examining frequentist intervals.

If the p.d.f. of the parameter is Gaussian and the MLE of the parameter near a physical boundary at $\mu = 0$, then $\mu_{\text{best}} = \max(0, x)$ and the likelihood ratio in simply becomes [2391]

$$ R(x, \mu) = \begin{cases} \exp(-(x - \mu)^2/2), & \text{for } x \ge 0, \\ \exp(x\mu - \mu^2/2), & \text{for } x < 0, \end{cases} \qquad (3.31) $$

where $x$ and $\mu$ are in units of $\sigma$, the width of the parabolic fit to $\Delta\chi^2(\mu)$. Then, the confidence interval can be derived by solving Eq. (3.30), with $P(x|\mu)$ being a Gaussian with mean $\mu$ and unit variance. If, however, the p.d.f. of the parameter is non-Gaussian and Wilks' theorem does not hold, a full Neyman construction of the confidence intervals using simulations for several input values of the parameter of interest is necessary, e.g., see Ref. [2393,2394] for details.

With the necessary modifications, this entire treatment can be generalized to multiple dimensions. For example, under the assumption of Gaussianity, the iso-$\chi^2$ curves for a two-dimensional case will be ellipses on the parameter space, whose angle of the major axis depends on the correlation of the two parameters of interest. Furthermore, if no physical boundary is close to the absolute minimum of $\chi^2$ and if we consider the limit of a large data set, it is possible to adopt a graphical approach to derive the 68% and 95% confidence regions, as we mentioned for the one-dimensional case. This time, one must search for $\Delta\chi^2_{\text{th}} = 2.3$ for 68% CL and $\Delta\chi^2_{\text{th}} = 5.99$ for 95% CL regions. Similarly to the one-dimensional case, if Wilks' theorem does not hold, one should adopt a full Neyman construction to extract the confidence regions.

### 3.6.3. Calculation of profile likelihoods

Computing a profile likelihood for a single parameter involves a series of M optimizations, typically with $M \sim O(10)$, in the $N-1$ dimensional parameter space. Furthermore, these cannot be reused for profiles in other parameters, unlike an MCMC, which simultaneously produces parameter constraints for all varied parameters. Since the optimizations themselves are moderately difficult, constructing a full profile likelihood analysis of several parameters can be very computationally demanding. In this subsection, we present some of the strategies adopted in the literature to accommodate the often high computational cost of computing profile likelihoods.

#### 3.6.3.1. Optimization strategies.
A converged set of MCMC chains of course contains information about the likelihood values at a large number of points in parameter space. Although MCMC algorithms are poor optimizers [1974], they can be used to estimate profile likelihoods by binning their points along the parameter $\theta$ being profiled. The bin around a fixed value of $\theta = \theta_0$ can be defined, for example, by a homogeneous central binning, as is done in Refs. [1944,1976,2395–2397], or by taking the $n$ points of the MCMC that have a value of $\theta$ closest to the fixed value $\theta_0$, as done in Ref. [2398], with $n$ a fixed fraction of the total amount of MCMC points. Due to the finite sampling of the MCMC, the largest likelihood value found in each bin will generally be smaller than its optimized value at $\theta_0$, which leads to noise in the profile likelihood curve. The great advantage of this approach, however, is that it has no computational cost if an MCMC is already given, so it can be used as an inexpensive initial test to check for marginalization and prior effects [1976], see also Ref. [2396]. Furthermore, such a profile can be used as the starting point for explicit optimizations at each fixed value of the profiled parameter.

When the binned MCMC is not adequately accurate, the profile likelihood can be computed by explicit optimization in the parameter subspace with the profiled parameter taking on a series of fixed values. Generally speaking, the most efficient optimization algorithms explore the parameter space for the highest likelihood value by using the gradient of the likelihood with respect to parameter space. Examples of such gradient-based optimizations algorithms often used in cosmology include the MIGRAD optimizer of the MINUIT package, which is a variable metric algorithm using the finite-difference estimates of the first derivative of the likelihood function [2399]. Refs. [2400–2404], amongst others, used MINUIT to construct profile likelihoods and found it to be computationally efficient.

In some cases, however, the cosmological likelihood functions are noisy, for example, due to approximation switching in the underlying theory code or insufficient precision settings, leading to inefficiency
in the gradient-based algorithms. For example, Refs. [710,1969,2397,2405] found that gradient-free methods outperformed the gradient-based ones. The most popular gradient-free method employed in the literature is simulated annealing [2406], which works by running an MCMC using the modified likelihood $\tilde{L} \equiv L^{1/T}$ with an iteratively decreasing temperature $T$. The decreasing temperature enhances the peak structures in the likelihoods and eventually traps the MCMC close to the best fit. Simulated annealing has been used for profile likelihood construction in Refs. [1959,1967–1970,1972,1977,2397,2398], amongst others. Unfortunately, simulated annealing is very sensitive to its hyperparameters, which include the initial point of the optimization, the covariance matrix of the proposal distribution of the MCMC, and the rate at which the temperature decreases. Broadly speaking, three approaches have been suggested to inform the initial point and covariance matrix:

* *Running fixed-parameter MCMCs to obtain information.* Refs. [1967–1969] run individual MCMCs for each fixed value of the profiled parameter and use the best fit and covariance matrices of their chains in the optimization. This method produces great initial points and local covariance matrices, although with the disadvantage of the added computational cost of running the additional MCMCs.

* *Using information from an existing MCMC.* The CAMEL [2407] and PROSPECT codes [2398] assume the user to have a converged MCMC in the full parameter space before constructing the profile. Using this, it uses the binned MCMC profile likelihood estimate discussed in the preceding subsection, and furthermore constructs individual covariance matrices from the points in each bin. This method usually gives great initial points and local covariances matrices, but can be sensitive to the quality of the initial MCMC.

* *Optimizing the profile sequentially and reusing information from previous optimizations.* The PROCOLI code [2397] also assumes the availability of a converged MCMC at the beginning and uses its maximum a posteriori point and covariance matrix to find the global best fit. PROCOLI then carries out optimizations at equidistant points along the profiled parameter direction sequentially, using the result of the previous iteration as the initial point of the next. This method yields great initial points, at the sacrifice of the parallelization of the optimizations.

### 3.6.3.2. Publicly available profile likelihood codes.
At the time of writing, there exist several publicly available tools for computing profile likelihoods in cosmology; here, we present them in order of the recency of their release.

* COBAYA FORK [2400]: Can be used with COBAYA [1918]. Initializes points from an existing MCMC as described above. Usable as a sampler inside COBAYA, making it easy to use if an MCMC has already been run with COBAYA.<fnref n="43" />

* PROCOLI [2397]: Can be used with MONTEPYTHON [1916,1917]. Uses the sequential strategy described above. Additionally allows for the extraction of individual $\chi^2$ values for each experiment used in the data.<fnref n="44" />

* PROSPECT [2398]: Can be used with MONTEPYTHON and COBAYA. Uses an assumed precomputed MCMC as initialization, as described in the last section. Introduces an automatic step size tuning scheme, and comes with a snapshot checkpoint scheme that enables efficient parallelization and inspections of the run.<fnref n="45" />

* CAMEL [2407]: Code written in C++. The minimization algorithm is based on MINUIT [2399].<fnref n="46" />

All of the above codes (except CAMEL) work directly with the cosmology-specifying parameter files of MONTEPYTHON or COBAYA, respectively, that have been used, for example, in a complementary Bayesian analysis. They therefore also work with the theory codes already available in these, such as CLASS [1949] and, in the case of COBAYA, CAMB [1948]. They also employ the GETDIST code [1952] for analysis.

While the brute-force computation of full triangle plots of profile likelihoods, analogous to the most common product of the Bayesian analysis, is computationally unfeasible, future advances in gradient-based inference and emulators of likelihood codes are expected to make extensive profile likelihood analyses more broadly available. For example, already available is the neural network emulator CON-nect<fnref n="47" /> [2408] has a run-mode that efficiently computes full triangle plots of profile likelihoods using gradient-based optimization, as explored in Ref. [2409].

### 3.6.4. Applications of profile likelihood in cosmology

* **EDE & NEDE:** Due to the complicated parameter structure of EDE (see Section 4.1.1), volume effects can have a strong impact on the constraints of this model. Using PLs, Refs. [1967,1968] showed that tight upper limits on EDE from CMB and LSS in Bayesian analyses, which questioned the ability of EDE to resolve the Hubble tension, are partially driven by prior volume effects, while in a frequentist analysis EDE presents a viable solution to the tension. Ref. [1969] uses PLs to explore the interplay between EDE and massive neutrinos. In Ref. [1970], both a Bayesian and PL analysis under an updated Planck CMB pipeline give tight upper limits on the fraction of EDE. Moreover, Ref. [1972] constrained the NEDE (see Section 4.1.2) model and showed that the frequentist analysis gives greater evidence of the model due to volume effects in the Bayesian posterior associated with the $\Lambda$CDM limit.

* **EFTofLSS:** Ref. [1959] studied the impacts of priors on the nuisance parameters of the EFTofLSS analysis method applied to full-shape power spectrum from galaxy surveys. They found that the priors usually employed to constrain the parameters to their physical regime were informative and had a significant effect on the constraining power of the analysis method.

* **Decaying dark matter:** Ref. [1977] studied a model where a fraction of the DM is allowed to decay into dark radiation. With two new parameters, the abundance of the decaying part of the CDM and its lifetime, series volume effects are found in the limits of parameter space that recover the $\Lambda$CDM model; these include the case of infinite lifetime and vanishing abundance. In particular, whereas the Bayesian analyses prefer either of these limits, [1977] shows that the best fit is associated with a narrow log-likelihood peak at a lifetime corresponding to decays just around recombination.

* **Phenomenological transition of dark matter to dark radiation:** Ref. [2398] tested the PROSPECT code on a model where a fraction of the DM transitions to dark radiation, characterized by three parameters: The fraction of DM that transitions, the central redshift of the transition, and a parameter that controls the shape of the transition [2410]. The $\Lambda$CDM limit is recovered in many limits of parameter space, leading to significant volume effects and parameterization dependencies that were diagnosed with profile likelihoods.

* **Lensing:** A complete study of the so-called $A_L$ tension in the $\Lambda$CDM + $A_L$ has been done in Ref. [727].

* **Neutrino:** One of the prototypical applications of profile likelihood in cosmology is for neutrinos. Since neutrino mass constraints are close to the boundary $\sum m_\nu > 0$ and the non-Gaussian posteriors

<footnote n="43">43 https://github.com/ggalloni/cobaya/tree/profile_sampler</footnote>

<footnote n="44">44 https://github.com/tkarwal/procoli</footnote>

<footnote n="45">45 https://github.com/AarhusCosmology/prospect_public</footnote>

<footnote n="46">46 http://camel.in2p3.fr/wiki/pmwiki.php</footnote>

<footnote n="47">47 https://github.com/AarhusCosmology/connect_public</footnote>
found MCMC data analysis, profile likelihood was used to obtain the confidence interval on the sum of the neutrino mass. This has been done in the context of WMAP5, Sloan Digital Sky Survey (SDSS), and HST in Ref. [2395] and for Planck data in Ref. [2404]. In Ref. [2411] and Ref. [2412], they study the effect of priors in constraining the number of neutrino species using cosmological data and no evidence for deviations from the standard number of neutrino using profile likelihood. In Ref. [2411] was also highlighted the link with $\Lambda_{\text{lens}}$ on the current limit. Refs. [1286,2413] explore the preference for “negative” neutrino masses with profile likelihoods.

*   **Number of relativistic species, $N_{\text{eff}}$:** Profile likelihood was used to test the possibility that extra relativistic species are present beyond the standard model prediction. A discrepancy between the posterior and profile likelihood confidence interval was found in $N_{\text{eff}}$ in Refs. [1973,1974], indicating volume effects when using Wilkinson Microwave Anisotropy Probe (WMAP7) and HST data. It was applied to Planck data in Ref. [1975].

*   **Gravitational wave energy density:** Constraints on the primordial GW density using CMB data and a study of the impact on cosmic string models has also been derived in Ref. [2414].

*   **Tensor-to-scalar ratio:** Profile likelihood was used to analyze the tensor-to-scalar ratio, $r$, with current data. In Ref. [2403] the most recent Planck and BICEP/Keck CMB data was used in order to gain a clearer perspective on the discrepancy between Bayesian and frequentist upper limits on $r$ previously found by the SPIDER collaboration [2394]. A similar analysis was done in Ref. [2400] combining CMB data from Planck PR4 and BICEP/Keck, and GW data from LIGO-Virgo-KAGRA. With the goal of testing the impact of MG theories on $r$, in Ref. [2415] combining data from CMB from Planck DR4 and BICEP/Keck and BAO data from 6dF Galaxy Survey, SDSS DR7 and eBOSS, the Bayesian and frequentist upper limits on $r$ were obtained.

*   **Inflation:** The advantages of using a profile likelihood, and more in general, of constructing frequentist confidence intervals with the FC prescription, are particularly evident when trying to constrain inflationary models. This typically allows for more informative constraints and physics insights, especially on models with numerous additional parameters, which are typically degenerate and hard to constrain simultaneously, especially with current data. An example is provided by the latest constraints from Planck and BICEP/Keck on axion-U(1) inflation [2402] and by forecasted constraints on axion-SU(2) inflation parameters for the future CMB satellite survey LiteBIRD [2393]. Another example is given by the attempt to constrain the spectral tilt of tensor primordial perturbations with the latest CMB and GW interferometers data [2400]. Profile likelihoods were also used to search for features in the primordial power spectrum imprinted in CMB and large-scale structure formation observables, using data from Planck and WiggleZ [2416].

*   **Brans-Dicke model with cosmological constant:** The Brans-Dicke model with a constant positive vacuum energy density, the so-called BD-ACDM model, has been confronted with observations multiple times in the last years [2417–2421]. In Ref. [2420] the authors showed that in the absence of the Planck 2018 high-$\ell$ CMB polarization and lensing data, but still under a very rich dataset including the Planck 2018 full temperature and low-$\ell$ polarization likelihoods, together with the state-of-the-art data on SNIa, BAO, CC and RSD, it is possible to loosen the $H_0$ tension in BD-ACDM, while keeping $\sigma_8 \sim 0.790$. In Ref. [1976] the author tested the impact of volume effects in BD-ACDM through the analysis of the profile distribution and showed that they do not play a major role in this model.

*   **Coupled quintessence:** Coupled quintessence was firstly studied in Refs. [2051,2422] considering a scalar-mediated interaction (fifth

force) between DM particles, with the scalar field playing the role of DE. Constraints on this model have been reported e.g., in Refs. [2423–2429]. Volume effects have been recently quantified in Ref. [1976] for the case of a Peebles-Ratra potential [295,2430] using the profile distributions built directly from the MCMC. Their impact on the results are found to be small.

*   **Neutrino-Dark Matter Interactions:** Ref. [2396] studied a model featuring scatter-like interactions between neutrinos and DM. Allowing the interaction strength and the total neutrino mass to vary, the model can have up to two parameters more than $\Lambda$CDM. A profile likelihood analysis was employed to test possible volume effects. In this case, no relevant volume effects were found. The profile likelihood supported the results derived by marginalizing over the parameter space, confirming a mild preference for non-vanishing interactions when small-scale CMB data are considered [2431,2432].

*   **Baryon Acoustic Oscillations:** Profile likelihood has been used, sometimes in combination with Bayesian methods, to quantify BAO errors in determining the fitting scale in surveys. This has been done in the context of many surveys, like the Sloan Digital Sky Survey (SDSS) BOSS [2433] and the eBOSS [2434], or DES [2435,2436]. In [2437], they extend this method to add noise and compare different approaches to determine the error in determining the BAO scale. Profile likelihood was also used to fit the Lyman-$\alpha$ (Ly$\alpha$) forest 3D correlation function [2438]. Using mocks, they found that profile likelihood is a good approximation for fitting the BAO peak from the Ly$\alpha$ forest correlation functions and is in full agreement with the Bayesian analysis (which is not the case for frequentist maximum likelihood estimators, which assume Gaussianity).

# 4. Fundamental physics

**Coordinator:** Vivian Poulin

Cosmology stands at a crossroads. The concordance $\Lambda$CDM model, which had shown a great level of success in describing a wide variety of cosmological observations, is being challenged by high precision data as reviewed in Section 2. However, the $\Lambda$CDM model has raised profound questions about the nature of its dominant components –DM and DE– as well as the mechanism at the origin of the primordial fluctuations (usually described by inflation) that remains largely unknown. Barring the existence of unknown systematic errors, these observational issues thus present a crucial opportunity to explore new physics beyond the $\Lambda$CDM paradigm. This section examines the collective efforts of the scientific community to interpret the implications of observations that challenge the $\Lambda$CDM model, with a particular focus on what the $H_0$ and $S_8$ tensions may reveal about fundamental physics.

At the core of the $\Lambda$CDM value of $H_0$ is the angular size of the sound horizon at recombination $\theta_s(z_*) = r_s(z_*)/D_\Lambda(z_*)$, that CMB data have determined at $\mathcal{O} \sim 0.1\%$ precision [192]. In the flat FLRW metric, the angular diameter distance to recombination $D_\Lambda(z_*)$ that carries information about $H_0$ can be computed as

$$D_\Lambda(z_*) = (1 + z_*)^{-1} \int_0^{z_*} \frac{dz}{H(z)}, \qquad (4.1)$$

where under flat $\Lambda$CDM, $H(z) = H_0 \sqrt{\Omega_r(1+z)^4 + \Omega_{m,0}(1+z)^3 + \Omega_\Lambda}$. Extracting the angular diameter distance from $\theta_s(z_*)$ requires knowledge of the sound horizon

$$r_s(z_*) = \int_{z_*}^{\infty} \frac{c_s(z) dz}{H(z)}, \quad \text{with} \quad c_s(z) = \frac{1}{\sqrt{3(1 + 4\rho_b(z)/3\rho_r(z))}}. \qquad (4.2)$$

A value of $r_s(z_*)$ can be computed within a model, given knowledge of the cosmological parameters. Within $\Lambda$CDM, Planck data have allowed us to determine $r_s(z_*) = 144.39 \pm 0.3$ Mpc (TTTEEE+lowE) from the complex structure of peaks and troughs [192]. The value of $H_0$ is then
adjusted in order to match the corresponding value of $\theta_s(z_*)$, which yields $H_0 = 67.27 \pm 0.6 \text{ km s}^{-1} \text{Mpc}^{-1}$. Similarly, the value of $S_8$ can be inferred from the linear matter power spectrum $P(k)$ given a set of cosmological parameters reconstructed from a fit to a given dataset as

$$S_8 \equiv \sqrt{\frac{\Omega_{m,0}}{0.3}} \times \left\{ \sigma_8 \equiv \int_0^\infty \frac{k^3}{2\pi^3} P(k) W_8^2(k) d \ln k \right\}, \quad (4.3)$$

where $W_8(k)$ is a window function describing a sphere of radius $R = 8 \text{Mpc/h}$ in Fourier space. Under the $\Lambda$CDM model fit to Planck data, it yields $S_8 = 0.834 \pm 0.016$.

The solutions to the $H_0$ discrepancy generally fall into two categories based on the requirement to leave the precisely measured value of $\theta_s(z_*)$ unaffected [709,710,713]. First, there are those that reduce the value of the sound horizon $r_s(z_*)$ such as to compensate for the shorter angular diameter distance implied by a larger $H_0$. Those typically correspond to early-time (pre-recombination) modifications and include for instance include, new relativistic species, EDE, or modified recombination scenarios. Second, late-time solutions modify the expansion history after recombination such that the angular diameter distance to recombination remains unaffected. These include e.g., phantom DE, IDE, or decaying DM. Both avenues offer advantages and disadvantages: the former scenarios may more easily adjust the shape of the expansion history that is constrained by BAO and SNIa (beyond the mere value of $H_0$), but will be strongly constrained by CMB observations, while the latter scenarios can more easily leave the CMB unaffected by exploiting the geometrical degeneracy, but are strongly constrained by late-time observations [710].

Resolving the $S_8$ tension, on the other hand, requires suppressing structure growth, which can be done through new DM interactions, evolving DE models, or modifications to general relativity. Yet, it remains to be established precisely what scales and what cosmic era must be affected in order to achieve concordance between all our LSS observations: is the tension solely appearing at small scales, that are potentially affected by large uncertainties related to non-linear and baryonic effects, or does it also appear on large-scales at late-times, where the physics is much more linear and under control? Answering these questions will drive the community effort to find a solution to the tension. Moreover, there is a clear interplay between the $H_0$ and $S_8$ tensions. First, simply because the definition of the physical scale $k$ involved in $S_8$ in Eq. (4.3) is trivially influenced by the value of $h$: as the power spectrum as a strong dependence on the scale $k$ (falling like $k^{-3}$) tension in $S_8$ may appear due to different values of $h$ in the analysis, fooling us in comparing different physical scales [2439–2441]. However, even when corrected (e.g., through the use of $R = 12 \text{ Mpc}$ as an absolute scale [2439]), the tension can be influenced by the new physics introduced to resolve the Hubble tension. This stems from the observation that distance indicators, that provide measurements of $\Omega_{m,0}$ and $H_0$, necessarily imply a measurement of the physical matter density $\omega_{m,0} \equiv \Omega_{m,0}h^2$ [804,2442,2443]. Hence, a larger value of $h$ will automatically imply a larger value of $\omega_{m,0}$, leading to earlier matter domination and thus a larger amplitude of fluctuations today.

Finding a solution to either or both tensions is an active area of research. As we illustrate within each section, understanding and resolving these tensions will not only improve the precision in measuring fundamental cosmological parameters but also offer a window into the fundamental physics governing the Universe [709]. Whether these tensions indicate systematic biases in data or genuine hints of new physics, they represent one of the most exciting frontiers in modern cosmology, potentially unlocking new insights into DM, DE, and the laws of gravity.

## 4.1. Early-time proposals

### 4.1.1. Early dark energy and variants

**Coordinator:** Laura Herold

**Contributors:** Adrià Gómez-Valent, Alexander Reeves, Alireza Talebian, David Benisty, Elisa G. M. Ferreira, Gabriela Garcia-Arroyo, Ivonne Zavala, Joan Solà Peracaula, Lu Yin, Luis Anchordoqui, Luz Ángela García, Matteo Forconi, Nikolaos Mavromatos, Rafaela Gsponer, and Vivian Poulin

#### 4.1.1.1. Introduction and motivation.
EDE models (for reviews see Refs. [2444–2446]) are a promising class of models proposed to resolve the Hubble tension with new physics. The central idea behind EDE is to introduce an additional component to the energy density of the Universe at early-times, which increases the expansion rate $H(z)$ before recombination $z^*$. This leads to a decrease of the (comoving) size of the sound horizon: $r_s(z^*) = \int_{z^*}^\infty c_s(z)/H(z)dz$, where $c_s$ is the sound speed in the baryon-photon plasma. To understand why this increases $H_0$ at present times, note that the angular size of the sound horizon $\theta_s = r_s/D_\Lambda$ is directly and precisely observed by the CMB, e.g., see Ref. [192]. Hence, to maintain the same size of $\theta_s$, the (comoving) angular diameter distance reduces, $D_\Lambda(z^*) = \int_0^{z^*} dz/H(z)$, and leads to an increase of $H_0$.

All EDE-type solutions to the Hubble tension have in common a DE-like rapid expansion with an equation of state $w \simeq -1$ around $z^*$, which peaks at some critical redshift $z_c$ and is followed by a phase in which the background energy density decays or dilutes (see Fig. 70, left). EDE models often take the form of a cosmological scalar, whose background dynamics follow the homogeneous Klein–Gordon equation: the field is initially frozen in its potential, such that the background energy density is constant. The fractional contribution of the field to the total energy density, $\rho_{\text{EDE}}/\rho_{\text{Tot}}$, increases over time in a manner similar to that of DE. Eventually, Hubble friction dropping below a critical value, or a phase transition changing the shape of the potential, releases the scalar, at which point the field becomes dynamical, and the background energy density dilutes away faster than matter. The contribution of EDE to the Hubble rate is therefore localized before recombination, where it efficiently acts to reduce the size of the sound horizon.

EDE-like models have already been studied before the appearance of the Hubble tension, e.g., in the context of quintessence models [2447–2450] or in the context of recurring DE [2451]. The (axion-like) EDE model has gained significant attention as a proposed solution to the Hubble tension [2452–2454]. To fully resolve the tension, a fractional EDE contribution of $\sim 10\%$ around the time of matter–radiation equality is required [2455], reflecting the percentage difference in $H_0$ estimates between Planck [192] and SH0ES [34]. Although EDE models present a promising class of solutions when compared to other proposed solutions [710,713,2456], the improvement in observational constraints from the CMB and LSS in the last few years poses challenges to this class of models.

#### 4.1.1.2. The axion-like EDE model.
The originally proposed EDE model [2452–2454] consists of a scalar field $\phi$ in a potential of the form

$$V(\phi) = V_0[1 - \cos(\phi/f)]^n, \quad (4.4)$$

where $V_0 = m^2 f^2$ with $m$ and $f$ being referred to as “mass” and “decay constant”, respectively. The index $n$ is typically fixed to $n = 3$ as this provides the best fit to the data [2445,2455]. The parameters $m$ and $f$, along with the initial value of the scalar field in the potential, $\theta_i = \phi_i/f$, need to be fit from the data. This “particle-physics parameterization” is often traded for the “phenomenological parameterization”, which consists of $f_{\text{EDE}}$, the maximum fraction of EDE, $f_{\text{EDE}} \equiv \rho_{\text{EDE}}(z_c)/\rho_{\text{Tot}}(z_c)$, at
<figure><img src="image_77.png" /><figcaption>**Fig. 70.** *Left:* Energy densities, $\rho(z)$, of different components of the Universe as indicated in the legend. EDE (green) plays a subdominant rule at all redshifts, $z$. The fractional energy density $\Omega_{\text{EDE}}(z)$ (bottom subplot) peaks at the critical redshift, $z_c$, and decays around recombination, $z^*$. *Right:* A comparison of constraints on $H_0$ within the axion-like EDE model under different data sets as indicated in the figure, where solid (dashed) lines represent Bayesian (frequentist) constraints (see text).
*Source:* The constraints are taken from Refs. [1968,1970,2454,2457,2458].</figcaption></figure>

the critical redshift, $z_c$, which is determined by the onset of oscillations, along with the initial value of the scalar field $\theta_i$.

Predictions for EDE cosmologies have been found using the effective fluid theory approach [2454] and by directly solving the linearized scalar field evolution equations [2455,2459]. Extensions of CLASS [1949] using the latter approach are publicly available.⁴⁸⁴⁹ CAMB [1948] includes EDE using the effective fluid approximation.<fnref n="50" />

The axion-like EDE model faces some theoretical challenges: the desired degeneracy between $f_{\text{EDE}}$ and $H_0$, and with that, the ability of EDE to solve the tension sensitively depends on the choice of the parameters $z_c$ and $\theta_i$. Since there is no a priori reason for these parameters to have specific values, this leads to a “fine-tuning” problem. Moreover, the shape of the potential with index $n=3$ is considered fine-tuned itself since it requires that lower-order terms in the axion-like potential are zero. Therefore, the axion-like EDE model is often seen as a phenomenological model of EDE. Nonetheless, there are approaches to derive the axion-like EDE model from fundamental theories such as string theory [2460–2462] or to motivate $n=3$ by higher-order instanton corrections [2463].

#### 4.1.1.3. EDE variants.
From the axion-like EDE model, we learned valuable lessons about the challenges of EDE as a solution to the Hubble tension. In order to have the maximal possible leverage to resolve the tension, we need a specific form of potential, EDE needs to become relevant around matter–radiation equality and it needs to decay sufficiently fast before recombination. These theoretical challenges, or “fine-tunings”, have inspired many variants of EDE, which will be presented here.

Several EDE models consist of **slow-rolling scalar fields** minimally coupled to the metric, with different shapes of potentials [2447,2464–2466]. Rock ‘n’ Roll EDE assumes a power law potential of the form $V(\phi) \sim \phi^{2n}$, which approximates the axion-like potential for small $\phi/f$ [2467]. Acoustic DE generalizes the EDE model by using a phenomenological fluid description of EDE parameterized by the sound speed $c_s^2$ and equation-of-state parameter $w_{\text{ADE}}$ [2468].

A class of axion-like EDE models arises in the **string-inspired Chern–Simons (CS) gravity** [2469], which characterizes the low energy limit of string theory [2470]. In the simplest of such models [2471–2474], it leads to the generic form of Running Vacuum

Model (RVM) inflation [2475–2478]. In the specific context of string-inspired CS gravity, it is called Stringy RVM (StRVM) [2473,2479], in which there is a pre-RVM-inflationary era dominated by a stiff-fluid of gravitational (Kalb-Ramond) axions that characterize the massless gravitational string multiplet. At the end of such eras, condensation of primordial chiral GWs leads to a condensation of the gravitational anomaly CS term that couples the model to string axions that are either of Kalb-Ramond type or arise from compactification [2480]. The condensed CS term leads to linear terms in the axion potential [2481, 2482], responsible for RVM inflation, while the nonperturbative instanton effects lead to additional periodic modulations. In the context of the StRVM [2471–2473], the coefficient of the linear axion term is not a constant, but an RVM-type functional of even powers of the Hubble parameter. This differentiates the model from other linear axion monodromy models that arise from brane compactification [2483]. In Refs. [2481,2482] it was shown that the CS condensate leads to a *metastable* EDE, as there are imaginary parts in the effective action, pointing to a finite-lifetime of the respective vacua. Notably, different forms of the StRVM can also help alleviate, in modern eras, the Hubble and growth-of-structure tensions [2484,2485]. In the framework of $\alpha$-attractors in inflation [2486–2488], a potential for EDE was proposed in Ref. [2489] that can reproduce various cases based on its functional form.

**Unified models of EDE and late DE** have been proposed in Refs. [2490–2492]. Ref. [2490] introduces a model where a scalar field can explain both EDE and late DE in a joined manner in the context of $\alpha$-attractors. In Ref. [2491], a unified model of EDE and late DE is introduced, which can also enhance the primordial GW spectrum at PTA scales. Moreover, the EDE scalar potential can arise from a D-brane moving in the compact dimensions of warped compactifications in the context of D-brane models of cosmology. Ref. [2492] used a dynamical system approach to analyze a unified model of EDE and late DE. Scaling EDE is characterized by a constant EDE fraction in both the radiation- and matter-dominated epochs. This behavior is naturally found in quintessence models with exponential potentials and dates back to the concept of scalar field DE [296]; see also Ref. [2493,2494]. The very tight constraints imposed by the CMB disfavor this type of models as a solution to the tension [2448,2450,2495,2496], even when allowing for a coupling between EDE and DM [2428]. In the case of tracker solutions, [2497] studied EDE using quintessence and K-essence in the tracker regimes.

The possibility of EDE field coupling to other particles was also explored. Refs. [2459,2498,2499] introduce a model with a coupling

<footnote n="48">⁴⁸ https://github.com/PoulinV/AxiCLASS</footnote>

<footnote n="49">⁴⁹ https://github.com/mwt5345/class_ede</footnote>

<footnote n="50">⁵⁰ https://github.com/cmbant/CAMB/tree/master</footnote>