# Towards Robustness Across Cosmological Simulation Models TNG, SIMBA, ASTRID, and EAGLE

Yongseok Jo, Shy Genel, Anirvan Sengupta, Benjamin Wandelt, Rachel

Somerville, Francisco Villaescusa-Navarro

► To cite this version:

Yongseok Jo, Shy Genel, Anirvan Sengupta, Benjamin Wandelt, Rachel Somerville, et al. Towards Robustness Across Cosmological Simulation Models TNG, SIMBA, ASTRID, and EAGLE. The Astrophysical Journal, 2025, 991 (1), pp.120. 10.3847/1538-4357/adec78. hal-04980315

HAL Id: hal-04980315

https://hal.science/hal-04980315v1

Submitted on 31 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

# Toward Robustness across Cosmological Simulation Models ILLUSTRISTNG, SIMBA, ASTRID, and SWIFT-EAGLE

Yongseok Jo<fnref n="1" />,<fnref n="2" /> ID, Shy Genel<fnref n="1" />,<fnref n="2" /> ID, Anirvan Sengupta<fnref n="3" />,<fnref n="4" />,<fnref n="5" /> ID, Benjamin Wandelt<fnref n="1" />,<fnref n="6" /> ID, Rachel Somerville<fnref n="1" /> ID, and
Francisco Villaescusa-Navarro<fnref n="1" />,<fnref n="7" /> ID

Received 2024 October 30; revised 2025 July 3; accepted 2025 July 4; published 2025 September 19

## Abstract

The rapid advancement of large-scale cosmological simulations has opened new avenues for cosmological and astrophysical research. However, the increasing diversity among cosmological simulation models presents a challenge to the robustness. In this work, we develop the Model-Insensitive ESTimator (MIEST), a machine that can robustly estimate the cosmological parameters, Ω_m and σ_8, from neural hydrogen maps of simulation models in the Cosmology and Astrophysics with MachinE Learning Simulations project—ILLUSTRISTNG, SIMBA, ASTRID, and SWIFT-EAGLE. An estimator is considered robust if it possesses a consistent predictive power across all simulations, including those used during the training phase. We train our machine using multiple simulation models and ensure that it only extracts common features between the models while disregarding the model-specific features. This allows us to develop a novel model that is capable of accurately estimating parameters across a range of simulation models, without being biased toward any particular model. Upon the investigation of the latent space—a set of summary statistics, we find that the implementation of robustness leads to the blending of latent variables across different models, demonstrating the removal of model-specific features. In comparison to a standard machine lacking robustness, the average performance of MIEST on the unseen simulations during the training phase has been improved by ~17% for Ω_m and 38% for σ_8. By using a machine learning approach that can extract robust, yet physical features, we hope to improve our understanding of galaxy formation and evolution in a (subgrid) model-insensitive manner, and ultimately, gain insight into the underlying physical processes responsible for robustness.

Unified Astronomy Thesaurus concepts: Cosmological parameters from large-scale structure (340); Cosmological parameters (339); Hydrodynamical simulations (767); Convolutional neural networks (1938); Robust regression (1949); Galaxy evolution (594); Cosmological evolution (336)

## 1. Introduction

The field of cosmological simulations has undergone significant advances in recent years, driven by both the rapid development of modern technology and the increasing sophistication of computational methods (A. Pillepich et al. 2018; R. Davé et al. 2019). As the number of cosmological simulations continues to grow, a major challenge has emerged: the problem of calibration (Y. Jo et al. 2023). Furthermore, with a variety of models available, there is currently no clear consensus on which simulation method or result is the most accurate or physical. This issue fundamentally arises from the significant divergence in simulation outcomes, where different subgrid models produce markedly disparate predictions, particularly for parameters on which they have not been calibrated. This uncertainty underscores the need for more robust methods to validate and compare simulations.

Two notable projects illustrate different approaches to this issue: the AGORA project (J.-h. Kim et al. 2014) and the

Cosmology and Astrophysics with MachinE Learning Simulations (CAMELS) project (F. Villaescusa-Navarro et al. 2021a) that is closely associated with the Learning the Universe (LtU) collaboration.<fnref n="8" /> The AGORA Project endeavors to achieve convergence across different simulation codes by implementing standardized conditions and rigorously comparing outcomes. In contrast, the CAMELS project acknowledges and leverages the intrinsic variability among simulations, concentrating on the simulation of a diverse array of cosmological models as they intrinsically exist.

However, in the context of the AGORA project, attaining convergence proves to be a complex endeavor, frequently necessitating meticulous calibration and the acceptance of discrepancies across various simulations. Furthermore, despite its reputation as an exemplary platform for the cross-validation of diverse simulations, the resultant converged simulations do not inherently ensure physical robustness. On the other hand, despite the extensive coverage offered by the suite of CAMELS simulations within the cosmological model space, it remains uncertain whether this space comprehensively encompasses the true structure of the Universe. These

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA; yj2812@columbia.edu</footnote>

<footnote n="2">²Columbia Astrophysics Laboratory, Columbia University, 550 West 120th Street, New York, NY 10027, USA</footnote>

<footnote n="3">³Department of Physics and Astronomy, Rutgers University, Piscataway, NJ 08854, USA</footnote>

<footnote n="4">⁴Center for Computational Quantum Physics, Flatiron Institute, New York, NY 10010, USA</footnote>

<footnote n="5">⁵Center for Computational Mathematics, Flatiron Institute, New York, NY 10010, USA</footnote>

<footnote n="6">⁶Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis boulevard Arago, 75014 Paris, France</footnote>

<footnote n="7">⁷Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544, USA</footnote>

<footnote n="8">⁸ https://learning-the-universe.org</footnote>
competing methodologies give rise to a pivotal question: “Should we prioritize the execution of a broader array of diverse simulations, or should our efforts be concentrated on the development of a convergent model capable of replicating a maximal number of observations?”

Concurrently, machine learning techniques are progressively being incorporated into cosmological simulations through various approaches, including the acceleration and enhancement of the simulation process (S. He et al. 2019; Y. Jo & J.-h. Kim 2019; D. Jamieson et al. 2023; T. Nguyen et al. 2024), as well as the estimation of cosmological and astrophysical parameters (N. S. M. de Santi et al. 2023b; C. Hahn et al. 2023; H. Shao et al. 2023; M. E. Lee et al. 2024). These endeavors have demonstrated significant progress and have become an essential component of cosmological simulations. Nevertheless, investigations that focus on the variety of cosmological simulations and their corresponding model spaces are relatively absent. In particular, the problem of robustness has emerged within the context of estimations. Robustness here refers to a phenomenon in which an estimator exhibits consistent performance across different models including those that are not included during the training phase. Despite the efforts with regards to the robustness of cosmological simulations (H. Shao et al. 2022; P. Villanueva-Domingo & F. Villaescusa-Navarro 2022; N. S. M. de Santi et al. 2023; N. Echeverri-Rojas et al. 2023), these attempts are limited by the dependence on manual identification of reliable physical quantities for parameter estimation.

In the meantime, the machine learning community is endeavoring to establish a rigorous methodology for achieving robustness, often referred to as domain generalization (K. Muandet et al. 2013; C. Finn & S. Levine 2018; Y. Li et al. 2018; H. Zhao et al. 2019; Z. Huang et al. 2020; G. Blanchard et al. 2021; J. Wang et al. 2023; Y. Zhang et al. 2023; K. Zhou et al. 2023). Domain generalization aims to train models that can generalize effectively to new, unseen domains without requiring access to data from those domains during training. In contrast to domain adaptation, which permits models to access target domain data (either labeled or unlabeled) to adjust their learning, domain generalization prohibits any such access. This constitutes a formidable yet highly pragmatic challenge, particularly in real-world scenarios where new domains may emerge continuously. Domain generalization usually includes some pivotal techniques such as domain shift, invariant feature learning, meta-learning, regularization, and ensemble learning.

The recent work in astrophysics has adopted domain shift to infer cosmological parameters from observational data, utilizing a machine trained on multiple simulations (J.-Y. Lee et al. 2024). The domain shift adopted in J.-Y. Lee et al. (2024) works in a way that introduces certain shifts in parameter space so that the parameters can be tuned to the true values. However, a significant limitation in this tuning process is the reliance on Sloan Digital Sky Survey–calibrated simulations for model calibration. Similarly, S. Andrianomena & S. Hassan (2025) presented a method to address covariate shift between data sets of the same observable with different distributions, improving the generalizability of neural networks for cosmological inference. Using H1 maps from IllustrisTNG and SIMBA, the authors employ adversarial training and optimal transport to adapt a pre-trained neural

network for out-of-distribution data without requiring labels. However, this does not constitute true generalization in principle, as the model has indirectly learned about the target data set through domain adaptation.

As part of the LtU collaboration that aims to reconstruct cosmological parameters, initial conditions, and astrophysical processes from upcoming observational data by leveraging machine learning to accelerate forward modeling, our primary objective is to address the robustness of multiple cosmological simulations using domain generalization in an adversarial fashion. In the line of efforts of having a diverse Universe, the suite of the CAMELS simulations now possesses several different cosmological simulation models, including ILLUSTRISTNG, SIMBA, ASTRID, and SWIFT-EAGLE that have reached $z=0$ with 1000 variations of cosmological and astrophysical parameters. Previous research has demonstrated that the different cosmological models exhibit unique signatures, especially in temperature maps and gas density maps (F. Villaescusa-Navarro et al. 2021a, 2022; V. Busillo et al. 2023; M. T. Tillman et al. 2023; I. Medlock et al. 2024), which subsequently influence the cross-estimation of cosmological parameters with a significant bias across different cosmological models; i.e., they are not robust (e.g., F. Villaescusa-Navarro et al. 2021b, 2021c; Y. Ni et al. 2023). Among these maps, this work focuses on neutral hydrogen (H1) maps. Neutral hydrogen is a key tracer of dark matter and large-scale structure in the Universe, making it a primary focus for several ongoing and future cosmological surveys, including major projects like the Square Kilometer Array, Hydrogen Epoch of Reionization Array, Low Frequency Array, the Vera C. Rubin Observatory’s Legacy Survey of Space and Time, the Nancy Grace Roman Space Telescope, SPHEREx, and Euclid. Moreover, H1 maps have been employed as tracers of cosmology within the realm of machine learning applications (e.g., S. Hassan et al. 2022).

With H1 maps, we build a Model-Insensitive ESTimator (MIEST) that estimates cosmological parameters in a robust manner that consistently performs on different cosmological simulation models. MIEST processes H1 maps using a convolutional neural network (CNN) and compress them into the latent variables, inspired by the Information Bottleneck (IB) method (N. Tishby et al. 1999). We seek to have these latent variables have limited information on the type of simulation that generated the original data. While A. A. Alemi et al. (2017) has provided a variation approach using a neural network to IB, we choose a slightly different objective function, which utilizes de-classification. In our approach, the latent variables are trained to get rid of classifiable information that is intrinsic to each simulation model, achieving robustness in MIEST. In this end, we aim at providing new insights into the configuration of simulations within the cosmological model space, thereby paving the way for more reliable comparisons and a clearer understanding of simulation results. Ultimately, we anticipate that our findings could suggest new directions for future cosmological simulations, helping to shape the field’s evolution toward greater accuracy and consistency.

The structure of this paper is as follows: In Section 2.1, we introduce the CAMELS simulations and the descriptions of ILLUSTRISTNG, SIMBA, ASTRID, and SWIFT-EAGLE that we adopt. In Section 2.2, we describe the structure and training procedure of MIEST in detail. Section 2.3 summarizes the
dimensionality reduction tools adopted to visualize the latent
space. In Section 3, we investigate the properties of the latent
space in response to implementation of robustness. In
Section 4, we provide the quantitative analysis on the
improved performance of the MIEST with robustness. In
Section 5, we discuss the definition of robustness, stability, and
convergence of MIEST, and comparisons between any two
simulations. Section 6 provides the summary of the MIEST.

**2. Methodology**

## 2.1. Cosmological Simulations

CAMELS<fnref n="9" /> is a suite of thousands of cosmological
simulations run with various cosmological models such as
ILLUSTRISTNG, SIMBA, ASTRID, and SWIFT-EAGLE (Y. Li
et al. 2018; F. Villaescusa-Navarro et al. 2021a). Each
simulation contains 256³ dark matter particles of mass
6.49 × 10⁷ (Ω_m − Ω_b)/0.251 h⁻¹ M☉ and 256³ gas cells with
an initial mass of 1.27 × 10⁷ h⁻¹ M☉ in a periodic box of a
comoving volume of (25 h⁻¹ Mpc)³. The cosmological and
astrophysical parameters of interest are Ω_m, σ_8, A SN1, A SN2,
A AGN1, and A AGN2 (refer to F. Villaescusa-Navarro et al. 2021a
for details). In this study, we adopt the LH suite for the training
data set. The 1000 simulations of the LH set are run with
Ω_m ∈ [0.1, 0.5], σ_8 ∈ [0.6, 1.0], A SN1 ∈ [0.25, 4.0],
A SN2 ∈ [0.5, 2.0], A AGN1 ∈ [0.25, 4.0], and
A AGN2 ∈ [0.5, 2.0] arranged in a latin hypercube. The
cosmological parameters are fixed in all simulations as
Ω_b = 0.049, h = 0.6711, n_s = 0.9624, M_ν = 0.0 eV, w = −1,
and Ω_K = 0.

