EuroSciPy 2026

Summary:
We present a likelihood-free inference framework, built entirely in Python (PyTorch, Optuna, NumPy/SciPy), that combines Gaussian process emulation with neural networks to jointly infer 28 cosmological and astrophysical parameters from galaxy cluster observations. For the first time, we achieve correlation coefficients above 0.97 for all cosmological parameters and above 0.90 for all astrophysical ones. The work has been published in The Astrophysical Journal (2025).

Detailed Description

1. Motivation and Scientific Context

Galaxy clusters are the most massive objects in the universe held together by gravity. Their properties are determined by two fundamentally different influences: cosmological parameters governing the large-scale evolution of the universe (such as the total matter density, the expansion rate, and the amplitude of matter fluctuations), and astrophysical processes occurring within them (such as star formation, supernova explosions, and energy injection from supermassive black holes).

In observational data, these two sets of effects are deeply entangled. Separating them is critical for two reasons. First, unresolved astrophysical uncertainties are a leading source of systematic error in cosmological measurements, contributing to ongoing discrepancies between different experiments measuring quantities like the expansion rate of the universe. Second, understanding feedback processes in clusters is essential for calibrating future cosmological simulation codes.

This challenge is becoming urgent. Telescopes such as eROSITA (X-ray), the Simons Observatory and CMB-S4 (microwave), and Euclid and the Vera Rubin Observatory (optical/infrared) are now delivering multiwavelength cluster observations at unprecedented sensitivity. The complexity of this data demands inference methods that can handle high-dimensional parameter spaces without restrictive assumptions about the form of the likelihood function.

2. Simulation Suite and Gaussian Process Emulation

Our work is built on the CAMELS-zoomGZ simulation suite: 768 hydrodynamic zoom-in simulations of galaxy clusters using the IllustrisTNG galaxy formation model, which is governed by 28 free parameters (5 cosmological, 23 astrophysical controlling star formation, winds, black hole accretion, and AGN feedback). The simulations sample this 28-dimensional space using a Sobol sequence for efficient coverage.

Since densely populating a 28-dimensional space with simulations is computationally prohibitive, we use the CARPoolGP emulator. This Gaussian process-based tool exploits correlations between pairs of simulations run with the same initial conditions but at different parameter space locations, combined with active learning to optimally place simulations. This allows CARPoolGP to produce low-variance emulations of averaged cluster profiles at any point in the 28-dimensional space.

From each cluster we extract five types of radial profiles: gas density, gas temperature, metallicity, X-ray surface brightness (0.5 to 2 keV), and the Compton-y parameter (thermal pressure along the line of sight). The emulator generates stacked versions of these profiles at arbitrary parameter space locations, producing the large training sets required for neural network training.

3. Neural Network Pipeline

The inference task is a regression problem: given a 1D vector of concatenated stacked profiles (29 to 148 values depending on which profile types are used), predict the posterior mean and standard deviation for each of the 28 parameters.

The architecture consists of fully connected blocks (linear layer, LeakyReLU activation, dropout), outputting 56 values: a mean and standard deviation per parameter. The loss function, following Jeffrey and Wandelt (2020), ensures outputs correspond to the first two posterior moments without assuming any posterior shape. This makes the approach entirely likelihood-free.

The pipeline is implemented in Python:

• PyTorch for model definition, training, and GPU-accelerated inference. The architecture is modular, with layers, neurons, and dropout as hyperparameters.
• Optuna for Bayesian hyperparameter optimization via Tree-Structured Parzen Estimation, running 1,000+ trials per experiment optimizing architecture, learning rate, weight decay, and dropout.
• NumPy/SciPy for all data preprocessing: profile extraction, concatenation, normalization, noise injection, radial cuts, and Sobol sequence sampling.
• Matplotlib for all visualizations including radar charts, bar charts, and predicted-vs-true plots.

Training uses Adam with batch size 256 for 1,000 epochs with early stopping. Data is split 70/15/15 into training, validation, and test sets. We found that 30,000 parameter space locations with five profile types reach near-saturation in accuracy, with marginal gains beyond that.

4. Key Results

With all five profile types concatenated, the network achieves correlation coefficients of 0.99 for matter density, baryon density, and Hubble constant, 0.97 for the spectral index and amplitude of matter fluctuations, and above 0.90 for all 23 astrophysical parameters. To our knowledge, this is the first time such accuracy has been achieved across the full 28-dimensional parameter space using cluster observables.

We attribute this to the fact that each parameter introduces distinct, non-degenerate features in the radial profiles. By varying one parameter at a time, we confirmed that no two parameters produce the same signature across all profile types, which is what allows the network to disentangle 28 simultaneous effects.

5. Robustness Analyses

We conducted four systematic studies to stress-test the results:

Noise sensitivity: Gaussian noise at 10% to 40% of the bin signal (signal-to-noise ratios from 10 to 2.5). Key parameters such as matter density, Hubble constant, and the IMF slope remain robust even at 40% noise. At a signal-to-noise of 10, typical for current X-ray observations, all parameters retain correlation coefficients above 0.7.

Radial cuts: Progressively truncating profiles from the outside reveals that the innermost 10% of the virial radius contains the bulk of the constraining information, with only modest performance loss.

Mass dependence: Extending to clusters from 10^13 to 3 x 10^14 solar masses, cosmological parameters remain well constrained across all masses. Astrophysical parameters show a mild decline at higher masses, consistent with more massive clusters approaching self-similar behavior.

Profiles vs. integrated quantities: Full radial profiles consistently outperform single integrated values, with the largest gap for temperature profiles. This demonstrates that spatial information within clusters carries significant constraining power lost when collapsing to a single number.

7. Relevance to the EuroSciPy Community

This work demonstrates a complete scientific Python pipeline from simulation-based data generation through Gaussian process emulation to neural network inference and visualization. The methodological core, combining GP emulation for scalable training data with likelihood-free neural network inference over a high-dimensional space, is transferable to any domain where forward simulations are expensive, parameter spaces are large, and likelihoods are intractable. Examples include climate modeling, fluid dynamics, and materials science. The talk will emphasize this transferable methodology alongside the scientific results.