2.1.1. *IllustridTNG*

The IllustrisTNG model utilizes the AREPO code (V. Springel 2010; R. Weinberger et al. 2020) to solve the coupled equations of gravity using TreePM and magnetohydrodynamics using a Voronoi moving mesh approach, and includes a comprehensive set of baryonic physics modules (R. Weinberger et al. 2017; A. Pillepich et al. 2018, for a full description). The subgrid models cover radiative cooling and heating (N. Katz et al. 1996; R. P. C. Wiersma et al. 2009), hydrogen self-shielding (A. Rahmati et al. 2013), star formation, and stellar evolution (V. Springel & L. Hernquist 2003; M. Vogelsberger et al. 2013).

In addition, galactic winds from stellar feedback are
implemented kinetically through temporarily hydrodynamami-
cally decoupled particles, ejected stochastically and isotropi-
cally from star-forming gas, based on V. Springel &
L. Hernquist (2003) with 10% thermal energy addition. Wind
speed and energy injection rate per unit star formation depend
on gas conditions (metallicity and dark matter velocity
dispersion), redshift, and two global normalization parameters.
Here, the energy injection rate and the wind speed are
modulated by the parameters A<sub>SN1</sub> and A<sub>SN2</sub>, respectively.

Supermassive black hole (SMBH) particles with initial mass $M_{\text{seed}} = 8 \times 10^5 M_{\odot}/h$ are seeded in halos with $M_{\text{FoF}} > 5 \times 10^{10} M_{\odot}/h$ (see R. Weinberger et al. 2017, for details of the SMBH model). SMBH particles are repositioned at the potential minimum and merge if within each other's feedback spheres. Gas accretion follows the spherical H. Bondi (1952) with the Eddington cap. SMBH feedback employs

thermal, kinetic, and radiative modes. Thermal and kinetic
modes operate separately, distinguished by the Eddington ratio
of the SMBH's accretion rate. The transition from high to low
mode generally occurs around $M_{\text{SMBH}} \sim 10^8 M_{\odot}$. The high
accretion rate thermal mode injects thermal energy into a
feedback sphere with an energy conversion efficiency of 0.02.
The astrophysical parameters $A_{\text{AGN1}}$ and $A_{\text{AGN2}}$ govern energy
injection in the kinetic feedback mode, dependent upon the
accretion rate and the total energy allocation in the low
accretion mode.

2.1.2. SIMBA

The SIMBA model utilizes the *N*-body/hydrodynamics
code GIZMO (P. F. Hopkins 2015) with the “Meshless Finite
Mass” hydrodynamics solver (R. Davé et al. 2019, for a full
description). Gravitational forces are computed using the
modified TreePM algorithm of the GADGET-III code
(V. Springel 2005), including adaptive gravitational softenings
for the gas, stellar, and dark matter components. The subgrid
models include radiative cooling and photoionization
(F. Haardt & P. Madau 2012; B. D. Smith et al. 2017),
hydrogen self-shielding (A. Rahmati et al. 2013), star
formation, and feedback and evolution (M. R. Krumholz &
N. Y. Gnedin 2011; R. Davé et al. 2016).

Stellar feedback drives galactic winds that are kinetically
implemented through hydrodynamically decoupled, two-
phase, metal-enriched outflows. In these outflows, 30% of
the wind particles are heated to a temperature determined by
subtracting the wind kinetic energy from the SNe energy. Star-
forming gas elements are ejected stochastically based on
predetermined values of the mass loading factor and wind
speed. The mass loading factor is derived from the FIRE
“zoom-in” simulations (P. F. Hopkins et al. 2014) and adjusts
according to the galaxy’s stellar mass, in line with D. Anglé-
s-Alcázar et al. (2017b). In CAMELS, A<sub>SN1</sub> is introduced to
regulate the overall normalization of the mass loading factor,
while A<sub>SN2</sub> adjusts the wind speed normalization for varying
parameters.

The black hole model in SIMBA includes gravitational
torque accretion and kinetic feedback (D. Anglés-Alcázar
et al. 2017a). Galaxies with stellar mass $M_* > 10^{9.5} M_\odot$ are
seeded with black holes of initial mass $M_{\text{seed}} = 10^4 M_\odot/h$
using the on-the-fly Friends-of-Friends (FoF) method, based
on the stellar mass threshold for seeding. Black hole particles
are relocated to the potential minimum's position within the
FoF host group if it lies within a distance of $<4 \times R_0$, where
$R_0$ represents the size of the black hole kernel that encloses the
nearest 256 gas elements. Black holes within $R_0$ of each other
merge instantaneously if their relative velocity is less than 3
times their mutual escape velocity. The accretion processes
governing black hole growth encompass the accretion of cold
gas as delineated by P. F. Hopkins & E. Quataert (2011), as
well as the accretion of hot gas described by H. Bondi (1952).
Black hole feedback consists of high mass loading outflows in
the radiative “QSO” mode and lower mass loading but faster
outflows at low Eddington ratios in the jet mode, both of which
eject gas elements bipolarly. Black holes with $M_{\text{BH}} > 10^{7.5} M_\odot$
that accrete below an Eddington ratio of 0.2 are switched to the
jet mode. $A_{\text{AGN1}}$ modulates the total momentum flux of both
modes, while $A_{\text{AGN2}}$ controls the maximum jet speed for the jet
mode (R. Davé et al. 2019, for a full review of the black hole
model).

<footnote n="9">⁹ https://www.camel-simulations.org</footnote>
### 2.1.3. ASTRID

The ASTRID model utilizes a new version of the MP-Gadget simulation code, a massively scalable version of the cosmological structure formation code Gadget-3 (V. Springel 2005), to solve gravity with TreePM and hydrodynamics with a smoothed particle hydrodynamics (SPH) method (full review S. Bird et al. 2022 and Y. Ni et al. 2023). The subgrid model includes radiative cooling and photoionization heating (N. Katz et al. 1996), metal line cooling with tracing the gas and stellar metallicities (M. Vogelsberger et al. 2014), a spatially uniform ionizing background (C.-A. Faucher-Giguère 2020), hydrogen self-shielding (A. Rahmati et al. 2013), and star formation (V. Springel & L. Hernquist 2003).

Winds driven by stellar feedback within galaxies are kinetically implemented through particles that are temporarily decoupled from hydrodynamics. These winds originate from newly formed star particles, which randomly select gas particles within their SPH smoothing length to transform into wind particles. A particle becomes recoupled after its density decreases by a factor of 10. Wind particles do not engage in or generate pressure but do receive mass return, undergo cooling, and participate in density estimations. In the CAMELS suite, $A_{SN1}$ controls the total energy injection rate per star formation, while $A_{SN2}$ regulates the speed of the SN wind.

The SMBH model adopted the BH seeding and dynamics prescriptions applied in BLUETIDES (Y. Feng et al. 2016). BHs with initial mass $M_{seed} = 5 \times 10^5 M_{\odot}/h$ are seeded in halos with $5 \times 10^{10} M_{\odot}/h$ by means of the on-the-fly FoF halo finding. BH particles are repositioned to the location of the local potential minimum at each active time step, and two BHs located within twice the gravitational softening length of each other are instantaneously merged (Y. Ni et al. 2022). BH accretion rates follow the Bondi–Hoyle–Lyttleton–like prescription (T. Di Matteo et al. 2005). Super-Eddington accretion periods are allowed, but are capped at twice the Eddington rate. BH bolometric luminosity is proportional to the accretion rate with a 0.1 mass-to-light conversion efficiency (N. I. Shakura & R. A. Sunyaev 1973). SMBH feedback uses a two-mode approach: thermal and kinetic, based on the Eddington ratio of the current BH accretion rate. The kinetic mode is activated contingent upon an Eddington threshold for black holes with $M_{BH} = 5 \times 10^8 M_{\odot}/h$ with an Eddington threshold cap of 0.05. $A_{AGN1}$ and $A_{AGN2}$ modulate the kinetic and thermal feedback efficiency, respectively, with $A_{AGN1}$ retaining the same physical meaning as in the TNG suite.

### 2.1.4. Swift-EAGLE

The SWIFT-EAGLE model utilizes the SPH and gravity code Swift (M. Schaller et al. 2018, 2024; J. Borrow et al. 2022). Swift is a parallel, open-source, versatile and modular code, with a range of hydrodynamics solvers, gravity solvers, and subgrid models for galaxy formation. In the CAMELS suite, the SPHENIX flavor of SPH is adopted, coupled with a modified version of the Evolution and Assembly of GaLaxies and their Environments (EAGLE) subgrid model for galaxy formation and evolution (see J. Schaye et al. 2015 and R. A. Crain et al. 2015). This includes element-by-element radiative cooling and heating rates from S. Ploeckinger & J. Schaye (2020), star formation (J. Schaye & C. Dalla Vecchia 2008), stellar evolution and enrichment from R. P. C. Wiersma et al. (2009), and single thermal-mode

feedback from massive stars and accreting active galactic nuclei (AGNs; see C. M. Booth & J. Schaye 2009; C. Dalla Vecchia & J. Schaye 2012; Y. M. Rosas-Guevara et al. 2015).

Given that SWIFT-EAGLE is a recent addition to the CAMELS suite and the accompanying release paper is currently in preparation (C. Lovell et al. 2025, in preparation), this paper provides a detailed description of each parameter for the benefit of the readers. $A_{SN1}$ represents the thermal energy injected in each SNII feedback event:

$$E_{\text{SNII}} = A_{\text{SN1}} \times 10^{51} \text{ erg}, \quad (1)$$

which falls within the range of $[2.5 \times 10^{50}, 4 \times 10^{51}]$. $A_{SN2}$ represents the metallicity dependence of the feedback fraction per unit stellar mass:

$$f_{\text{th}} = f_{\text{th,min}} + \frac{f_{\text{th,max}} - f_{\text{th,max}}}{1 + \left(\frac{z}{0.1Z_{\odot}}\right)^{n_Z} \left(\frac{n_{\text{H,birth}}}{n_{\text{H,0}}}\right)^{-n_n}}, \quad (2)$$

where

$$n_Z = A_{\text{SN2}}^2. \quad (3)$$

$A_{AGN1}$ represents the scaling of the black hole Bondi accretion rate:

$$\dot{m}_{\text{accr}} = \alpha \frac{4\pi G^2 m_{\text{BH}}^2 \rho}{(c_s^2 + v^2)^{3/2}}, \quad (4)$$

where

$$\alpha = A_{\text{AGN1}}. \quad (5)$$

$A_{AGN2}$ represents the temperature jump of gas particles in AGN feedback events.

$$\Delta T_{\text{AGN}} = A_{\text{AGN2}}^2 \times 10^{8.28}, \quad (6)$$

which ranges from $4.76 \times 10^7$ to $7.62 \times 10^8$.

### 2.1.5. H I Maps

In this study, we employ projected 2D neutral hydrogen (H I) density maps as input to allow the machine to estimate the cosmological parameters, utilizing the CAMELS Multi-field Dataset (F. Villaescusa-Navarro et al. 2022). The H I field represents the spatial distribution of the neutral hydrogen density. To create 2D maps, it is necessary to read the locations and neutral hydrogen masses of all gas particles within a given simulation. The neutral hydrogen mass of each gas particle is determined using the self-shielding fitting formula outlined in A. Rahmati et al. (2013), while star-forming particles are assumed to be completely self-shielded from external radiation, as described in F. Villaescusa-Navarro et al. (2018). The value stored in each pixel is the total neutral hydrogen density from all gas particles that contribute to that specific location, or more explicitly

$$\frac{1}{Q} \sum_{i} M_{\text{H},i} \int_{x} W_{g,i}(\mathbf{x} - \mathbf{r}_{g,i}) d\mathbf{x}, \quad (7)$$

where $W_{g,i}$ is the window function, $M_{\text{H},i}$ is the neutral hydrogen mass, and $\mathbf{r}_{g,i}$ is the position of gas particle i. The sum and integral are over all gas particles and the area of the pixel, respectively. Q denotes the pixel area in 2D maps.

Two-dimensional maps are constructed by initially extracting the spatial coordinates and attributes of a specified field
from the simulation data. Each simulation comprises 15 slices: five in each of the XY, XZ, and YZ planes, respectively. Particles located within a slice of the volume (25 × 25 × 5)(h⁻¹Mpc)³ are selected, such that slices oriented in the same projection direction are nonoverlapping. The particle radius R is determined by calculating the distance to the 32nd nearest gas particle, and the window function of the particles, W₂D(r, θ), is mapped onto a 2D plane, where

$$W_{3D}(r, \theta, z) = \begin{cases} \frac{3}{4\pi R^3} & \text{if } |r| < R, \\ 0 & \text{otherwise} \end{cases} \quad (8)$$

and W₂D(r, θ) = ∫z W₃D(r, θ, z)dz. Lastly, the attributes of the particles are superimposed on a grid of 256 × 256 using Equation (7). This procedure involves uniformly sampling 1000 tracers within each particle radius, subsequently assigning each tracer to the corresponding pixel. This numerical approximation method asymptotically approaches the precise outcome as the number of tracers increases. In summary, the maps are generated at a redshift of z = 0 and consist of 256 × 256 pixels spanning an area of 25 × 25(h⁻¹Mpc)². A total of 15,000 maps are produced for each simulation suite, with 15 maps per simulation in 1000 different simulations.

The correlation of the H I maps with cosmological parameters—Ω_m and σ_8—overall resembles the correlation between the large-scale structures of dark matter and cosmological parameters. For instance, larger Ω_m results in more pronounced structures as more matter is available for gravitational collapse, while a higher σ_8 leads to more tightly clustered, clumpy structures due to enhanced fluctuations. However, a notable difference from dark matter is that neural hydrogen is more sensitive to subgrid models or astrophysical phenomena such AGNs and supernovae. Subsequently, the disparity in subgrid models of cosmological simulations are imprinted on H I maps through the variation of the cosmological parameters. For example, in the case of AGN feedback, increasing the matter content in the Universe (high Ω_m) results in the formation of dense clusters and galaxies, which facilitates the growth of SMBHs and subsequently leads to a greater fraction of ionized diffuse gas in the Universe due to the jet feedback from these black holes. This can produce effects comparable to the increase in the AGN parameters. Consequently, the distinct subgrid models employed in cosmological simulations can significantly influence the notable disparity in H I maps across cosmological simulation models (refer to Section 5.3 for the detailed discussion).

## 2.2. Neural Network: MIEST<fnref n="10" />

The objective of this study is to eliminate information that distinguishes between simulation models, while simultaneously extracting information pertinent to the cosmological parameters. To achieve this, we utilize methods of compression and de-classification by combining the approaches of deep variational IB (A. A. Alemi et al. 2017) and generative adversarial network (I. Goodfellow et al. 2014). Here, de-classification refers to the procedure in which the machine imposes penalties on the classification mechanism to inhibit the identification of simulation models. Hereafter, our machine is referred as to MIEST.

Figure 1 illustrates the structure of the MIEST. The MIEST is composed of four parts:

1. The CNN $\mathcal{F}$ takes 2D maps $\mathbf{X}$ as input and manipulate them into 1D arrays;

2. The encoder $\mathcal{G}$ takes 1D arrays from the CNN and processes (usually compress) their information to latent variables $Z$;

3. The regressor $\mathcal{R}$ outputs the mean and standard deviation of cosmological parameters $\theta$ taking $Z$ as input;

4. The classifier $\mathcal{C}$ predicts labels of simulation models, $i$, and this information is used to train the CNN and the Encoder to omit the classifiable information.

Algorithm 1 elucidates the comprehensive training procedure. The training process comprises two distinct phases. First, we train the CNN, the encoder, and the regressor using the objective function $L_1$ (refer to Equation (9) in Section 2.2.1). Second, only the classifier is updated utilizing the cross-entropy loss $H$. The rationale for the separate training phases is motivated by our objective to leverage the classifier for the de-classification or removal of classifiable information within the CNN and predominantly the decoder, rather than optimizing for classifier performance. To this end, it is imperative that we possess a well-trained classifier readily available for the de-classification process. By repeatedly applying this procedure over multiple iterations with the same training data, the encoder is trained to generate more compact and compressed latent variables such that they are un-classifiable across different simulations, simultaneously enhancing the accuracy of the regressor.

### Algorithm 1. Model Insensitive Estimation

**Input:** Training set $\mathbf{X}$, iterations $N$

**Neural Network:** CNN $\mathcal{F}(\mathbf{X})$, Encoder $\mathcal{G}(\mathbf{Y})$, Classifier $\mathcal{C}(\mathbf{Z})$, $\mathcal{R}(\mathbf{Z})$

**Variable:** latent $Z$, simulation label $i$

**Output:** Mean of cosmological parameters $\bar{\theta}$, Standard deviation of cosmological parameters $\sigma(\theta)$

**for** $n = 1$ to $N$ **do**

**Phase I: Regression**

Process $\mathbf{X}$ through the CNN: $Y = \mathcal{F}(\mathbf{X})$

Compress $\mathbf{Y}$ through the encoder: $Z = \mathcal{G}(\mathbf{Y})$

Predict cosmological parameters $(\theta, \sigma) = \mathcal{R}(\mathbf{Z})$

Optimize $\mathcal{F}$, $\mathcal{G}$, and $\mathcal{R}$ using gradient decent with $\mathcal{L}_1 = J_1 + J_2 - \beta H + \gamma I$

**Phase II: Classifier For De-classification**

Process $\mathbf{X}$ through the CNN: $Y = \mathcal{F}(\mathbf{X})$

Compress $\mathbf{Y}$ through the encoder: $Z = \mathcal{G}(\mathbf{Y})$

Predict simulation labels $i = \mathcal{C}(\mathbf{Z})$

Optimize $\mathcal{C}$ using gradient decent with $\mathcal{L}_2 = H$

**end**

### 2.2.1. Equations

The objective function $\mathcal{L}_1$ in Algorithm 1 is written as

$$\mathcal{L}_1 = J_1 + J_2 - \beta H + \gamma I, \quad (9)$$

where

$$J_1 = \int (\theta - \mathcal{R}_{\theta}(z))^2 p(\theta, z) d\theta dz, \quad (10)$$

$$J_2 = \int ((\theta - \mathcal{R}_{\sigma}(z))^2 - \hat{\sigma}(z))^2 p(\theta, z) d\theta dz, \quad (11)$$

<footnote n="10">¹⁰ https://github.com/YongseokJo/MIEST.git</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 1. Structure of the neural network. The neural network comprises a convolution neural network (CNN) $\mathcal{F}$, an encoder $\mathcal{G}$, a classifier $\mathcal{C}$, and a regressor $\mathcal{R}$.</figcaption></figure>

$$H = -\int p(z, i) \log \mathcal{C}(z, i) dzdi, \qquad (12)$$

$$I = I(X; Z) = D_{KL}(P_{X,Z}||P_X \otimes P_Z). \qquad (13)$$

Here, the latent variables are given through the CNN and the encoder as $Z = \mathcal{G}(\mathcal{F}(X))$, and the regressor outputs two values —the mean and the standard deviation—as $\mathcal{R}(z) = (R_{\theta}(z), R_{\sigma}(z))$.

The brief descriptions of each loss term are as follows: $J_1$ and $J_2$ are the mean square errors of the mean and variance of the parameters, which are used for the estimation of the mean and variance of the parameters (N. Jeffrey & B. D. Wandelt 2020). The cross-entropy loss $H$ is used for the de-classification to discard the distinguishable information in the latent variables as in Figure 1. Lastly, the mutual information $I$ is implemented for compression of the input data into latent variables motivated by the deep variational IB (A. A. Alemi et al. 2017), where the IB method is a framework in information theory used to extract the most relevant information about a target variable from an input data while compressing the input data as much as possible.

### 2.2.2. Hyperparameter and Optimization

We use Optuna for tuning the hyperparameters within the neural network and the ones used in training such as the structure of the neural networks, dimension of latent variable, dropout, learning rate, and weight decay (T. Akiba et al. 2019). Optuna is an open-source hyperparameter optimization framework designed to enhance machine learning models by automating the search for optimal hyperparameters. It uses a define-by-run paradigm, allowing dynamic search space specification during execution, accommodating complex configurations. Optuna primarily employs Bayesian optimization, which we adopt in this work, but also offers grid

and random search strategies. Additionally, it uses pruning techniques to stop unpromising trials early, saving computational resources, and offers visualization tools for optimization analysis.

The objective in the optimization with Optuna is to minimize the regression error $J_1$ and minimize the normalized AUC score<fnref n="11" /> for the de-classification on the validation data. The normalized AUC score is defined as $|\text{AUCscore} - 0.5|$, of whose values range from 0–0.5, where 0 indicates nonclassifiability and 0.5 represents perfect classification performance. Thus, the optimization process is designed to enhance the estimation accuracy of cosmological parameters while inhibiting the classifier’s ability to differentiate among the simulation models.

The free parameters subject to optimization consist of the number of hidden layers, dropout, weight decay, learning rate, and the dimensionality of the latent space. Specifically, the dimensionality of the latent space, as defined by $\text{dim}(Z)$, is sampled from a uniform probability density ranging from 50–2000. This range has been chosen based on empirical evidence. In this study, the machines trained with and without robustness exhibit a difference in the dimensionality of the latent space, specifically 1375 and 793 dimensions, respectively. This variation in dimensionality may be associated with the quantity of physical information within the input-output pairs, depending on constraints such as the degree of robustness. However, we will not explore this relationship further in the work, as it lies beyond the scope of our current investigation.

In this work, we choose the one that is most accurate on the test set of the simulations, on which the machine is trained, as our best model. Presented in the result sections—Sections 3 and 4—are from the best model. Further investigation on variance will be discussed in Section 5.2

## 2.3. Dimensionality Reduction

Owing to the inherent high dimensionality of the latent space, discerning its implications presents a significant challenge. Accordingly, we employ dimensionality reduction techniques to facilitate a more detailed analysis and visualization of the latent space.

### 2.3.1. Uniform Manifold Approximation and Projection

Uniform Manifold Approximation and Projection (UMAP) is a nonlinear dimensionality reduction technique based on manifold learning and topological data analysis (L. McInnes et al. 2020). The goal of UMAP is to reduce the dimensionality of data while preserving as much of the local and global structure as possible. The algorithm first constructs a weighted k-nearest neighbor (k-NN) graph, where distances between points in the high-dimensional space are converted into a fuzzy topological representation using a Riemannian metric. Once this high-dimensional graph is constructed, UMAP proceeds to optimize it into a lower-dimensional space by minimizing a cross-entropy between the two topological representations, thereby ensuring that the relationships and distances from the

<footnote n="11">¹¹ The Area Under the Curve (AUC) score measures binary classification performance by representing the area under the ROC curve, which plots the true-positive rate against the false-positive rate at various thresholds. The AUC score ranges from 0–1, where 1 indicates perfect performance, 0.5 represents random guessing, and values below 0.5 suggest worse-than-random performance.</footnote>
original space are preserved in the low-dimensional embedding.

UMAP consists of several parameters, including “n neighbors,” which controls the size of the local neighborhood used for manifold approximation. Smaller values of “n neighbors” focus on preserving local structure, whereas larger values help capture more global structures. Another important parameter is “min dist,” which determines how closely points are packed together in the low-dimensional space, effectively balancing the preservation of local versus global structures. Additionally, UMAP supports various distance metrics, such as Euclidean and cosine, for computing distances between points in the input space. The selection of these parameters varies between different figures, in alignment with the specific objectives of our analysis. The parameters are mentioned in the accompanying text.

### 2.3.2. Kernel Principal Component Analysis

Kernel principal component analysis (KernelPCA) is an extension of traditional principal component analysis (PCA) that allows for the extraction of principal components in a higher-dimensional feature space (B. Schölkopf et al. 1997). KernelPCA employs the kernel trick to map the data into a higher-dimensional space using a kernel function, such as the radial basis function or polynomial kernel. Importantly, this mapping is done implicitly, with the kernel trick computing the dot product between data points in the transformed space without explicitly performing the transformation.

After the kernel function generates the kernel matrix, representing pairwise similarities between data points in the transformed space, the matrix is centered by subtracting the mean to ensure that the data has zero mean in the new feature space. PCA is then applied to this centered kernel matrix, where the eigenvectors corresponding to the largest eigenvalues are selected as the principal components. These components capture the most significant variance in the data after transformation through a kernel.

In general, UMAP is preferred for its scalability, speed, and ability to preserve local and global structure, making it ideal for visualization and exploratory analysis of large, high-dimensional data sets. On the other hand, kernelPCA excels in cases where nonlinear data structures are crucial, and specific kernel choices can exploit known relationships, but it is limited by scalability issues and requires careful tuning of kernel parameters. Hence, we utilize UMAP for broad applications due to its versatility, while implementing kernelPCA with the linear kernel specifically to interrogate the behavior of cosmological parameters within the latent space, where established relationships are presumed.

## 3. Latent Space

The aim of this work is to build a robust machine that performs consistently on different models. However, before diving directly into machine performance, we first investigate the anatomy of MIEST and then investigate how robustness in MIEST works. To this end, we investigate the configuration of the latent space that encapsulates the summary statistics derived from simulations.

The latent space refers to an abstract representation of data that can be obtained as output from the encoder that compresses the simulation maps through CNN (see

Section 2.2). By construction, MIEST is designed to process the information from the simulations into latent variables so that the latent variables are simultaneously expressive and robust. Hence, understanding the behaviors and characteristics of the latent space is crucial to analyze the mechanism or performance of MIEST and ultimately to grasp the intricate relationship between different simulations.

However, the dimension of the latent space ranges from $O(10^2)$ to $O(10^3)$. Substantial dimensionality poses considerable challenges for latent-space analysis. To address this, we employ dimensionality reduction techniques that yield insightful assessments of high-dimensional data (refer to Section 2.3). The reduced latent space is denoted as $\tilde{Z}_{UMAP}$ and $\tilde{Z}_{KPCA}$, respectively. Throughout this paper, we choose three dimensions to which the latent space is reduced. First of all, there are at least two axes reserved for $\Omega_m$ and $\sigma_8$; it is not necessary that the two $\Omega_m$ and $\sigma_8$ axes be orthogonal, but their combinations will be. Theoretically, there can be a third axis that represents simulation models or classifications of simulations. However, it is uncertain whether the classification axis should exist and whether the axis can be an orthogonal basis that does not overlap with other axes.

After examinations with multiple different dimensions for reduction with UMAP and KernelPCA, we have found that dimensions exceeding three cannot make any significant differences in visualization; the extra dimensions are redundant with respect to $\Omega_m$, $\sigma_8$, and classification. Additionally, no axis forms an orthogonal basis for classification, although it consists of a signature of classifications; presumably, it is overlapped with other information. Hence, we empirically select the three dimensions. Notice that these results are from UMAP and KernelPCA, implying that there may theoretically still exist an orthogonal basis for classification.

### 3.1. How MIEST Achieves Robustness

In this section, we use MIEST trained on ILLUSTRISTNG and SIMBA with and without robustness to investigate how robustness works within the framework of MIEST. Shown in Figure 2 are 2D projections of the 3D latent space $\tilde{Z}_{UMAP}$ of the whole data set<fnref n="12" /> that is dimension reduced through UMAP. For UMAP, we use the hyperparameters of 20 neighbors, 0.1 minimum distance, and the “Euclidean” metric.<fnref n="13" /> The upper panels exhibit latent variables $\tilde{Z}_{UMAP}$ of MIEST trained *without* robustness, while the lower panels show *with* robustness. The machines *with* and *without* robustness are trained, respectively, and end up having different dimensions of the latent space. Hence, we conduct the dimensionality reduction separately for each machine. In other words, there is no direct connection between the upper and lower panels, except that we use the same data set to train the machine.

The projections in the first column and second column are color-coded with the true values $\Omega_m$ and $\sigma_8$, respectively, where darker colors indicate higher values of the parameters. With this in mind, ILLUSTRISTNG (blue) and SIMBA (red) have similar trends in terms of cosmological parameters in the latent space along the axes, regardless of robustness. This suggests that the two axes—denoted as regression axes 0 and 1

<footnote n="12">¹² This includes both the training set and test set.</footnote>

<footnote n="13">¹³ From the experiments, the “canberra” metric is better at capturing clustering. However, the axes for the cosmological parameters become interwoven. For better visualization, we choose “Euclidean” for this overview plot.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Two-dimensional projections of the 3D-reduced latent variable $\tilde{Z}_{UMAP}$ using UMAP. Since UMAP naturally reduces the dimension of the data, the axes above are totally arbitrary given by UMAP. However, one can deduce that due to the nature of the data, there must be at least two axes that represent $\Omega_m$ and $\sigma_8$, respectively. We denote them as regression axes 0 and 1. Each data point is color-coded by their true parameter values $\Omega_m$ (first column) and $\sigma_8$ (second column). This shows that each regression axis is a combination of $\Omega_m$ and $\sigma_8$. Arguably, we assume that the third axis should be the classification axis that represents the simulation model of the data. Here, the reason they are sparse or porous is that we use a test set. Although the “No Robustness” case looks like ILLUSTRISTNG and SIMBA are blended together, it is the structure that two sheets of ILLUSTRISTNG and SIMBA are spanned really close to each other almost in parallel. One can confirm a clearer impact of robustness with the 3D construction of the latent variable space. In addition, it is manifested in the zoomed-in plot (see Figure 3). For 3D visualization using the entire data set, check out the movies of “No Robustness” (youtu.be/_6hqrcqTQ_c) and “With Robustness” (youtu.be/Xd5SIorPA1s).</figcaption></figure>

are correlated with $\Omega_m$ and $\sigma_8$, respectively. One must notice that the axes from the dimensionality reduction are completely arbitrary (that is, we have not provided any information on cosmological parameters). We empirically choose to use the three dimensions as the final reduced dimension. We anticipate that due to the intrinsic mutual information between the inputs and outputs, there must be two axes reserved for the cosmological parameters $\Omega_m$ and $\sigma_8$. Precisely, there should exist two axes, corresponding to $\Omega_m$ and $\sigma_8$, which can be constructed as a combination of the basis vectors within the reduced latent space $\tilde{Z}_{UMAP}$. For this reason, we denote two axes that are correlated with cosmological parameters as “regression axes 0 and 1” and the other axis as “axis 2.”

The structural characteristics of “No Robustness” are analogous to networks of blood vessels or large-scale structures, being composed of sheets and clumps. ILLUSTRISTNG and SIMBA represent two closely juxtaposed thin sheets. Significant separation is observed at elevated values of $\Omega_m$ and $\sigma_8$. Upon meticulous examination, slight misalignments in their clumpy structures become evident, indicating classificationability. We also investigate the classificationability of latent variables Z derived from two distinct simulations utilizing the Area Under the Curve (AUC) score. An AUC score of 0.97 is obtained, demonstrating that the latent variables Z of the two simulations are clearly classifiable.

Conversely, “With Robustness” exhibits smooth oval or spherical distributions that are consistent across both the ILLUSTRISTNG and SIMBA data sets, despite the UMAP

hyperparameters remaining identical to those in the “No Robustness,” which generally influence the structural features of the reduced space $\tilde{Z}_{UMAP}$. The distinction between the two simulations becomes visually indistinguishable in the reduced space $\tilde{Z}_{UMAP}$. Achieving an AUC score (see footnote 11) of 0.53, it successfully mitigates the classification of simulations, marking a progression toward robustness also in the full latent space $Z$. However, it should be noted that the correlation between cosmological parameters and the regression axes weakens compared to the “No Robustness” case overall. This potentially signifies a loss of cosmological information within the latent variables Z, as discussed in Section 4. An important caveat is that the observed clumpiness or any characteristic structures in the reduced latent space $\tilde{Z}_{UMAP}$ are likely contingent upon the choice of dimensionality reduction methodology, and hence are not physically meaningful.

We turn to the effect of robustness on the test set that the MIEST has never seen. Figure 3 shows the latent variables $\tilde{Z}_{UMAP}$ divided into three parts with respect to the values of $\Omega_m = [0.1, 0.25], [0.25, 0.35], [0.35, 0.5]$—denoted lower, intermediate, and upper, respectively. Each column shows colored data points in latent space $\tilde{Z}_{UMAP}$ for each $\Omega_m$ range. The latent space $\tilde{Z}_{UMAP}$ is again reduced by UMAP using hyperparameters of 20 neighbors, minimum distance of 0.5, and the “canberra” metric. The selection of the “canberra” metric accentuates the clustering of simulations more effectively. The axes presented have been rotated in an arbitrary
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Two-dimensional projections of the 3D-reduced latent variable space of the test set of ILLUSTRISTNG (blue) and SIMBA (red) via UMAP. The axes are rotated arbitrarily for better visualization. Each column shows latent variables from three ranges of Ω_m: lower, intermediate, and upper, with data points color-coded by their own ranges. The upper panels represent the “No Robustness” case, where ILLUSTRISTNG and SIMBA are classifiable. The lower panels show the “With Robustness” case, where classification is less evident.</figcaption></figure>

manner to improve the visualization of the clustering between ILLUSTRISTNG and SIMBA data sets.

"No Robustness" (upper panels) shows how ILLUSTRISTNG and SIMBA are configured in latent space $\tilde{Z}_{UMAP}$. The most noticeable separation between ILLUSTRISTNG and SIMBA is observed in the upper $\Omega_m$ (right panel), as also seen in Figure 2. The ILLUSTRISTNG (blue) and SIMBA (red) are proximal but distinctly segregated on a global scale. Meanwhile, clustering is also evident in the lower and intermediate $\Omega_m$, though it occurs on a more localized level, with some instances of intermixing noted. The AUC scores are 0.93 (lower), 0.98 (intermediate) and 0.99 (upper), indicating the distinctive features of the H1 maps for ILLUSTRISTNG and SIMBA manifest most in the higher $\Omega_m$. In contrast, the "With Robustness" panels (lower) exhibit no evident clustering or segregation, evidenced by AUC scores of 0.50 (lower), 0.49 (intermediate), and 0.53 (upper). Therefore, robustness implementation effectively integrates simulations.

## 3.2. Effect of Robustness in Unseen Simulations

We focus our investigation on the influence of robustness on data that the machine has not encountered previously. We introduce the unseen data sets, ASTRID and SWIFT-EAGLE, into the same machine utilized in the previous section, which has been trained on ILLUSTRISTNG and SIMBA. Figure 4 illustrates the configurations of the latent space $\tilde{Z}_{kPCA}$ for ASTRID (top) and SWIFT-EAGLE (bottom), superimposed on the latent variables of ILLUSTRISTNG and SIMBA. In contrast to the previous section, kernelPCA with the linear kernel is employed instead of UMAP, often yielding more interpretable results wherein the axes exhibit greater correlation with cosmological parameters (refer to Section 2.3.2). The objective

of this section is to scrutinize how ASTRID and SWIFT-EAGLE are arranged within the latent space $\tilde{Z}_{kPCA}$ spanned by MIEST trained on ILLUSTRISTNG and SIMBA, with a particular emphasis on robustness, especially in the context of cosmological parameters of interest for regression. Despite minor intricacies that are challenging to discern visually, the variations in ASTRID and SWIFT-EAGLE, in response to robustness, exhibit a similar pattern.

Without robustness (upper panels in each plot), ASTRID (green) and SWIFT-EAGLE (purple) are placed at the edge of the sheets that ILLUSTRISTNG (blue) and SIMBA (red) spans. The coverage of the unseen data (ASTRID and SWIFT-EAGLE) over the training data (ILLUSTRISTNG and SIMBA) is much less than 50% on average. This can introduce bias in the regression results since the locus of the latent space $\tilde{Z}_{kPCA}$ can roughly translate into the predicted values. However, since what we are looking at is the reduced space, the actual regression occurs in higher dimensions in a more complex manner. The third column that is supposedly related to classification of simulations demonstrates the evident distinction between ASTRID and the training data, whereas SWIFT-EAGLE seems to share more similarity with the training data, yet it is classifiable. We will revisit the similarity between simulations in Section 5.3

The incorporation of robustness results in a significant enhancement in coverage for both ASTRID and SWIFT-EAGLE (lower panels in each plot) compared to the case without robustness. Consequently, the majority of data points now align with the plane that ILLUSTRISTNG and SIMBA span on the reduced latent space $\tilde{Z}_{kPCA}$. Nonetheless, there remain biases and regions where ASTRID and SWIFT-EAGLE barely populate. In the case of $\Omega_m$ for SWIFT-EAGLE with robustness (lower-left corner), the distribution of the purple dots in the
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Two-dimensional projections of the 3D-reduced latent space $\tilde{Z}_{\text{kPCA}}$ using KernelPCA with the linear kernel. The MIEST machine is trained on ILLUSTRISTNG (blue) and SIMBA (red) and tested on ASTRID (green) and SWIFT-EAGLE (purple). The presented latent spaces are constructed from the whole data set of ILLUSTRISTNG and SIMBA including the training set and the test set, whereas the latent variables of ASTRID and SWIFT-EAGLE have not been seen by the MIEST machine. Each data point is color-coded by its true parameter value of $\Omega_m$ (first column) or $\sigma_8$ (second column). Due to robustness, the latent variables of ASTRID and SWIFT-EAGLE are spanned more evenly in the “With Robustness” case (lower panels) than in “No Robustness” (upper panels).</figcaption></figure>

high $\Omega_m$ region is notably sparse, with a discernible bias toward the middle and lower $\Omega_m$ values. This finding is also consistent with the prediction results depicted in Figure 5, which illustrates the predictions for SWIFT-EAGLE with robustness (details in Section 4). This indicates that these are

potential sources of systematic bias and inaccuracies in the prediction. Overall, the introduction of robustness demonstrates the capacity to enhance predictive accuracy by augmenting the congruence between training simulations and unobserved simulations. However, residual biases and
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Prediction of $\Omega_m$ vs. truth with and without robustness (lower and upper). The MIEST is trained on the H1 maps of ILLUSTRISTNG (first column) and SIMBA (second column) and tested on that of ASTRID (third column) and SWIFT-EAGLE (fourth column). Each plot has a relative error $\delta$, $\chi^2$, bias $b$, and $\alpha$. Here, $\alpha$ measures how the machine can cover the parameter space with less bias and sufficient confidence intervals (refer to Section 4).</figcaption></figure>

potential inaccuracies, which presumably stems from intrinsic properties of simulations, persist to a certain extent.

## 4. Performance

In this section, we dive into the robustness performance in an explicit manner. Figures 5 and 6 illustrate the predictions of $\Omega_m$ and $\sigma_8$, respectively, versus truth with and without robustness using ILLUSTRISTNG, SIMBA, ASTRID, and SWIFT-EAGLE. ILLUSTRISTNG (first column) and SIMBA (second column), on which the machine is trained, demonstrate the precision of the machine with respect to the robustness. We also provide in the figures the values of four metrics that measure the performance of MIEST: relative error $\delta$, chi-square $\chi^2$, bias $b$, and generalization $\alpha$. $\chi^2$ is defined by $\langle\theta_{\text{true}} - \hat{\theta}_{\text{pred}}|/\hat{\sigma}_{\text{pred}}\rangle$, where $\langle \cdot \rangle$ and $\hat{\sigma}_{\text{pred}}$ are the average over the test set the standard deviation given by MIEST (refer to Equations (9) and (10)). The generalization $\alpha$ is defined as $\langle b/\hat{\sigma}_{\text{pred}}\rangle$. The generalization $\alpha$ is designed to achieve a minimal value when the confidence interval of the prediction covers a broader parameter space with minimal bias; in other words, the model exhibits less bias and larger standard deviation compared to the bias. However, this metric does not directly measure the predictive power.

In the absence of robustness, the relative errors of the model for ILLUSTRISTNG (SIMBA) are quantified as 4.3%(5.3%) for $\Omega_m$ and 4.2%(4.9%) for $\sigma_8$. The slight difference in accuracy between two simulations originates from the intrinsic amount of information that each simulation carries, as discussed in Section 5.3. However, the introduction of robustness results in performance degradation, increasing the relative errors to 4.5% (5.6%) for $\Omega_m$ and 7.4%(7.9%) for $\sigma_8$. During the implementation of robustness, the machine is compelled to disregard

information exclusive to individual simulations, leading to an unavoidable decline in performance. However, the both values of bias and $\alpha$ decrease with robustness from $4 \times 10^{-4}$ to $-2 \times 10^{-4}$ (4.3 $\times 10^{-3}$ to $-2.7 \times 10^{-3}$) for bias and $1.5 \times 10^{-3}$ to $-1.2 \times 10^{-3}$ (0.14 to $-0.13$) for $\alpha$, indicating that the machine yields a more uniform prediction across the parameter space than without robustness. This suggests that the incorporation of robustness not only has the potential to enhance the generalization for unseen data but also augments the intrinsic robustness for the train data itself.

In the case of the unseen simulations, ASTRID (third column) and SWIFT-EAGLE (fourth column) show relative errors of 14.1% and 18.3% for $\Omega_m$ and 20.2% and 12.9% for $\sigma_8$, respectively. The substantial precision discrepancy between the train data and the unseen data can be attributed to the configurations of the latent space $\tilde{Z}_{\text{KPCA}}$ depicted in the “No Robustness” cases of Figure 4, which elucidates the evident bias present in ASTRID and SWIFT-EAGLE as opposed to ILLUSTRISTNG and SIMBA. However, by incorporating robustness, the relative errors of ASTRID and SWIFT-EAGLE have been improved by 23% and 11% for $\Omega_m$ and 47% and 29% for $\sigma_8$, respectively. Furthermore, $\chi^2$ is reduced by 2 ~ 3 times, with biases decreased by a factor of 7 for ASTRID and a factor of 4 for SWIFT-EAGLE, and the values of $\alpha$ are enhanced approximately four-fold for both. These improvements correlate with the alterations in the latent space $\tilde{Z}_{\text{KPCA}}$ attributed to robustness (see “With Robustness” of Figure 4). This substantiates the pivotal role of robustness in increasing performance.

In addition to performance enhancement, two noteworthy observations can be made. First, the improvement manifests itself with greater prominence in $\sigma_8$ compared to $\Omega_m$. In the absence of robustness, the predictive accuracy of $\sigma_8$ within
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Prediction of $\sigma_8$ vs. truth with and without robustness (lower and upper). The MIEST is trained on the H I maps of ILLUSTRISTNG (first column) and SIMBA (second column) and tested on that of ASTRID (third column) and SWIFT-EAGLE (fourth column). Each plot has a relative error $\delta$, $\chi^2$, bias $b$, and $\alpha$. Here, $\alpha$ measures how the machine can cover the parameter space with less bias and sufficient confidence intervals (refer to Section 4).</figcaption></figure>

both ASTRID and SWIFT-EAGLE exhibits substantial bias, largely yielding values in the range of 0.8–1.0 for $\sigma_8$, irrespective of the actual true value. However, the implementation of robustness remarkably mitigates such biases. The predictions have become increasingly aligned with the median values on average, indicating a reduction in bias, and exhibit a higher degree of correlation with the actual true values. On the other hand, the improvement in $\Omega_m$ is relatively minor compared to $\sigma_8$. Therefore, two conclusions can be drawn: (a) In the absence of robustness, $\sigma_8$ is substantially biased to their own simulations. (b) Nevertheless, $\sigma_8$ encompasses a substantial amount of both common and distinguishable information compared to $\Omega_m$, since $\Omega_m$ shows limited potential for improvement with respect to robustness.

The second noteworthy observation is that the prediction of $\Omega_m$ in SWIFT-EAGLE is significantly biased regardless of the presence of robustness. This underscores the inherent divergence in the physical models, as shown in the bottom-left panel of Figure 4, which illustrates the substantial absence of data points for SWIFT-EAGLE (purple) in the high $\Omega_m$ region both with and without robustness. The limitations associated with attaining robustness and divergence of the simulation models will be discussed in greater detail in Section 5.3.

We now investigate how performance can change with respect to the degree of robustness. As discussed in Section 2.2, the magnitude of de-classification and compression can be modulated by the hyperparameters $\beta$ and $\gamma$, respectively, in the loss function (see Equation (9)). De-classification and compression together play a significant role in robustness. Figure 7 shows how the relative error $\delta$ (left), bias $b$ (middle), and standard deviation $\sigma_{\text{pred}}$ (right) of $\Omega_m$ (top) and $\sigma_8$ (bottom) changes with respect to the degree of robustness, demonstrating qualitatively the relationship

between robustness and hyperparameters. There are four different levels of robustness on the x-axis: zero (($\beta, \gamma$) = (0, 0)), low (($\beta, \gamma$) = ($10^{-3}, 10^{-4}$)), mid (($\beta, \gamma$) = ($10^{-2}, 10^{-3}$)), and high (($\beta, \gamma$) = ($10^{-1}, 10^{-2}$)).

We can categorize the simulations into those used for training (ILLUSTRISTNG and SIMBA) and those reserved for testing (ASTRID and SWIFT-EAGLE). For the simulations employed in the training phase, the relationship between relative error and robustness is monotonic. An increase in robustness correlates with a reduction in predictive power due to loss of information through compression and de-classification processes. In principle, such a trend becomes manifested in the pursuit of robustness, as increments in $\gamma$ and $\beta$ result in gradient descent being increasingly influenced by compression $I$ and de-classification $H$, rather than the regression objective function $J_1$ (refer to Section 2.2.1). This phenomenon is conceptually understandable, as it involves the omission of information that, while distinguishable, may still exhibit correlation with cosmological parameters.

In the case of ASTRID and SWIFT-EAGLE, the accuracy in relation to the degree of robustness progressively enhances until reaching a turning point. This suggests the existence of one or more optimal points characterized by maximal robustness, necessitating an exhaustive parameter search to identify. Nevertheless, we have not conducted such exploration due to the complexity and high dimensionality inherent in the hyperparameters, which encompass those pertinent to the optimization of the neural network. The turnaround occurs relatively earlier for $\Omega_m$ than $\sigma_8$, indicating that $\Omega_m$ retains less common information or is more sensitive to compression and de-classification processes, thereby resulting in its premature degradation relative to $\sigma_8$. Furthermore, biases exhibit analogous patterns. Conversely, the standard deviations
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Relative errors of $\Omega_m$ (to panels) and $\sigma_8$ (bottom panels) vs. the degree of robustness. The MIEST is trained on the H I maps of ILLUSTRISTNG (blue dashed) and SIMBA (red dotted) and tested on that of ASTRID (green solid) and SWIFT-EAGLE (purple dotted-dashed).</figcaption></figure>

demonstrate unexpected outcomes. Typically, the loss of information, as evidenced by the relative error in ILLUSTRISTNG and SIMBA, would be associated with an increase in standard deviations. However, the observed standard deviations reveal a decline followed by abrupt increases at certain points.

# 5. Discussion

## 5.1. On Robustness

The definition of robustness varies depending on the context of the problem. This particular study is dedicated to the development of a machine learning model that demonstrates consistent performance across various cosmological simulations, particularly in the context of estimating cosmological parameters from H I maps. To this end, we have tuned the neural network to remove distinct information but preserve common information that is, importantly, correlated with cosmology. In this sense, achieving robustness in this work can be defined as finding common information, or latent space, that is both physical and universal to any simulation, with the important assumption that all simulations have physical integrity at a certain level. In this context, the latent variables are not required to be directly observable physical quantities, but rather should encompass any physical features of the simulations.

Imagine the feature space $\mathcal{H}_i$ of the physical features $h_i$ that can be found in simulation *i*. There can be a subset of features $C_i$ that are correlated with the parameter of interest $\theta$. That is, $C_i = \{h_i | I(h_i; \theta) > 0\}$, where the mutual information $I(X;Y) = D_{\text{KL}}(P_{(X,Y)}||P_X \otimes P_Y)$ and $D_{\text{KL}}$ is the Kullback-Leibler divergence. In the conventional, regular estimator, the objective is to find an optimal subset of features $M_i \subset C_i$. Then, our goal is to find a subset $\mathcal{R} = \{h|h \in \mathcal{A}, \mathcal{A} = \bigcap_i C_i\}$. Here, the subscript *i* is no longer needed in $\mathcal{R}$ since the features in $\mathcal{R}$ should be universal, i.e., in common between all

simulations. In other words, this subset represents the collection of robust physical features. We can also define a subset $\mathcal{R}_{ab} = \{h|h \in \mathcal{A}, \mathcal{A} = \bigcap_{i=\{a,b\}} C_i\}$. $\mathcal{R}_{ab}$ represents the subset of the robust features *h* that are common only in simulations *a* and *b*. Note that since the robust subset $\mathcal{R}$ comprises $C_i$, whose input features (in this work, H I map) are correlated with the target parameters (here, cosmological parameters), robustness is conditioned on given inputs and outputs; e.g., our model is only robust on cosmological parameters with respect to H I maps.

With this, we consider the case where we train a model on simulations [$a...m$] with robustness. Ideally, we can have a robust model that has learned $\mathcal{R}_{[a...m]}$ to the fullest. By definition, the simulations are not distinguishable on $\mathcal{R}_{[a...m]}$. In other words, on $\mathcal{R}_{[a...m]}$, the cross-entropy $H(P_a, P_b)$ should be the same as the entropy $H(P_a)$ and $H(P_b)$ for any two simulations *a* and *b* (i.e., $H(P_a, P_b) = H(P_a) = H(P_b)$). In practice, however, it is nearly impossible to obtain solely $\mathcal{R}_{[a...m]}$ as whole. Rather, during training, it is more likely to get $\mathcal{R}'_{[a...m]} \cup \mathcal{H}'_{[a...m]}$ where $\mathcal{R}'_{[a...m]} \subset \mathcal{R}_{[a...m]}$ and an arbitrary subset $\mathcal{H}'_{[a...m]} = \bigcup_{i=\{a...m\}} \mathcal{H}'_i$ for $\mathcal{H}'_i \subset \mathcal{H}_i$ that causes noise. Then, the leakage $\mathcal{H}'_{[a...m]}$ eventually lead to classifiability, namely $H(P_a, P_b) > H(P_a)$ since the features $h_i' \in \mathcal{H}_i'$ are distinguishable across different simulations. This phenomenon is observed in the AUC score trained with robustness.

In addition, we have used relative errors to measure the performance of MIEST. For instance, the relative error on ASTRID is related to the mutual information between the intersection—$\mathcal{R}_{\text{[IllustrisTNG,SIMBA]}} \cap \mathcal{R}_{\text{[IllustrisTNG,SIMBA,Astrid]}}$— and cosmological parameters. In contrast to the AUC score that accounts for distinguishability, the relative error on the unseen simulations has a theoretical upper bound, which is determined by the mutual information between the intersection—$\mathcal{R}_{\text{[IllustrisTNG,SIMBA]}} \cup C_{\text{Astrid}}$—and cosmological parameters.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Probability of a simulation model *a* given latent variables of a certain model—*Z*<sub>b</sub>, *P*(*a*|*Z*<sub>b</sub>), and AUC score *with* and *without* robustness. *P*(*a*|*Z*<sub>b</sub>) represents the probability of MIEST, which is trained on ILLUSTRISTNG and SIMBA, identifying a given input as either ILLUSTRISTNG or SIMBA. With *P*(IllustrisTNG|*Z*<sub>b</sub>) ~ 0.5 ~ *P*(SIMBA|*Z*<sub>b</sub>), robustness ensures that the MIEST identifies the input to be either ILLUSTRISTNG or SIMBA regardless of its origin.</figcaption></figure>

The intersection represents the degree of overlap among ILLUSTRISTNG, SIMBA, and ASTRID, reflecting the inherent properties of simulations irrespective to MIEST.

Furthermore, during the training phase, the fact that MIEST lacks information about ASTRID results in unstable outputs. In practice, $\mathcal{R}'$ of the MIEST—a subset of $\mathcal{R}[\text{IllustrisTNG,SIMBA}]$—is determined in a stochastic manner due to inherent randomness or imperfections in the training processes. The inherent stochasticity in the training processes introduces randomness into the intersection between $\mathcal{R}[\text{IllustrisTNG,SIMBA}]$ and $\mathcal{R}[\text{IllustrisTNG,SIMBA,Astrid}]$, which in turn determines the performance on ASTRID. This phenomenon arises since MIEST lacks information about where ASTRID should lie in $\mathcal{R}[\text{IllustrisTNG,SIMBA}]$. Consequently, this randomness inevitably results in inherent instability and divergence in the outcomes of MIEST or robustness studies, which is discussed in detail next.

## 5.2. Convergence and Stability of Robustness

In this section, we investigate the convergence as well as the stability of the robustness of our models. In this context, convergence is described as the phenomenon that the metrics accounting for robustness—e.g., AUC score, and performance on ASTRID and SWIFT-EAGLE—progressively approach a stable point as the self-performance of MIEST on the training simulations improves. Similarly, stability refers to the variance of these robustness metrics in comparison to those of the machine without robustness. For instance, if the variance of the machine with robustness is so large to the extent that it overlaps with the machine without robustness, then it may be considered highly unstable. Throughout this section, we use MIEST that is trained on ILLUSTRISTNG and SIMBA.

Figure 8 displays the changes in distinguishability of the latent variables *with* (bottom) and *without* (top) robustness

using AUC score (right column) and evidence probabilities (other columns). The evidence probabilities in Bayesian statistics refer to the probability of observing the data given all possible values of the parameters within a model. In this context, a probability of simulation models is defined as the evidence. We present the outputs of over 100 MIESTs, each trained individually and selected based on a mean relative error of $\Omega_m$ and $\sigma_8$ below 6%, represented by the blue dots. The evidence probability $P(a|Z_b)$ represents the probability of an arbitrary simulation model *a* given the latent variables of an arbitrary model, $Z_b$. In other words, this value represents the probability of a given latent variable $Z_b$—generated by the H I maps of simulation *b*—being classified as *a*. To compute this, we train a random forest classifier using latent variables *Z* and simulation labels *a*. The random forest subsequently estimates the probability of simulation models as a function of the latent variables.

"No Robustness" (upper panels) shows an AUC score of ~0.84 and evidence probabilities of $P(A|A)$ ~ 0.75 and $P(B|A)$ ~ 0.25, where $P(A|A)$ and $P(B|A)$ represent the cases where two simulations are the same or are different, respectively. This indicates that the "MIEST" without robustness can classify the simulations, yet there is some overlap of around 25% in terms of the evidence probability. The standard deviations of these values are approximately over 10%, which is not very stable. For example, some cases demonstrate a probability of 0.4 or 0.6, indicating that ILLUSTRISTNG and SIMBA are not distinguishable to some extent. Furthermore, there are no evident signs of convergence in the probabilities and the AUC scores toward smaller mean relative errors. Consequently, training MIEST without robustness yields unstable and divergent outputs in terms of distinguishability between ILLUSTRISTNG and SIMBA.
<table><caption>**Table 1**
Summary of the Relative Error $\delta$, Bias $b$, and Generalization $\alpha$ on ASTRID and SWIFT-EAGLE</caption><thead><tr><th rowspan="2">Robust.</th><th rowspan="2"></th><th colspan="2">ASTRID</th><th colspan="2">SWIFT-EAGLE</th></tr><tr><th>Without</th><th>With</th><th>Without</th><th>With</th></tr></thead><tbody><tr><td rowspan="3">Ω<sub>m</sub></td><td>δ</td><td>15.7 ± 2.8</td><td><b>14.1 ± 2.6</b></td><td><b>16.9 ± 2.0</b></td><td>17.9 ± 2.4</td></tr><tr><td>10b</td><td>1.5 ± 0.3</td><td><b>1.1 ± 0.4</b></td><td>1.5 ± 0.3</td><td><b>1.2 ± 0.5</b></td></tr><tr><td>α</td><td>6.1 ± 1.6</td><td><b>4.9 ± 1.8</b></td><td>6.7 ± 1.5</td><td><b>4.75 ± 2.4</b></td></tr><tr><td rowspan="3">σ<sub>8</sub></td><td>δ</td><td>14.5 ± 1.9</td><td><b>12.8 ± 0.5</b></td><td>12.7 ± 1.0</td><td><b>11.8 ± 1.0</b></td></tr><tr><td>10b</td><td>1.3 ± 0.3</td><td><b>0.2 ± 0.1</b></td><td>0.9 ± 0.3</td><td><b>0.2 ± 0.1</b></td></tr><tr><td>α</td><td>3.1 ± 0.6</td><td><b>0.2 ± 0.1</b></td><td>2.1 ± 0.6</td><td><b>0.15 ± 0.10</b></td></tr></tbody></table>

**Note.** The values presented are averages across multiple MIESTs trained with and without robustness to account for the effect of randomness in robustness. Lower numbers indicate superior performance. The superior values among “No Robustness” and “With Robustness” are highlighted in bold.

Meanwhile, “With Robustness” (lower panels) shows consistent results across different MIESTs, demonstrating the AUC score of ~0.51 and $P(A|A) \sim 0.51$, and $P(B|A) \sim 0.49$. This evidently indicates the robustness succeeds in removing the distinct model-specific information so that they are indistinguishable in the latent variables. Moreover, the standard deviations of the probabilities are <1% while that of the AUC score is ~2%, which is a huge improvement from “No Robustness.” This also can evidently be considered convergent. Overall, the implementation of robustness results in the stable, convergent MIEST.

We now move on to the performance on the unseen simulations—ASTRID and SWIFT-EAGLE. Table 1 summarizes the mean of the relative errors $\delta$, biases $b$, and generalization $\alpha$ and their standard deviations for $\Omega_m$ and $\sigma_8$ (refer to the first paragraph of Section 4 for the definitions). For ASTRID, it is evident that robustness improves the performance in every metric by more than 10% in terms of the mean values. However, for $\Omega_m$ of ASTRID, the standard deviations are more than 10%, leading to overlaps between “No Robustness” and “With Robustness.” This suggests that the improvements in $\Omega_m$ exhibit instability, implying that despite certain improvements, the influence of randomness is still significant. In contrast, the improvements in $\sigma_8$ are more significant and stable. Notably, the generalization performance is enhanced by a factor of 15, in line with Figure 6.

The case of SWIFT-EAGLE displays a similar tendency to ASTRID, albeit with relatively fewer overall improvements. A key difference in $\Omega_m$ is the inconsistency where the incorporation of robustness results in a higher relative error, deteriorating predictive performance; however, the robustness leads to improvements in the bias and generalization, despite being unstable. This demonstrates that while robustness in MIEST results in a loss of predictive power accompanied by increased scatter, it in return enhances generalization, enabling a more uniform coverage of the parameter space with heightened levels of confidence. This phenomenon is illustrated in Figure 4, which demonstrates that robustness enhances the coverage of SWIFT-EAGLE, albeit at the cost of inducing misalignments.

The convergence and (in-)stability arise from both the inherent characteristics of various simulations and the influence of stochastic elements in the training process. Ideally, if one can train a machine such that it removes

all of the distinct features, yet retaining common features, i.e., the latent space of the machine represents exactly $\mathcal{R}_{\text{IllustrisTNG,SIMBA}}$ in the previous section, this machine is considered convergent and stable, providing consistent results. In practice, however, due to stochasticity and incompleteness of the training process, the latent space results in representing only a subset of $\mathcal{R}_{\text{IllustrisTNG,SIMBA}}$, leading to potential issues in convergence and stability. Furthermore, the performance on ASTRID (SWIFT-EAGLE) is determined by the overlap between the intersection of ILLUSTRISTNG and SIMBA, $\mathcal{R}_{\text{IllustrisTNG,SIMBA}}$, and ASTRID (SWIFT-EAGLE). However, since the training lacks the information about ASTRID and/or SWIFT-EAGLE, the extent of the overlap remains completely unaccounted for in the training process; i.e., it is determined stochastically. Consequently, the performance on unseen simulations tends to be unstable and inherently unable to converge, posing significant hurdles in the robustness studies.

### 5.3. Different Combinations Of Simulations

Thus far, we have conducted this study using the machine only trained on ILLUSTRISTNG and SIMBA. In this section, we extend our scope to all four simulations and investigate differences between all possible pairs of simulations. To this end, we use three approaches; (1) We train a MIEST machine on two simulations without robustness to assess the degree to which the two simulations are distinguishable; (2) We individually train CNNs on each simulation and evaluate their performances on the other three simulations to obtain a comprehensive understanding of cross-simulation performance; (3) We perform a comparative analysis of the HI maps derived from the four simulations based on the results above.

First, we train MIEST without robustness on every combination of any two simulations out of the four simulations—ILLUSTRISTNG, SIMBA, ASTRID, and SWIFT-EAGLE, which gives us six combinations in total. Due to the lack of robustness, the MIEST is trained in a way that both compresses the maps and increases the mutual information between latent variables and cosmological parameters. Since the de-classification is not imposed, there is no need for MIEST to remove the classifiable information. For instance, this is the case in Section 3 where MIEST without robustness shows that the latent variables from ILLUSTRISTNG and SIMBA are apparently classifiable. Thus, we then adopt an evidence probability of the model to measure the distinguishability from a binary classifier. In this process, we train around 100 MIEST for each combination and calculate the mean and the standard deviation to get rid of the randomness in the result. Second, we individually train simple CNNs on each simulation and then evaluate the cross-performances by testing each CNN on the three simulations other than the one used for training. To mitigate the effects of randomness, we again repeat it over 100 times, and obtain the mean and the standard deviation for each combination.

Table 2 summarizes the probability of a simulation model A given latent variables of a certain model—$Z_B$, $P(A|Z_B)$, and AUC score without robustness. $P(A|Z_B)$ represents the probability of MIEST that is trained on Simulation A and Simulation B identifying a given input as either one of those simulations. Here, we present the results for all six combinations of the four simulations. Additionally, the corresponding latent space representation is shown in Figure 9. Table 2
<table><caption>**Table 2**
Summary of Probability of a Simulation Model A Given Latent Variables of a Certain Model—Z<sub>B</sub>, P(A|Z<sub>B</sub>), and AUC Score without Robustness</caption><thead><tr><th>Simulation A</th><th>Simulation B</th><th>P(B|Z<sub>A</sub>)</th><th>P(A|Z<sub>A</sub>)</th><th>P(B|Z<sub>B</sub>)</th><th>P(A|Z<sub>B</sub>)</th><th>AUC Score</th></tr></thead><tbody><tr><td>ILLUSTRISTNG</td><td>SIMBA</td><td>0.28 ± 0.08</td><td>0.72 ± 0.08</td><td>0.72 ± 0.08</td><td>0.28 ± 0.09</td><td>0.83 ± 0.07</td></tr><tr><td>ILLUSTRISTNG</td><td>ASTRID</td><td>0.45 ± 0.02</td><td>0.55 ± 0.02</td><td>0.51 ± 0.02</td><td>0.49 ± 0.02</td><td>0.58 ± 0.04</td></tr><tr><td>ILLUSTRISTNG</td><td>SWIFT-EAGLE</td><td>0.45 ± 0.04</td><td>0.55 ± 0.04</td><td>0.52 ± 0.04</td><td>0.48 ± 0.04</td><td>0.59 ± 0.07</td></tr><tr><td>SIMBA</td><td>ASTRID</td><td>0.30 ± 0.06</td><td>0.70 ± 0.06</td><td>0.66 ± 0.07</td><td>0.34 ± 0.07</td><td>0.81 ± 0.08</td></tr><tr><td>SIMBA</td><td>SWIFT-EAGLE</td><td>0.29 ± 0.06</td><td>0.71 ± 0.06</td><td>0.71 ± 0.06</td><td>0.29 ± 0.06</td><td>0.82 ± 0.05</td></tr><tr><td>ASTRID</td><td>SWIFT-EAGLE</td><td>0.47 ± 0.02</td><td>0.53 ± 0.02</td><td>0.52 ± 0.02</td><td>0.48 ± 0.02</td><td>0.57 ± 0.04</td></tr></tbody></table>

Note. P(A|Z<sub>B</sub>) represents the probability of MIEST, which is trained on Simulation A and Simulation B, identifying a given input as either Simulation A or Simulation B. With P(B|Z<sub>A</sub>) ~ 0.5 ~ P(A|Z<sub>B</sub>), robustness ensures that the MIEST identifies the input to be either Simulation A or Simulation B regardless of its origin. This analysis indicates that ILLUSTRISTNG, ASTRID, and SWIFT-EAGLE share more common features, whereas SIMBA is an outlier. The latent spaces of this result can be found in Figure 9.

remarkably exhibits consistent results in terms of the relationships between the four simulations. Specifically, ILLUSTRISTNG, ASTRID, and SWIFT-EAGLE share more common features, whereas SIMBA is an outlier. All combinations of two simulations among ILLUSTRISTNG, ASTRID, and SWIFT-EAGLE show probabilities P(B|Z<sub>A</sub>) of ~0.5 and AUC scores of ~0.6, while all combinations that include SIMBA exhibit probabilities P(B|Z<sub>A</sub>) of ~0.3 and AUC scores of ~0.8 (see Figure 9 for an extensive latent-space comparison of all four simulations).

One interesting observation is that despite the consistency of these trends, there exists a notable asymmetry. For instance, P(IllustrisTNG|Z<sub>Astrid</sub>) = 0.49 while P(Astrid|Z<sub>IllustrisTNG</sub>) = 0.45, and P(SIMBA|Z<sub>Astrid</sub>) = 0.34 while P(Astrid|Z<sub>SIMBA</sub>) = 0.30. Despite the comparable sizes of the variances (~0.02 for the ILLUSTRISTNG-ASTRID pair and ~0.06 for the SIMBA-ASTRID pair), the presence of these asymmetric features suggests that certain simulations are relatively more distinctive, making them more challenging to blend effectively. In this spirit, Σ<sub>A</sub>P(A|Z<sub>B</sub>) can be a proxy for measuring how general or agreeable the simulation B is with other simulations. Then, Σ<sub>A</sub>P(A|Z<sub>Astrid</sub>) shows the highest value of 0.43, whereas Σ<sub>A</sub>P(A|Z<sub>SIMBA</sub>) exhibits the lower value of 0.29. In other words, this implies that ASTRID possesses more common features and is more agreeable, while SIMBA stands out as the most distinctive simulation model.

In consideration of this, we proceed to analyze the projected H1 maps to gain a direct insight into the analysis above. Figures 10 and 11 display projections of the H1 maps from the four simulations with variations of cosmological parameters. The initial conditions of the simulations are designed to yield the same large-scale structure across different simulations. Figure 10 depicts the H1 maps with different values of Ω<sub>m</sub> with the other parameters fixed—Ω<sub>m</sub> = [0.1, 0.2, 0.3, 0.4, 0.5] from the top, while Figure 11 illustrates the H1 maps with different values of σ<sub>8</sub> with the other parameters fixed—σ<sub>8</sub> = [0.6, 0.7, 0.8, 0.9, 1.0] from the top.

The large-scale structure of all of the presented H1 maps—commonly known to be determined by physics of dark matter and the initial conditions of cosmological simulations—shows a great similarity across different simulations. However, appreciable differences in the general features are observable with ease. The primary distinction is the amount of ionized diffuse gas. In every map presented here, it is readily apparent that SIMBA contains a significantly greater amount of ionized diffuse gas—dark blue regions—than the other simulations,

with this disparity becoming more prominent as Ω<sub>m</sub> and σ<sub>8</sub> increase. This phenomenon is also completely in agreement with Figure 2, which demonstrates significantly increasing distances between ILLUSTRISTNG and SIMBA in the latent space Z̃<sub>UMAP</sub> as Ω<sub>m</sub> and σ<sub>8</sub> increase. On the other hand, ILLUSTRISTNG, ASTRID, and SWIFT-EAGLE exhibit similarity in terms of hot diffuse gas. This similarity accounts for the indistinguishability among ILLUSTRISTNG, ASTRID, and SWIFT-EAGLE, while also elucidating the distinguishability with SIMBA, shown in Table 2.

This disparity is also demonstrated in M. T. Tillman et al. (2023), who studied the impact of galactic feedback on Lyα forest statistics, which is characterized by column densities of H1, resulting from parameter variations in the SIMBA and ILLUSTRISTNG feedback models. The study highlights that in SIMBA, the AGN jet feedback mode effectively transports energy to the diffuse intergalactic medium, and stellar feedback is crucial in regulating SMBH growth and feedback, emphasizing the need to constrain both stellar and AGN feedback. In contrast, the ILLUSTRISTNG simulations show that AGN feedback variations do not influence the Lyα forest, but changes in the stellar feedback model lead to subtle effects (see M. T. Tillman et al. 2024; for a further discussion on SIMBA). This result elucidates the intrinsic discrepancies in the galactic feedback models in ILLUSTRISTNG and SIMBA, which results in the distinguishability.

On the other hand, the indistinguishability—e.g., among ILLUSTRISTNG, ASTRID, and SWIFT-EAGLE—does not necessarily result in better cross-performance. It is worth noting that the cross-performance depends on mutual information between the cosmological parameters and the H1 maps, whereas classification occurs on each map. Thus, the cross-performance is sensitive to how structures, such as the clustering of galaxies, changes with respect to the cosmological parameters, whose relationship is not readily discernible to the naked eye. For instance, SWIFT-EAGLE in comparison to ILLUSTRISTNG and ASTRID displays smoother structures, while ILLUSTRISTNG and ASTRID exhibit denser clumping, albeit this density is not always pronounced. In addition, such differences are not solely driven by the cosmology but also driven by the baryon physics even without varying the subgrid model parameters. For example, increasing the matter content in the Universe (high Ω<sub>m</sub>) results in the formation of dense clusters and galaxies, which facilitates the growth of SMBHs and subsequently leads to a greater fraction of hot diffuse gas in the Universe due to the
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Latent-space visualization using UMAP with the “canberra” metric from the six MIESTs trained on all possible combinations of two simulations without robustness. The same machines used in generating Table 2 are employed (see Section 5.3 for details).</figcaption></figure>

jet feedback from these black holes. Consequently, this can
produce effects comparable to the increase in the AGN
parameters. In other words, the estimation of cosmological
parameters can be significantly influenced by the subgrid
models utilized in cosmological simulations. This underscores
the necessity for robustness in these simulation models.

6. Summary

The robustness of cosmological parameter inference across
different simulation frameworks has garnered significant
attention in recent research. For instance, F. Villaescusa-Nav-
arro et al. (2021a) demonstrated that total matter density maps
exhibit a high degree of robustness: neural networks trained on
<figure><img src="image_13.png" /><figcaption>**Figure 10.** H I maps from ILLUSTRISTNG, SIMBA, ASTRID, and SWIFT-EAGLE with variation of Ω_m with the other parameters fixed. The initial conditions of the simulations are designed to yield the same large-scale structure across different simulations. Ω_m has been set to 0.1 (first row), 0.2 (second row), 0.3 (third row), 0.4 (fourth row), and 0.5 (fifth row). The simulation run with Ω_m = 0.3 (middle) is the fiducial CAMELS simulation.</figcaption></figure>

these maps from one simulation suite (e.g., ILLUSTRISTNG or SIMBA) can accurately recover cosmological parameters when tested on maps from the other. In addition, H. Shao et al. (2022) employed graph neural networks trained on halo

catalogs to infer cosmological parameters. Their findings suggest that velocity fields offer more robust estimations of parameters like Ω_m and σ_8 compared to halo positions. However, these quantities are not essentially direct observable.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** H I maps from ILLUSTRISTNG, SIMBA, ASTRID, and SWIFT-EAGLE with variation of $\sigma_8$ with the other parameters fixed. The initial conditions of the simulations are designed to yield the same large-scale structure across different simulations. $\sigma_8$ has been set to 0.6 (first row), 0.7 (second row), 0.8 (third row), 0.9 (fourth row), and 1.0 (fifth row). The simulation run with $\sigma_8$ = 0.6 (middle) is the fiducial CAMELS simulation.</figcaption></figure>

However, when it comes to baryon maps such as gas temperature maps, which can be reconstructed from direct observables (e.g., X-ray and thermal SZ observations), discrepancies become more pronounced. Y. Ni et al. (2023)

investigated the estimation of cosmological parameters using total mass and gas temperature maps derived from different simulations, including ILLUSTRISTNG, SIMBA, and ASTRID. Despite the consistency in $\Omega_m$ estimation from the total mass
maps with an appreciable bias in $\sigma_8$, their study revealed significant inconsistencies in parameter estimations from gas temperature maps across these simulations, underscoring the challenges in achieving robustness when relying on baryon maps, which are proxies for direct observable and sensitive to subgrid models (see Figures 9 and 11 in Y. Ni et al. 2023).

To address these challenges, we develop the MIEST that can robustly estimate the cosmological parameters $\Omega_m$ and $\sigma_8$ from H I maps of various cosmological simulation models in the CAMELS simulations. The cosmological simulation models employed in this study include ILLUSTRISTNG, SIMBA, ASTRID, and SWIFT-EAGLE. Our primary objective is to train MIEST on ILLUSTRISTNG and SIMBA with and without robustness, followed by testing on ASTRID and SWIFT-EAGLE to evaluate the robustness of the model. Here, an estimator is considered robust if it possesses a consistent predictive power across all simulations, including those used during the training phase.

To achieve robustness, we construct the MIEST by integrating the approaches of deep variational IB (A. A. Alemi et al. 2017) and generative adversarial network (I. Goodfellow et al. 2014). The structure of MIEST comprises a CNN, an encoder, a regressor, and a classifier. The CNN and an encoder image-process and compress the H I maps into latent variables, which can be deemed as summary statistics of these maps. With the latent variables, the MIEST can predict cosmological parameters via the regressor. During the training phase, the MIEST aims to discard the classifiable information specific to each simulation through de-classification process with the classifier, while preserving the information about the cosmological parameters (refer to Section 2.2 for details). Subsequently, the hyperparameters and training procedures are optimized using optuna to improve the performance of the machine.

We first investigate the latent space—a set of summary statistics within neural networks that captures the essential patterns or features of complex data—of the MIEST to understand the mechanism through which MIEST attains robustness and to assess the efficacy of the MIEST (see Section 3). We utilize methods for reducing dimensionality that provide meaningful evaluations of high-dimensional data, condensing the latent space to three dimensions. On the 3D latent space derived from the test set of ILLUSTRISTNG and SIMBA, MIEST without robustness generates two separate sheets of latent variables originating from ILLUSTRISTNG and SIMBA, respectively. The two axes of the sheets correspond to $\Omega_m$ and $\sigma_8$. However, the introduction of robustness results in these two sheets transforming into two intermingled spheres, showing no sign of separation or distinctiveness, while preserving the two axes for $\Omega_m$ and $\sigma_8$. Consequently, the removal of the classifiable information about simulation models leads to the blending of latent variables.

We turn to examine how robustness facilitates the incorporation of the training simulations with the unseen simulations. To this end, we focus on the overlapping area in the latent space between the training simulations—ILLUSTRISTNG and SIMBA—and the unseen simulations—ASTRID and SWIFT-EAGLE. The implementation of robustness increases coverage of ASTRID and SWIFT-EAGLE on ILLUSTRISTNG and SIMBA, leading to improved alignment of the axes for $\Omega_m$ and $\sigma_8$. This potentially indicates better performance on ASTRID and SWIFT-EAGLE. Nonetheless, biases persist, particularly in specific regions such as high $\Omega_m$ areas for SWIFT-EAGLE with

robustness, suggesting that systematic biases and inaccuracies still exist despite the improvements. Consequently, the incorporation of robustness helps align training and unseen simulations, but some persistent biases are likely due to the inherent properties of the simulations.

Subsequent to the qualitative analysis, we transition to a quantitative assessment of the robustness and performance for MIEST (refer to Section 4). For an estimator trained on IllustrisTNG and SIMBA without robustness, the mean relative errors of $\Omega_m$ and $\sigma_8$ are high for ASTRID (~17%) and SWIFT-EAGLE (~15%) compared to ILLUSTRISTNG (~4%) and SIMBA (~5%), showing clear bias. Incorporating robustness improves performance of ASTRID and SWIFT-EAGLE considerably, reducing the mean relative errors of $\Omega_m$ and $\sigma_8$ by ~35% (ASTRID) and ~20% (SWIFT-EAGLE) and decreasing biases by factors of 7 (ASTRID) and 4 (SWIFT-EAGLE). Although not shown here, preliminary results we have obtained indicate that similar improvements are observed in other fields such as gas temperature and electron density. These improvements, attributed to enhanced coverage and alignment in the latent space, underscore the importance of robustness in improving consistent performance over various simulations.

In addition to the enhancement in performance, two observations can be highlighted. First, the improvement is more significant in $\sigma_8$ compared to $\Omega_m$. Without robustness, the accuracy of $\sigma_8$ in both ASTRID and SWIFT-EAGLE exhibits considerable bias, yielding values between 0.8 and 1.0 for $\sigma_8$, regardless of the actual true value. Nonetheless, the incorporation of robustness significantly alleviates these biases, showing a higher correlation with the true values. Conversely, the improvement in $\Omega_m$ is relatively minor compared to $\sigma_8$. The second is that the prediction of $\Omega_m$ in SWIFT-EAGLE is significantly biased regardless of the presence of robustness. This underscores the inherent divergence in the physical models, as shown in the latent space that illustrates the substantial absence of data points in the high $\Omega_m$ region both with and without robustness.

We then investigate the changes of performance of MIEST in response to the different degrees of robustness. The incorporation of robustness increases the errors for ILLUSTRISTNG and SIMBA but improves the accuracy for ASTRID and SWIFT-EAGLE. Specifically, the accuracy of $\sigma_8$ is more significantly improved by the introduction of robustness than $\Omega_m$. Furthermore, the performance variations with different robustness levels are analyzed using hyperparameters $\beta$ and $\gamma$ in Equation (9) in Section 2.2.1, which control the degrees of compression and de-classification. For the training simulations, increased robustness leads to reduced predictive power due to information loss from compression and de-classification, respectively. In the unseen simulations, the accuracy gained through robustness shows an initial increase, followed by a decline after reaching an optimal point of maximum robustness. The turning point for $\Omega_m$ occurs earlier than that for $\sigma_8$, indicating that $\Omega_m$ is more sensitive to compression and de-classification.

The primary strength of robustness lies in its ability to mitigate biases arising from specific physical or subgrid modeling assumptions—an essential advantage given the considerable uncertainties and lack of consensus in this domain. Cosmological hydrodynamic simulations naturally allow for direct comparisons with observational data, yet they inherently depend on complex physical models and subgrid prescriptions, which themselves often become substantial
sources of systematic bias. Without robustness, the MIEST achieved a relative error of approximately 4%–5%; however, introducing robustness increased this uncertainty to roughly 15%. This degradation in precision is inevitable due to the limited common information available across diverse simulation models. Nevertheless, it crucially ensures that the *ideal* MIEST remains free from systematic modeling biases and, in principle, allows for direct applicability to observational data without concerns about such biases.

## Acknowledgments

The Flatiron Institute is supported by the Simons Foundation. This work was supported by the Simons Collaboration on “Learning the Universe,” by NSF grant AST-2108678 and by Simons Foundation grant PD-CMB-00011009.

## ORCID iDs

Yongseok Jo https://orcid.org/0000-0003-3977-1761

Shy Genel https://orcid.org/0000-0002-3185-1540

Anirvan Sengupta https://orcid.org/0000-0002-1080-0744

Benjamin Wandelt https://orcid.org/0000-0002-5854-8269

Rachel Somerville https://orcid.org/0000-0002-6748-6821

Francisco Villaescusa-Navarro https://orcid.org/0000-0002-4816-0455

## References

Akiba, T., Sano, S., Yanase, T., Ohta, T., & Koyama, M. 2019, in Proc. 25th ACM SIGKDD Int. Conf. Knowledge Discovery & Data Mining, KDD '19 (New York: Association for Computing Machinery), 2623

Alemi, A. A., Fischer, I., Dillon, J. V., & Murphy, K. 2017, Int. Conf. Learning Representations (OpenReview) https://openreview.net/forum?id=HyxQzBceg

Andrianomena, S., & Hassan, S. 2025, Ap&SS, 370, 14

Anglés-Alcázar, D., Davé, R., Faucher-Giguère, C.-A., Özel, F., & Hopkins, P. F. 2017a, MNRAS, 464, 2840

Anglés-Alcázar, D., Faucher-Giguère, C.-A., Kereš, D., et al. 2017b, MNRAS, 470, 4698

Bird, S., Ni, Y., Di Matteo, T., et al. 2022, MNRAS, 512, 3703

Blanchard, G., Deshmukh, A. A., Dogan, U., Lee, G., & Scott, C. 2021, J. Mach. Learn. Res., 22, 1

Bondi, H. 1952, MNRAS, 112, 195

Booth, C. M., & Schaye, J. 2009, MNRAS, 398, 53

Borrow, J., Schaller, M., Bower, R. G., & Schaye, J. 2022, MNRAS, 511, 2367

Busillo, V., Tortora, C., Napolitano, N. R., et al. 2023, MNRAS, 525, 6191

Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 450, 1937

Dalla Vecchia, C., & Schaye, J. 2012, MNRAS, 426, 140

Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827

Davé, R., Thompson, R., & Hopkins, P. F. 2016, MNRAS, 462, 3265

de Santi, N. S. M., Shao, H., Villaescusa-Navarro, F., et al. 2023, ApJ, 952, 69

de Santi, N. S. M., Villaescusa-Navarro, F., Abramo, L. R., et al. 2025, JCAP, 2025, 082

Di Matteo, T., Springel, V., & Hernquist, L. 2005, Natur, 433, 604

Echeverri-Rojas, N., Villaescusa-Navarro, F., Chawak, C., et al. 2023, ApJ, 954, 125

Faucher-Giguère, C.-A. 2020, MNRAS, 493, 1614

Feng, Y., Di-Matteo, T., Croft, R. A., et al. 2016, MNRAS, 455, 2778

Finn, C., & Levine, S. 2018, arXiv:1710.11622

Goodfellow, I., Pouget-Abadie, J., Mirza, M., et al. 2014, in Advances in Neural Information Processing Systems 27 (NIPS 2014), The 28th Int. Conf. on Neural Information Processing Systems, ed. Z. Ghahramani et al. (Red Hook, NY: Curran Associates, Inc.) https://dl.acm.org/doi/10.5555/2969033.2969125

Haardt, F., & Madau, P. 2012, ApJ, 746, 125

Hahn, C., Lemos, P., Parker, L., et al. 2023, arXiv:2310.15246

Hassan, S., Villaescusa-Navarro, F., Wandelt, B., et al. 2022, ApJ, 937, 83

He, S., Li, Y., Feng, Y., et al. 2019, PNAS, 116, 13825

Hopkins, P. F. 2015, MNRAS, 450, 53

Hopkins, P. F., Kereš, D., Oñorbe, J., et al. 2014, MNRAS, 445, 581

Hopkins, P. F., & Quataert, E. 2011, MNRAS, 415, 1027

Huang, Z., Wang, H., Xing, E. P., & Huang, D. 2020, in Computer Vision - ECCV 2020: 16th European Conf., Glasgow, UK, August 23-28, 2020 Proceed., Part II (Berlin: Springer), 124

Jamieson, D., Li, Y., de Oliveira, R. A., et al. 2023, ApJ, 952, 145

Jeffrey, N., & Wandelt, B. D. 2020, arXiv:2011.05991

Jo, Y., Genel, S., Wandelt, B., et al. 2023, ApJ, 944, 67

Jo, Y., & Kim, J.-h. 2019, MNRAS, 489, 3565

Katz, N., Weinberg, D. H., & Hernquist, L. 1996, ApJS, 105, 19

Kim, J.-h., Abel, T., Agertz, O., et al. 2014, ApJS, 210, 14

Krumholz, M. R., & Gnedin, N. Y. 2011, ApJ, 729, 36

Lee, J.-Y., Kim, J.-h., Jung, M., et al. 2024, ApJ, 975, 38

Lee, M. E., Genel, S., Wandelt, B. D., et al. 2024, ApJ, 968, 11

Li, Y., Tian, X., Gong, M., et al. 2018, in Computer Vision ECCV 2018: 15th European Conf., Munich Germany, September 8-14, 2018, Proc., Part XV (Berlin: Springer), 647

McInnes, L., Healy, J., & Melville, J. 2020, arXiv:1802.03426

Medlock, I., Nagai, D., Singh, P., et al. 2024, ApJ, 967, 32

Muandet, K., Balduzzi, D., & Schölkopf, B. 2013, in Proc. 30th Int. Conf. Machine Learning (ICML 2013), Vol. 28, (PMLR), 10, https://proceedings.mlr.press/v28/muandet13.html

Nguyen, T., Modi, C., Yung, L. Y. A., & Somerville, R. S. 2024, MNRAS, 533, 3144

Ni, Y., Di Matteo, T., Chen, N., Croft, R., & Bird, S. 2022, ApJL, 940, L49

Ni, Y., Genel, S., Anglés-Alcázar, D., et al. 2023, ApJ, 959, 136

Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077

Ploeckinger, S., & Schaye, J. 2020, MNRAS, 497, 4857

Rahmati, A., Pawlik, A. H., Raičević, M., & Schaye, J. 2013, MNRAS, 430, 2427

Rosas-Guevara, Y. M., Bower, R. G., Schaye, J., et al. 2015, MNRAS, 454, 1038

Schaller, M., Gonnet, P., Draper, P. W., et al. 2018, SWIFT: SPH With Interdependent Fine-grained Tasking, Astrophysics Source Code Library, ascl:1805.020

Schaller, M., Borrow, J., Draper, P. W., et al. 2024, MNRAS, 530, 2378

Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521

Schaye, J., & Dalla Vecchia, C. 2008, MNRAS, 383, 1210

Schölkopf, B., Smola, A. J., Müller, K. R., et al. 1997, in Proc. 7th Int. Conf. Artificial Neural Networks (ICANN '97) (Berlin: Springer), 583

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Shao, H., de Santi, N. S. M., Villaescusa-Navarro, F., et al. 2023, ApJ, 956, 149

Shao, H., Villaescusa-Navarro, F., Villanueva-Domingo, P., et al. 2022, arXiv:2209.06843

Smith, B. D., Bryan, G. L., Glover, S. C. O., et al. 2017, MNRAS, 466, 2217

Springel, V. 2005, MNRAS, 364, 1105

Springel, V. 2010, ARA&A, 48, 391

Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289

Tillman, M. T., Burkhart, B., Tonnesen, S., Bird, S., & Bryan, G. L. 2025, ApJ, 980, 72

Tillman, M. T., Burkhart, B., Tonnesen, S., et al. 2023, AJ, 166, 228

Tishby, N., Pereira, F., & Bialek, W. 1999, in Proc. 37th Annual Allerton Conf. Communication, Control and Computing (Urbana-Champaign, IL: Univ. Illinois Press), 368

Villaescusa-Navarro, F., Anglés-Alcázar, D., Genel, S., et al. 2021a, ApJ, 915, 71

Villaescusa-Navarro, F., Anglés-Alcázar, D., Genel, S., et al. 2021b, arXiv:2109.09747

Villaescusa-Navarro, F., Genel, S., Angles-Alcazar, D., et al. 2021c, arXiv:2109.10360

Villaescusa-Navarro, F., Genel, S., Anglés-Alcázar, D., et al. 2022, ApJS, 259, 61

Villaescusa-Navarro, F., Genel, S., Castorina, E., et al. 2018, ApJ, 866, 135

Villanueva-Domingo, P., & Villaescusa-Navarro, F. 2022, ApJ, 937, 115

Vogelsberger, M., Genel, S., Sijacki, D., et al. 2013, MNRAS, 436, 3031

Vogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518

Wang, J., Lan, C., Liu, C., et al. 2023, ITKDE, 35, 8052

Weinberger, R., Springel, V., Hernquist, L., et al. 2017, MNRAS, 465, 3291

Weinberger, R., Springel, V., & Pakmor, R. 2020, ApJS, 248, 32

Wiersma, R. P. C., Schaye, J., Theuns, T., Dalla Vecchia, C., & Tornatore, L. 2009, MNRAS, 399, 574

Zhang, Y., Deng, B., Li, R., Jia, K., & Zhang, L. 2023, arXiv:2301.12643

Zhao, H., des Combes, R. T., Zhang, K., & Gordon, G. J. 2019, arXiv:1901.09453

Zhou, K., Liu, Z., Qiao, Y., Xiang, T., & Loy, C. C. 2023, ITPAM, 45, 4396