Celia: Catalog Embeddings Learned from cosmologIcal simulAtions

Introduction

Catalogs of galaxy positions and properties are precise probes of both the growth of cosmic structure and the physics of galaxy formation. The positions of galaxies trace the matter distribution, and their masses, sizes and stellar populations reflect how gas cools and forms stars within dark-matter halos (Wechsler & Tinker 2018). Modern surveys deliver such catalogs over large cosmological volumes. Spectroscopic surveys map the three-dimensional distribution of millions of galaxies

(DESI Collaboration et al. 2026; Abdul Karim et al. 2025), and imaging surveys measure weak gravitational lensing and galaxy clustering ([1]; DES Collaboration et al. 2026; [55]; Wright et al. 2025). Wide-area imaging with Euclid and the Rubin Observatory extends the reach of these measurements (Euclid Collaboration et al. 2025; Ivezić et al. 2019; Prat et al. 2023). Studies of nearby satellites and high-redshift galaxies probe the same galaxy-formation processes over different environments and epochs (Geha et al. 2017; Eisenstein et al. 2026).

These data constrain cosmology and galaxy formation, whose effects are entangled on nonlinear scales. There, stellar and active galactic nucleus (AGN) feedback redistribute the gas, changing both galaxy populations and the matter distribution measured by lensing [110] (Chisari et al. 2019). The resulting suppression of the matter power spectrum depends on the feedback model [111] (Schneider & Teyssier 2015) and affects comparisons between lensing and cosmic microwave background constraints (Amon & Efstathiou 2022; McCarthy et al. 2023). These scales also contain information beyond two-point statistics (Villaescusa-Navarro et al. 2020). Simulation-based analyses now use non-linear galaxy clustering and weak-lensing maps in survey data (Hahn et al. 2023; Jeffrey et al. 2025), making the treatment of baryonic physics an integral part of cosmological inference (Huang et al. 2019; Aricò et al. 2023). Hydrodynamical simulations impose physical priors on this entanglement, which makes them valuable for understanding galaxy formation and for constraining cosmology.

Hydrodynamical simulations follow structure formation and galaxy evolution together (Vogelsberger et al. 2020; Crain & van de Voort 2023). IllustrisTNG, EAGLE, SIMBA, Astrid and FLAMINGO implement star formation, black-hole growth and feedback through different subgrid prescriptions (Pillepich et al. 2018a,b; Schaye et al. 2015; Davé et al. 2019; Bird et al. 2022; Ni et al. 2022; Schaye et al. 2023). Selected parameters are calibrated against observables such as the galaxy stellar mass function and cluster gas fractions (Crain et al. 2015; Pillepich et al. 2018a; Kugel et al. 2023). From the same initial conditions, different models predict different galaxy masses and structures (Scannapieco et al. 2012). Their predictions also differ for the matter power spectrum and for the response of dark matter to baryons [111] (Gebhardt et al. 2026). CAMELS and DREAMS sample this model dependence by varying cosmological or dark-matter parameters together with feedback parameters across large simulation ensembles (Villaescusa-Navarro et al. 2021, 2023; Rose et al. 2025). CAMELS repeats these variations with several galaxy-formation models [70], providing catalogs in which the effects of cosmology and feedback can be studied across different physical implementations.

Parameter inference illustrates both the information available and the challenge of transferring it between models. Simulation-based inference connects simulations to data without requiring a tractable likelihood [43] (Cranmer et al. 2020; Alsing et al. 2019). Neural compression has been used to retain cosmological information while marginalizing over baryonic effects (Villaescusa-Navarro et al. 2022b). Parameter information has been extracted from two-dimensional maps (Villaescusa-Navarro et al. 2022a), individual galaxies and small galaxy populations (Villaescusa-Navarro et al. 2022c; Echeverri-Rojas et al. 2023; Chawlak et al. 2024; Liao et al. 2026), and luminosity functions and colors (Lovell et al. 2025). Transfer between simulation models depends on the input properties. Positions and velocities have supported robust inference of $\Omega_{\mathrm{m}}$ across hydrodynamical codes [109] (Shao et al. 2023), and internal galaxy properties can introduce stronger model dependence (Villaescusa-Navarro et al. 2022c; Echeverri-Rojas et al. 2023; Shao et al. 2023). Domain adaptation, domain generalization and training across several models improve this transfer (Roncoli et al. 2023; Lee et al. 2024; Jo et al. 2025).

Inference and emulation both need a description of an entire catalog, including its spatial arrangement, the properties of its members and the relations among them. Set and point-cloud networks provide such a description through pooling and attention (Zaheer et al. 2017; Qi et al. 2017; Vaswani et al. 2017; Lee et al. 2019; Zhao et al. 2021). Graph and point-cloud models already infer cosmological parameters from halo and galaxy catalogs [109] (Villanueva-Domingo & Villaescusa-Navarro 2022; Makinen et al. 2022; Anagnostidis et al. 2022; Shao et al. 2023; Chatterjee & Villaescusa-Navarro 2025; Baldwin et al. 2026). Generative models emulate catalogs, evaluate their likelihoods or populate dark-matter simulations with galaxies (Cuesta-Lazaro & Mishra-Sharma 2024; Nguyen et al. 2026; Pandey et al. 2025). Benchmarks such as CosmoBench bring several of these tasks together (Huang et al. 2025). Each of these networks addresses one task with one set of inputs, and what it extracts from the catalogs serves that task alone. A representation that could be reused across these analyses would make more of the information in each simulation available without learning a new catalog encoder for every new application.

Pretraining separates representation learning from downstream task (Bommasani et al. 2021; Brown et al. 2020). Self-supervised objectives learn from the data themselves by reconstructing masked inputs (Devlin et al. 2019; He et al. 2022; Pang et al. 2022), predicting representations of masked regions (Assran et al. 2023), matching augmented views (Grill et al. 2020), or aligning paired modalities (Radford et al. 2021). Such reuse is valuable for hydrodynamical simulations, whose computational cost limits their number, volume and diversity. Each CAMELS Latin-hypercube set used here contains 1,000 simulations of a $(25\,h^{-1}\,\mathrm{Mpc})^3$ box (Villaescusa-Navarro et al. 2021). Pretraining on N-body simulations has already reduced the number of hydrodynamical simulations needed for an inference task (Saoulis et al. 2025). In astronomy, pretrained image and spectrum representations support multiple measurements of galaxy properties (Hayat et al. 2021; Parker et al. 2024; Parker et al. 2025), and related approaches learn from simulated galaxies, cosmological fields and physical simulations (Xia et al. 2025; Akhmetzhanova et al. 2024; Kannan et al. 2025; McCabe et al. 2025; Mukhopadhyay et al. 2026). Ting (2026) reviews these developments.

Shared representations can organize the differences between physical models. To measure the baryonic suppression of the matter power spectrum, Lin et al. (2026) train a two-dimensional latent model jointly on the IllustrisTNG, Astrid and Swift-EAGLE suites. Pretrained representations also transfer between data sets and scientific tasks. Models of particle sets in collider physics can be reused across tasks (Golling et al. 2024; Mikuni & Nachman 2025), and Omni-Cosmos applies such a model to cosmological inference from halo and galaxy point clouds (Mikuni et al. 2026). In 21-cm cosmology, SKATR learns summaries that support parameter inference and can be transferred, with a frozen encoder, to a different simulator (Ore et al. 2025; Pietschke et al. 2026). For galaxy catalogs, the corresponding question is whether one representation can retain the responses to both cosmology and feedback across different galaxy-formation models, and make that information accessible to simple downstream readouts.

Here we present Celia, a self-supervised foundation model for simulated galaxy catalogs. We pretrain its encoder on the subhalo catalogs of the IllustrisTNG, SIMBA and Astrid suites of CAMELS and on their N-body counterparts, without simulation-parameter labels in the pretraining objective. The encoder learns the relations among subhalos and their properties by reconstructing the properties hidden from a masked view and by predicting the representations that a teacher network assigns to the same subhalos in the complete view. After pretraining, the frozen encoder produces a fixed-length sum mary of each catalog, and parameter inference requires only a linear readout. We test how well such readouts recover the cosmological and feedback parameters, how the summary is transferred to Swift-EAGLE, a code outside the training suites, and what the principal axes of the embedding reveal about each simulation suite. Section 2 presents the simulations, encoder, training and readouts. Section 3 reports the inference, transfer and representation analyses. Sections 4 and 5 discuss the implications and summarize the conclusions.

Method

Figure 1 outlines the method. We train an encoder (§ 2.3) on the subhalo catalogs of three CAMELS simulation suites (§ 2.1 and § 2.2) with a self-supervised objective that at no point uses labeled simulation parameters (§ 2.4). After training, each catalog is summarized by the covariance of the object representations that the frozen encoder assigns to its subhalos (§ 2.5). Linear readouts of this summary measure which parameters the representation encodes, and a second linear model assigns every prediction an error bar (§ 2.6). Appendices A–C give the input fields and the details of the encoder, its training and the readouts, and Appendix D tests the calibration of the error bars.

Simulations

We use the first generation of the Cosmology and Astrophysics with MachinE Learning Simulations (CAMELS; [70] (Villaescusa-Navarro et al. 2021; Villaescusa-Navarro et al. 2023)). Every hydrodynamical simulation follows $256^3$ dark-matter particles and initially $256^3$ gas particles from $z = 127$ to $z = 0$ in a periodic box of comoving side length $L = 25\,h^{-1}\,\mathrm{Mpc}$. The cosmology is flat $\Lambda$CDM with $\Omega_{\mathrm{b}} = 0.049$, $h = 0.6711$, $n_{\mathrm{s}} = 0.9624$, $w = -1$ and massless neutrinos. We use the snapshots at $z = 0$.

The four suites differ in their hydrodynamics code and galaxy-formation model (Table 1). The IllustrisTNG suite uses the moving-mesh code AREPO (Springel 2010) with the IllustrisTNG model (Weinberger et al. 2017; Pillepich et al. 2018a), and SIMBA uses the meshless finite-mass method of GIZMO (Hopkins 2015) with the SIMBA model (Davé et al. 2019). Astrid uses the smoothed-particle hydrodynamics code MP-Gadget (Feng et al. 2018) with the Astrid model ([69, 70] (Bird et al. 2022)), and Swift-EAGLE uses the smoothed-particle hydrodynamics code SWIFT (Schaller et al. 2024) with a galaxy-formation model based on EAGLE (Schaye et al. 2015). The models differ in their prescriptions for star formation, stellar feedback, black-hole growth and AGN feedback, and therefore populate halos of the same mass with different galaxies, gas and black holes.

	IllustrisTNG	SIMBA	Astrid	Swift-EAGLE
Code	AREPO	GIZMO	MP-Gadget	SWIFT
Hydrodynamics	Moving mesh	Meshless finite mass	SPH	SPH
$A_{\mathrm{SN1}}$	Wind energy per unit SFR	Wind mass loading	Wind energy per unit SFR	Energy per SN II event
$A_{\mathrm{SN2}}$	Wind speed	Wind speed	Wind speed	Metallicity dependence of stellar feedback energy
$A_{\mathrm{AGN1}}$	Kinetic-mode energy per unit accretion	Quasar- and jet-mode momentum flux	Kinetic-mode energy per unit accretion	Bondi accretion rate
$A_{\mathrm{AGN2}}$	Kinetic-mode ejection speed (burstiness)	Jet speed	Thermal-mode energy per unit accretion	Heating temperature of AGN feedback
Role	Training, readouts	Training, readouts	Training, readouts	Readouts only
LH catalogs	700 / 100 / 200	700 / 100 / 200	700 / 100 / 200	700 / 100 / 200

Each suite has a Latin-hypercube (LH) set of 1,000 simulations that sample six parameters jointly, and every simulation has its own random seed for the initial conditions. Two cosmological parameters vary, $\Omega_{\mathrm{m}} \in[0.1, 0.5]$ and $\sigma_8 \in[0.6, 1.0]$, together with four dimensionless amplitudes $A_{\mathrm{SN1}}$, $A_{\mathrm{SN2}}$, $A_{\mathrm{AGN1}}$ and $A_{\mathrm{AGN2}}$ that rescale the stellar and AGN feedback of each model relative to its fiducial value $A = 1$. The feedback parameters carry the same names in all suites but scale different processes (Table 1). $A_{\mathrm{AGN1}}$, for example, sets the energy of kinetic-mode AGN feedback in IllustrisTNG and Astrid, the momentum flux of quasar- and jet-mode feedback in SIMBA, and the Bondi accretion rate of black holes in Swift-EAGLE. $A_{\mathrm{SN1}}$ and $A_{\mathrm{AGN1}}$ vary within $[0.25, 4]$ and $A_{\mathrm{SN2}}$ and $A_{\mathrm{AGN2}}$ within $[0.5, 2]$, except $A_{\mathrm{AGN2}}$ in Astrid, which varies within $[0.25, 4]$.

Sections 3.5 and 3.7 use two further sets of IllustrisTNG, SIMBA and Astrid. In a one-parameter (1P) set, one parameter at a time takes its fiducial value and two values on either side of it, and all other parameters stay at their fiducial values. All 1P simulations share the random seed of their initial conditions. The cosmic-variance (CV) set keeps the fiducial parameters $\Omega_{\mathrm{m}} = 0.3$, $\sigma_8 = 0.8$ and $A = 1$, and its simulations differ only in the random seed of their initial conditions. For IllustrisTNG, we leave out the CV simulations that were run with a different code revision. The IllustrisTNG 1P simulations at the two values next to the fiducial one and at $\Omega_{\mathrm{m}} = 0.1$ were run with other code revisions, and we keep them in all 1P analyses. The 1P and CV catalogs are embedded in the same way as the LH catalogs of their suite.

A fixed random permutation of the LH indices assigns 700 simulations of each suite to training, 100 to validation and 200 to testing, and every suite uses the same permutation. We used the validation catalogs to compare encoder configurations during development and to set the confidence tiers of § 3.6. The readouts are fitted on the training catalogs, and the test catalogs enter only the evaluation. The self-supervised training uses the 2,100 training simulations of IllustrisTNG, SIMBA and Astrid and their 2,100 N-body counterparts. Swift-EAGLE is a held-out set. Its catalogs enter neither the training nor the normalization of the inputs. Its LH set is split in the same way and serves only to evaluate the frozen encoder (§ 3.4).

Catalogs and input fields

Each input object is a gravitationally bound subhalo identified by SUBFIND (Springel et al. 2001; Dolag et al. 2009) in the CAMELS catalogs at $z = 0$ (Villaescusa-Navarro et al. 2023), together with the properties of its host friends-of-friends (FoF) group. The subhalo-to-group relation supplies the host environments, and the first subhalo of each group identifies its central. We keep centrals and satellites with dark-matter mass

where $\Omega_{\mathrm{m}}^{\max} = 0.5$ is the largest matter density of the hypercube and $\rho_{\mathrm{crit},0} = 2.775 \times10^{11}\,h^2\,\mathrm{M}_{\odot}\,\mathrm{Mpc}^{-3}$ is the critical density. The threshold is the mass of 50 particles in an N-body simulation with $\Omega_{\mathrm{m}} = \Omega_{\mathrm{m}}^{\max}$, so every selected subhalo contains at least 50 dark-matter particles in both the hydrodynamical and the N-body catalogs. The threshold is a fixed mass, so the selection does not change with cosmology through the particle mass. The same threshold applies to all suites.

The inputs use a subset of the fields of the SUBFIND and FoF catalogs (Table 3 and Table 4 in Appendix § 7). Besides its position, each subhalo is described by 30 properties of itself and of its host group and by a flag that marks whether it is the central subhalo of its group. Some properties are taken directly from catalog fields, such as the stellar mass or $V_{\max}$, and others are derived from them, such as the mass fractions of the host and the position and velocity of the subhalo relative to its host. A binary positivity indicator accompanies each of 13 properties that can vanish, such as the stellar mass and the star-formation rate, so that a measured zero differs from a missing value. A subhalo thus enters the encoder with 53 input features. Three of them are the coordinates of the position, and the properties give 36 continuous features, with three components for each vector property. The 13 positivity indicators and the central

Catalog field	Dim.	Quantity	Transformation	Input group
Geometry and kinematics
SubhaloPos	3	Position	$(\mathbf{r}\bmod L)/L$	—
SubhaloVel	3	Velocity	$\mathbf{v}/10^{3}$	velocity
Subhalo dark matter and structure
SubhaloMassType[1]	1	Dark-matter mass	$\log_{10}q+\delta_M$	DM mass
SubhaloHalfmassRadType[1]	1	DM half-mass radius	$\log_{10}q$	half-mass radii
SubhaloHalfmassRad	1	Total half-mass radius $R_{1/2}$	$\log_{10}q$	half-mass radii
SubhaloVmax	1	Maximum circular velocity $V_{\max}$	$\log_{10}q$	dynamics
SubhaloVmaxRad	1	Radius of $V_{\max}$	$\log_{10}q$	dynamics
SubhaloVelDisp	1	Velocity dispersion	$\log_{10}q$	dynamics
$\lVert\mathrm{SubhaloSpin}\rVert$	1	Specific angular momentum	$\log_{10}q$	dynamics
$\lVert\mathrm{SubhaloPos}-\mathrm{SubhaloCM}\rVert$	1	Center offset $d$	$\ln(1+d/R_{1/2})$	center offset
Subhalo baryons
SubhaloMassType[4]	2	Stellar mass	$\log_{10}q$, $1_{q>0}$	stellar
SubhaloHalfmassRadType[4]	1	Stellar half-mass radius $R_{1/2,\star}$	$\log_{10}q$	stellar radius
SubhaloStarMetallicity	2	Stellar metallicity within $2R_{1/2,\star}$	$\log_{10}q$, $1_{q>0}$	stellar metallicity
SubhaloMassType[0]	2	Gas mass	$\log_{10}q$, $1_{q>0}$	gas
SubhaloGasMetallicity	2	Gas metallicity within $2R_{1/2,\star}$	$\log_{10}q$, $1_{q>0}$	gas metallicity
SubhaloGasMetallicitySfr	2	Metallicity of star-forming gas	$\log_{10}q$, $1_{q>0}$	SF-gas metallicity
SubhaloSFR	2	Star-formation rate	$\log_{10}q$, $1_{q>0}$	star formation
SubhaloSFRinRad	2	SFR within $2R_{1/2,\star}$	$\log_{10}q$, $1_{q>0}$	star formation
Subhalo black holes
SubhaloBHMass	2	Black-hole mass	$\log_{10}q$, $1_{q>0}$	black hole
SubhaloBHMdot	2	Black-hole accretion rate	$\log_{10}q$, $1_{q>0}$	BH accretion

Catalog field	Dim.	Quantity	Transformation	Input group
$\sum_k$ GroupMassType[$k$]	1	FoF mass $M_{\mathrm{FoF}}$	$\log_{10}q+\delta_{M,h}$	FoF mass
Group_M_Crit200	1	Halo mass $M_{200c}$	$\log_{10}q+\delta_{M,h}$	$M_{200c}$
Group_M_Crit500	1	Halo mass $M_{500c}$	$\log_{10}q+\delta_{M,h}$	$M_{500c}$
$\lVert\mathrm{GroupPos}-\mathrm{GroupCM}\rVert$	1	Center offset $d_h$	$\ln(1+d_h/R'_{200c})$	host center offset
GroupMassType[$k\ne1$]	1	Dark-matter fraction $f'_{\mathrm{DM}}$	$\log_{10}f'_{\mathrm{DM}}$	host DM fraction
GroupMassType[0]	2	Gas fraction	$\log_{10}(q/M'_{\mathrm{FoF}})$, $1_{q>0}$	gas
GroupMassType[4]	2	Stellar fraction	$\log_{10}(q/M'_{\mathrm{FoF}})$, $1_{q>0}$	stellar
GroupBHMass	2	Black-hole mass ratio	$\log_{10}(q/M'_{\mathrm{FoF}})$, $1_{q>0}$	black hole
GroupSFR	2	Star-formation rate	$\log_{10}q$, $1_{q>0}$	star formation
SubhaloPos − GroupPos	3	Position relative to the host	$\Delta r/R'_{200c}$	host-relative
SubhaloVel − GroupVel	3	Velocity relative to the host	$\Delta v/V'_{200c}$	host-rel. velocity
GroupFirstSub	1	Central flag	$1_{q=i}$, subhalo index $i$	host-relative

flag make up the remaining 14 features, which are binary. Masses, radii, rates and characteristic velocities such as $V_{\max}$ enter as logarithms, and all features except the positions are standardized with statistics of the training catalogs.

Every input feature carries an availability flag and a validity flag. A feature is unavailable when a catalog does not contain the underlying field, which is the case for the baryonic and black-hole properties of the N-body catalogs. It is invalid when its transformation is undefined, as for the logarithm of a vanishing mass, and the positivity indicator of such a mass remains valid. The encoder receives no label of the suite or of the kind of simulation. An N-body catalog differs from a hydrodynamical one only through its unavailable baryonic features and the properties of its subhalos.

The mass of a simulated subhalo is a multiple of the particle mass, which scales with $\Omega_{\mathrm{m}}$. We therefore omit the particle counts and perturb the base-10 logarithm $\ell_M$ of every mass that contains dark matter,

so that the encoder cannot read $\Omega_{\mathrm{m}}$ from the discreteness of the masses. Each offset is drawn once per object and stays the same in all training views and in the final embedding. The three host masses share one offset, and the host scales and fractions are computed from the perturbed masses.

Encoder

The encoder takes a view of a catalog with up to 1,024 subhalos as input. It returns a contextual state $\boldsymbol{c}_i \in\mathbb{R}^{256}$ for each subhalo, and these states serve as the object representations (Figure 1a). It builds a token for each subhalo from its visible features, adds spatial context from a periodic grid and from attention between all pairs of subhalos, and contextualizes the tokens with a Set Transformer (Lee et al. (2019)). All its operations are equivariant to permutations of the subhalos, so reordering the subhalos of a view reorders their states in the same way.

Let $\boldsymbol{f}_i \in\mathbb{R}^{53}$ collect the input features of subhalo $i$, and let $a_{ij}$, $v_{ij}$ and $o_{ij}$ be the availability, validity and visibility of feature $j$, where $o_{ij}=0$ marks an entry that is unavailable, invalid or hidden during training (§ 2.4.1). An entry $f_{ij}$ is readable when $r_{ij}=a_{ij}v_{ij}o_{ij}=1$. The features other than the position form 21 input groups of closely related features, such as the three components of the velocity or a property and its positivity indicator (Appendix § 8.§ 1). For group $g$, the vector $\boldsymbol{s}_{ig}=|g|^{-1}\sum_{j\in g}(a_{ij},v_{ij},o_{ij},r_{ij})$ holds the fractions of its features that are available, valid, visible and readable, and each subhalo is first embedded as

where $W_g \in\mathbb{R}^{256\times|g|}$ projects the values of group $g$ without bias, $S_g \in\mathbb{R}^{256\times4}$ embeds $s_{ig}$, $b_g$ is a learned bias and $\boldsymbol{e}_0$ a learned base token. Unreadable values are set to zero before the projection, so a hidden or missing value reaches neither the token nor its gradient. The token $\boldsymbol{t}_i \in\mathbb{R}^{256}$ is obtained from $\boldsymbol{u}_i$ by layer normalization followed by a residual two-layer multilayer perceptron (MLP) with GELU activation.

Two modules add spatial context at different resolutions. The grid module supplies aggregated context, the environment of each subhalo on the scale of a grid cell ($1.6\,h^{-1}\,\mathrm{Mpc}$) and above. The attention module supplies pairwise context from the exact separations between subhalos, which resolve satellites and neighbors below this scale. In the first, a linear map projects each layer-normalized token to 32 channels $h_i$, which are deposited with cloud-in-cell weights $w_{ik}$ onto the cells $k$ of a periodic $16^3$ grid. Each cell holds the weighted mean of the deposited channels and an occupancy term,

Three-dimensional convolutions with circular padding process the grid, and trilinear interpolation with the same weights $w_{ik}$ returns the processed grid to the subhalos, where it is added to the tokens.

The second module is one layer of dense attention between all pairs of subhalos, biased by their separation. For fractional positions $\boldsymbol{x}_i$,

so that $\rho_{ij}$ is the minimum-image separation in units of $L$. For attention head $h$, the query-key score is

where $L_* \equiv L/(1\,h^{-1}\,\mathrm{Mpc})=25$ and $n_p=256$, so that the arguments express the separation in units of the box, of $1\,h^{-1}\,\mathrm{Mpc}$ and of the initial mean particle spacing. A small MLP shared by the heads returns the biases $B_h$. Each spatial module is followed by a residual feed-forward block (Appendix § 8.§ 2).

Four induced set-attention blocks (Lee et al. (2019)) then contextualize the tokens. The tokens of the $N_v\leq1,024$ subhalos of a view, after both spatial modules, form the rows of the state matrix $H\in\mathbb{R}^{N_v\times256}$, and each block updates it as

where $I\in\mathbb{R}^{64\times256}$ holds 64 learned inducing vectors and $\operatorname{MAB}(X,Y)$ is a multi-head attention block in which the rows of $X$ attend to those of $Y$, so that its output has the shape of $X$. The inducing vectors gather information from all subhalos of the view and return it to each of them. A final layer normalization of the rows of $H$ gives the contextual states $\boldsymbol{c}_i\in\mathbb{R}^{256}$.

Self-supervised training

The encoder is trained as a student that sees a masked view of a catalog (Figure 1b). The student reconstructs the hidden properties of each subhalo and predicts the representations that a teacher network assigns to the partially hidden subhalos in the complete view. Predicting the teacher’s representations of hidden content follows the joint-embedding predictive architecture of Assran et al. (2023). The target state is given by a slowly updated teacher to keep the representation from collapsing to a constant.

Views and masking

For every training catalog, we draw 64 views of up to 1,024 subhalos from those above the mass threshold, and successive epochs cycle through them. Each view is rotated by one of the 24 rotations of the cube (the combination of 90° rotations around three box axes), translated by a uniform periodic shift and permuted. The rotations act on the positions and on all vector fields. The student input, the teacher input and the reconstruction targets all use the same transformed view.

The student input then loses part of its fields randomly via grouped feature masking. Four fields are never drawn for masking. The positions, the host-relative positions and the central flag describe the set and its host structure, and the host dark-matter fraction is nearly fixed by the other host fractions. The positions are never hidden, and the other three fields are hidden only with their whole host in step 3 or by the rules of step 4. All other fields are masked in families. A family consists of an input group and the groups that describe the same component or are computed from it. The velocity and the velocity relative to the host form one family, and the black-hole mass and the black-hole accretion rate form another. The three host masses $M_{\mathrm{FoF}}$, $M_{200\mathrm{c}}$ and $M_{500\mathrm{c}}$ also form one family. Of the 11 families, the four baryonic ones are built around the stellar mass, the gas mass, the star-formation rate and the black-hole mass (Appendix § 8.§ 1). The mask of a view is drawn in four steps.

1. With probability 0.25, the view of a hydrodynamical catalog loses all four baryonic families. The student then sees an N-body-like description of the catalog, and the hidden baryonic fields remain reconstruction targets.

2. Each remaining family is hidden independently for all subhalos of the view, with probability 0.5 for a non-baryonic family and 1/3 for a baryonic family. Together with step 1, every family is hidden with probability 0.5, since $0.25 + 0.75 \times1/3 = 0.5$.

3. Each host group is chosen with probability 0.2, and all its members in the view lose all their fields except the positions. Hiding whole hosts keeps the host fields of a hidden subhalo from being read from the other members of its host.

4. Two rules hide fields that could be computed from visible ones. The host dark-matter fraction is hidden whenever the host gas, stellar or black-hole fraction is hidden. Whenever the host masses are hidden, the host-relative positions and velocities, the host center offset and the central flag are hidden for all members of the host, since two members would otherwise give $R_{200\mathrm{c}}$ or $V_{200\mathrm{c}}$ and hence $M_{200\mathrm{c}}$ (Appendix § 8.§ 3). This rule hides the host-relative velocities without their family, so the velocities stay visible.

The subhalo total mass equals the sum of the dark-matter, gas, stellar and black-hole masses. It is not an input, so that the baryonic mass stays hidden whenever the baryonic families are hidden. In the student input, a hidden entry has its availability, validity and visibility flags set to zero. A field missing from the catalog carries the same flags, so the student cannot tell which entries were hidden. Teacher input and reconstruction target selection use the original flags, and the teacher sees every available entry of the same view.

Objectives and optimization

The reconstruction head is a two-layer MLP with GELU activation. It maps each student state $c_i$ to predictions $\hat{f}_i$ of the continuous features and to logits for the positivity indicators. The targets are all available and valid entries hidden from the student, except the host-relative positions and the central flag. They include the baryonic features hidden in step 1, so the encoder learns to predict the baryonic properties of a subhalo from its dark-matter properties and its environment. Binary targets that a visible indicator already determines are not scored (Appendix § 8.§ 4).

A continuous target is scored in standardized units with the smooth-$L_1$ penalty of its residual $\Delta= \hat{f}_{bij} - f_{bij}$, where $j$ indexes the feature, $i$ the subhalo and $b$ the view in the batch. The penalty is

The 13 positivity indicators are scored with the binary cross-entropy $\ell_{\mathrm{BCE}}$ between the predicted logit and the original indicator value of 0 or 1. With $\ell_j$ the penalty of feature $j$, the loss averages first over the scored entries $\mathcal{I}_j$ of feature $j$ in the batch, then over the scored features $\mathcal{J}_g$ of a loss group $g$, and finally over the loss groups $\mathcal{G}_+$ with scored entries,

The targets fall into 20 loss groups of physically related features, so that groups with many features do not dominate the loss.

The teacher is a copy of the encoder whose weights follow an exponential moving average (EMA) of the student weights, $\theta_{\mathrm{T}} \leftarrow m\theta_{\mathrm{T}} + (1-m)\theta_{\mathrm{S}}$, and receive no gradient ((Grill et al. 2020; Assran et al. 2023)). It encodes the complete view and returns the states $c_i^{\mathrm{T}}$. The student predicts these states from a masked view in two attention steps. First, 32 learned query vectors attend to embeddings of the student states $c_i$ and then to one another. These catalog slots summarize the masked view in 32 vectors. Second, each subhalo $i \in\mathcal{R}$ reads its prediction from the slots, where $\mathcal{R}$ holds the subhalos that keep at least one visible group other than the position and have at least one hidden group. The query is the token $t_i$ of the subhalo, which is built from its visible fields. It attends with eight heads to the slots, and a final layer normalization of the result gives the prediction $\hat{c}_i$. The loss is the mean cosine distance

with the teacher states held fixed. Subhalos that keep nothing but their positions enter only the reconstruction loss. The total loss is

This value is the ratio of the gradient norms of $\mathcal{L}_{\mathrm{rec}}$ and $\mathcal{L}_{\mathrm{lat}}$ with respect to the encoder weights at the start of training, so that the two terms contribute gradients of similar size.

Each batch holds 32 training catalogs of one hydrodynamical suite together with their N-body counterparts, and successive batches cycle through the three suites. We optimize the student with AdamW (Loshchilov & Hutter 2017) for 18,000 steps with a warmup–stable–decay learning-rate schedule, and the EMA momentum $m$ rises from 0.996 to 1 during the decay (Appendix § 7). All results use the student weights at the end of this single training run.

Catalog representation

After training, we freeze the student model and make every available field visible (Figure 1c). To encode a full catalog, we partition its $N$ subhalos above the mass threshold by a fixed random permutation into blocks of at most 1,024 and encode each block as one view. A last block with fewer than 1,024 subhalos is filled up with other subhalos of the catalog, which serve as context, and only the outputs of its own subhalos are kept. Each subhalo thus contributes one state $\mathbf{c}_i$.

We summarize the catalog by the centered covariance of its states,

As an average over the subhalos, $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ do not depend on the order of the subhalos and are normalized by their count. The symmetric vectorization $\boldsymbol{\xi}=\operatorname{svec}(\boldsymbol{\Sigma})\in\mathbb{R}^{32,896}$ keeps the diagonal and multiplies the upper off-diagonal entries by $\sqrt{2}$, so that it preserves the Frobenius inner product of covariance matrices.

Every readout is fitted on a set $\mathcal{T}$ of labeled catalogs. On these catalogs alone, we center the vectors $\boldsymbol{\xi}$ and compute their principal components, keeping all those with eigenvalues above $10^{-12}$ of the largest. With the mean $\bar{\boldsymbol{\xi}}_{\mathcal{T}}$ and the matrix $P_{\mathcal{T}}$ of retained components, the catalog embedding is

where the scalar $s_{\mathcal{T}}$ sets the mean variance of the retained components over $\mathcal{T}$ to one, so that the regularization strength below is in units of this variance (Appendix C).

Linear readouts and error bars

The readout targets are

For each target $y$, the readout minimizes

over the weights $\mathbf{w}$ and the intercept $b$. A local readout is fitted on the 700 training catalogs of one suite, and a pooled readout on the 2,100 training catalogs of IllustrisTNG, SIMBA and Astrid together, with its predictions evaluated separately for each suite. The regularization strength $\alpha$ is chosen by a three-fold cross-validation on the training catalogs (Appendix C). Over the test catalogs of a suite, we report the coefficient of determination $R^2=1-\sum_n(y_n-\hat{y}_n)^2/\sum_n(y_n-\bar{y})^2$, the root-mean-square error and the Pearson correlation $r$.

Every prediction $\hat{y}(\mathbf{z})$ comes with an error bar $\hat{y}(\mathbf{z})\pm q_\beta\hat{\sigma}(\mathbf{z})$ for a nominal coverage $\beta$ and with a confidence score $p_{\mathrm{conf}}$. Both require a model $\hat{\sigma}(\mathbf{z})$ of the error size and multipliers $q_\beta$ that calibrate the error bars to the coverage $\beta$. Like the readout, they are fitted on the training catalogs and then applied unchanged to the test catalogs.

The size of the readout error varies from catalog to catalog. A catalog with few massive black holes, for example, constrains the AGN feedback less than one with many. The error model therefore predicts the error size of each catalog from its embedding.

The error model is fitted to the out-of-fold residuals $r_n=y_n-\hat{y}^{\mathrm{oof}}_n$, in which every training catalog is predicted by a readout that did not see its label, so that the residuals represent the errors on new catalogs. We write the error as $r=\sigma(\mathbf{z})\varepsilon$, with $\varepsilon$ independent of $\mathbf{z}$. The variable

satisfies $\mathbb{E}[Z\mid z]\simeq\ln\sigma^2(z)+\kappa$, with the constant $\kappa=\mathbb{E}[\ln\varepsilon^2]$ (Harvey 1976). Here $V_0$ is the variance of $y$ over the training catalogs, and the floor $\eta$ avoids $\ln0$. A ridge regression of $Z_n$ on $z_n$ gives $\hat{Z}(z)$ and the error size $\hat{\sigma}(z)=\exp[\hat{Z}(z)/2]$, which is proportional to $\sigma(z)$. Its regularization strength is chosen on the same grid and folds as for the readout.

The multipliers are computed out of fold as well. For each fold, an error model fitted to the residuals of a cross-validation within the other two folds predicts $\hat{\sigma}_n$ for the catalogs of the held-out fold. With the normalized residuals $R_n=|r_n|/\hat{\sigma}^{\mathrm{oof}}_n$ of the $N_{\mathcal{T}}$ catalogs in $\mathcal{T}$ sorted as $R_{(1)}\leq\cdots\leq R_{(N_{\mathcal{T}})}$,

This quantile follows split-conformal and cross-validation-based prediction intervals (Lei et al. 2018; Barber et al. 2021). Applied to the readout fitted on all training catalogs, it calibrates the coverage empirically, and Appendix D tests this calibration on the test catalogs. The unknown factor $e^{\kappa/2}$ in $\hat{\sigma}$ cancels in the product $q_\beta\hat{\sigma}$, so the distribution of $\varepsilon$ need not be known.

The confidence score compares the 68% error bar with the spread of the parameter over the training catalogs,

For Gaussian errors, the 68% half-width is close to $\sigma$, so $p_{\mathrm{conf}}\simeq1-\sigma^2(z)/V_0$ is a per-catalog analogue of $R^2$, and its mean over a set of catalogs approximates $1-\mathrm{MSE}/V_0$.

Results

We now examine what the catalog embedding $\mathbf{z}$ contains. Unless stated otherwise, all results use the frozen encoder and the 200 test catalogs of each suite. We first visualize the latent space (§ 3.1) and measure the parameter information that a linear map extracts from it, with all input features (§ 3.2) and with the restricted inputs of earlier supervised studies (§ 3.3). We then apply the encoder to a simulation code that it never saw (§ 3.4) and relate the geometry of the embedding to the physical effect of each parameter (§ 3.5). Section 3.6 connects the confidence of the readout with the black-hole population of each catalog, and Section 3.7 follows the numbers of massive galaxies and black holes along the principal axes of the embedding.

Visualization of the latent space

Figure 2 projects the embeddings of the 600 test catalogs onto the first two principal components of the 2,100 training catalogs of the three suites, which carry 37.4% and 16.9% of the total variance. The suites occupy separate regions of this plane. The first component separates SIMBA from IllustrisTNG and Astrid, and the second separates Astrid from IllustrisTNG. Note that the encoder never receives the identity of the code during training. It learns this separation automatically from the catalogs themselves, via the intrinsic differences between how the three galaxy-formation models populate their halos with galaxies, gas and black holes ([70] (Villaescusa-Navarro et al. 2021)).

Within each suite, the catalogs are ordered by $\Omega_{\mathrm{m}}$, and in IllustrisTNG and Astrid also by a supernova-feedback amplitude. The IllustrisTNG catalogs form an elongated sequence along which $\Omega_{\mathrm{m}}$ decreases from the compact lower-left end to the extended upper-right end, and $A_{\mathrm{SN1}}$ changes across the sequence. In SIMBA and Astrid, $\Omega_{\mathrm{m}}$ changes smoothly across the point cloud with the same trend. In Astrid, $A_{\mathrm{SN2}}$ decreases from the lower to the upper end of the cloud, and $\sigma_8$ varies across it more weakly. The AGN amplitudes show no clear trend in this plane. A linear readout of the full embedding recovers $A_{\mathrm{AGN1}}$ in SIMBA with $R^2 = 0.94$ and $A_{\mathrm{AGN2}}$ in Astrid with $R^2 = 0.996$ (§ 3.2), so the information on AGN feedback lies along directions of lower variance.

Linear readout of the simulation parameters

For each suite and parameter, we fit a ridge regression on the embeddings of the 700 training catalogs, and apply it to the 200 test catalogs. Figure 3 compares the predictions with the true values, and Table 2 lists the test $R^2$ and root-mean-square error (RMSE). A linear map recovers the cosmological and supernova-feedback parameters with $R^2 \ge0.97$ in all three suites, and the AGN parameters with a precision that depends on the code.

Suite	Parameter	$R^2$	RMSE	$\langle p_{\mathrm{conf}}\rangle$
IllustrisTNG	$\Omega_{\mathrm{m}}$	0.999	0.003	0.999
	$\sigma_8$	0.986	0.014	0.988
	$A_{\mathrm{SN1}}$	0.999	0.011	0.999
	$A_{\mathrm{SN2}}$	0.991	0.016	0.992
	$A_{\mathrm{AGN1}}$	0.762	0.176	0.751
	$A_{\mathrm{AGN2}}$	0.627	0.105	0.727
SIMBA	$\Omega_{\mathrm{m}}$	0.999	0.003	0.999
	$\sigma_8$	0.991	0.010	0.995
	$A_{\mathrm{SN1}}$	0.999	0.008	0.999
	$A_{\mathrm{SN2}}$	0.975	0.028	0.967
	$A_{\mathrm{AGN1}}$	0.945	0.084	0.953
	$A_{\mathrm{AGN2}}$	0.836	0.070	0.909
Astrid	$\Omega_{\mathrm{m}}$	0.999	0.004	0.999
	$\sigma_8$	0.990	0.012	0.992
	$A_{\mathrm{SN1}}$	0.996	0.021	0.999
	$A_{\mathrm{SN2}}$	0.997	0.009	0.999
	$A_{\mathrm{AGN1}}$	0.615	0.230	0.591
	$A_{\mathrm{AGN2}}$	0.996	0.021	0.997

Cosmological parameters. — The matter density is recovered almost exactly, with $R^2 = 0.999$ in all three suites and an RMSE of 0.003 to 0.004, below 1% of the prior range

$0.1 \le\Omega_{\mathrm{m}} \le0.5$. This $\Omega_{\mathrm{m}}$ information is likely carried by two components: the standard large-scale clustering signal in the galaxy positions ((Hahn et al. 2023)) as well as the $\Omega_{\mathrm{m}}$-modulated correlations in individual galaxy properties. The latter have previously been shown empirically in the CAMELS simulations ((Villaescusa-Navarro et al. 2022c; Hahn et al. 2024)) and explained via the mechanism that galaxy properties constrain the universal baryon fraction $\Omega_{\mathrm{b}/\Omega_{\mathrm{m}}}$, which when $\Omega_{\mathrm{b}}$ is held constant in CAMELS, provides precise constraints on $\Omega_{\mathrm{m}}$ ((Lue et al. 2025; Liao et al. 2026)). For $\sigma_8$, $R^2 = 0.986$–$0.991$ and the RMSE is $0.010$–$0.014$. The readouts reach nearly the same precision in the three codes, which form galaxies differently. A readout fitted on Swift-EAGLE labels reaches the same precision in a code that the encoder never saw (§ 3.4). The $\sigma_8$ constraints are similar to those of $\Omega_{\mathrm{m}}$, combining information from large-scale clustering as well as ensemble statistics of galaxy populations, driven by differences in merger histories and environments ((Chawak et al. 2024; Lovell et al. 2025)).

Supernova feedback. — The two supernova-feedback amplitudes are recovered with $R^2 \ge0.97$ in every suite, $A_{\mathrm{SN1}}$ with $R^2 = 0.996$–$0.999$ and $A_{\mathrm{SN2}}$ with $0.975$–$0.997$. Galactic winds regulate the stellar mass, gas content and metallicity of nearly every galaxy, so their strength leaves an imprint on most of the objects that the encoder sees. In the one-parameter variations, raising a supernova parameter to its largest value also moves the catalogs furthest in every suite. (§ 3.5).

AGN feedback. — The precision for the AGN amplitudes depends on the code and on the feedback mode that they scale. In SIMBA, $A_{\mathrm{AGN1}}$ and $A_{\mathrm{AGN2}}$ reach $R^2 = 0.94$ and $0.84$, and in Astrid the thermal-mode amplitude $A_{\mathrm{AGN2}}$ reaches $0.996$. The three parameters with the lowest $R^2$ are $A_{\mathrm{AGN1}}$ and $A_{\mathrm{AGN2}}$ in IllustrisTNG ($0.76$ and $0.63$) and $A_{\mathrm{AGN1}}$ in Astrid ($0.62$), and all three scale the kinetic feedback mode. This mode operates mainly in massive black holes. In IllustrisTNG, the Eddington ratio below which a black hole switches to kinetic feedback grows with black-hole mass about a pivot of $10^8\,M_\odot$ ((Weinberger et al. 2017)), and black holes in Astrid enter the kinetic mode above $5 \times10^8\,h^{-1}\,M_\odot$ ([70]). Only a small fraction of the subhalos in a $25\,h^{-1}\,\mathrm{Mpc}$ box host such black holes, so these amplitudes act directly on few of the objects in a catalog. § 3.6 shows that the information on them grows with the number of black holes above these masses.

Confidence score. — Each prediction carries an error bar and a confidence score $p_{\mathrm{conf}}$, both computed from the embedding alone (§ 2.6). The score is one minus the squared ratio of the half-width of the 68% error bar to the standard deviation of the parameter over the training catalogs (equation (19)). Like $R^2$, it approaches 1 when the error bar is small compared with the spread of the parameter and 0 when it is as wide as that spread.

The shading in Figure 3 shows $p_{\mathrm{conf}}$ for every test catalog, and Table 2 lists its mean. It is close to 1 for the cosmological and supernova parameters and for $A_{\mathrm{AGN2}}$ in Astrid. For the other AGN parameters, it varies from catalog to catalog, and the catalogs with low $p_{\mathrm{conf}}$ scatter most widely about the diagonal. In SIMBA, the five catalogs with the largest errors in each AGN parameter all lie in the quarter of catalogs with the lowest $p_{\mathrm{conf}}$. § 3.6 relates $p_{\mathrm{conf}}$ to the accuracy of the readout and to the black-hole population of each catalog, and Appendix § 10 tests the calibration of the error bars.

Restricted inputs and comparison with supervised networks

The frozen encoder also embeds catalogs in which only some features are available, because it was trained to represent catalogs with randomly hidden inputs. We use this property to compare the embedding with the supervised graph neural networks (GNNs) of two earlier studies, which infer $\Omega_{\mathrm{m}}$ from restricted galaxy catalogs. For each input set, we select galaxies by their stellar content as described below, without the dark-matter threshold of § 2.2, hide all other features, embed the training and test catalogs with the unchanged encoder, and fit a new linear readout on the training catalogs. De Santi et al. [109] trained GNNs on the 3D positions and line-of-sight velocities of galaxies with $M_\star> 1.95 \times10^8\,h^{-1}\,M_\odot$ and tested them on catalogs of the same suite. With the same input, the linear readout of our embedding reaches $R^2 = 0.925$, $0.958$ and $0.905$ in IllustrisTNG, SIMBA and Astrid (Figure 4), compared with $0.911$, $0.928$ and $0.835$ for their GNNs.

Villanueva-Domingo & Villaescusa-Navarro (2022) used positions alone, positions with stellar mass, stellar half-mass radius and stellar metallicity, and these with $V_{\max}$ added, in IllustrisTNG and SIMBA. We give the encoder the same three input sets for all subhalos with at least 20 star particles (Figure 5). With positions alone, their GNNs reach $R^2 = 0.83$ and 0.85 and our readout 0.71 and 0.81. With the stellar properties, the embedding reaches $R^2 = 0.96$ and 0.97, comparable to the 0.96 and 0.94 of their GNNs, and with $V_{\max}$ added the embedding and the GNNs both reach $R^2 = 0.98$ in each suite.

Once the positions are complemented by the line-of-sight velocity (Figure 4) or by the stellar properties (Figure 5), the linear readout of a single encoder trained without labels reaches a precision comparable to that of the GNNs, which were each trained for one parameter and one input set. § 4.2 discusses the case of positions alone. The precision grows with every property that becomes visible. In IllustrisTNG, the RMSE of $\Omega_{\mathrm{m}}$ falls from 0.064 with positions alone to 0.024 with the stellar properties and 0.016 with $V_{\max}$ added, and in SIMBA from 0.048 to 0.019 and 0.015. With all input features of every subhalo, the full catalogs give 0.003 in both suites, and the same embedding carries information on all six parameters.

Transfer to an unseen simulation code

The Swift-EAGLE suite combines the SWIFT code with a galaxy-formation model based on EAGLE (§ 2.1), and its catalogs were used in neither the self-supervised training nor the input normalization. We embed them with the frozen encoder and read them out in two ways (Figure 6).

Zero-shot readout. — We apply the pooled readout of the training suites to the Swift-EAGLE catalogs. The $\Omega_{\mathrm{m}}$ predictions correlate with the true values with $r = 0.985$. They follow a line of slope 0.77 and lie above the true values by 0.04 on average, which gives $R^2 = 0.825$. The $\sigma_8$ predictions lie above the true values by 0.32 on average, with $r = 0.79$. The $A_{\mathrm{SN1}}$ predictions follow the true values with a slope of 1.00 and $r = 0.89$ and lie below them by 0.25 in $\log_{10} A_{\mathrm{SN1}}$, which gives $R^2 = 0.20$. The embedding of the new code thus varies with $\Omega_{\mathrm{m}}$, $\sigma_8$ and $A_{\mathrm{SN1}}$ along directions that the readout of the training suites recognizes, and each of these predictions needs a new calibration. The training suites share their response directions for the same three parameters (§ 3.5). The other feedback parameters do not transfer ($r \le0.28$). In Swift-EAGLE, $A_{\mathrm{SN2}}$, $A_{\mathrm{AGN1}}$ and $A_{\mathrm{AGN2}}$ scale the metallicity dependence of the stellar feedback energy, the Bondi accretion rate and the heating temperature of AGN feedback, and in the training codes the same amplitudes scale other processes (Table 1).

Readout with Swift-EAGLE labels. — With labels, the embedding of Swift-EAGLE is as informative as those of the training suites. A linear readout fitted on its 700 training catalogs reaches $R^2 = 0.999$, 0.992, 0.998 and 0.979 for $\Omega_{\mathrm{m}}$, $\sigma_8$, $A_{\mathrm{SN1}}$ and $A_{\mathrm{SN2}}$ on its 200 test catalogs. For $A_{\mathrm{AGN1}}$ and $A_{\mathrm{AGN2}}$, it reaches 0.988 and 0.887, higher than in IllustrisTNG and SIMBA.

Label efficiency. — Figure 7 shows how many labeled catalogs these readouts need. For each $n < 700$, we draw 100 random sets of $n$ training catalogs, choose $\alpha$ and fit the readout on each set, and evaluate each fit on the test catalogs. Besides $R^2$, we show the Pearson correlation $r$ between predicted and true values, which measures how well the predictions follow the true values up to a linear calibration. For Swift-EAGLE, 20 labeled catalogs give $r \ge0.94$ for $\Omega_{\mathrm{m}}$, $\sigma_8$ and $A_{\mathrm{SN1}}$, and 100 give $r \ge0.96$ for all parameters except $A_{\mathrm{AGN2}}$ ($r = 0.85$). With 20 labels, $\Omega_{\mathrm{m}}$ reaches $R^2 \ge0.94$ and $\sigma_8$ reaches $R^2 \ge0.87$ in all four codes. Of the cosmological and supernova parameters, $A_{\mathrm{SN2}}$ needs the most additional labels in Swift-EAGLE, with $R^2 = 0.72$ for 20 labels against 0.77 to 0.95 in the training suites. The Swift-EAGLE curves fall within or close to the range of the training suites, so the new code needs about as many labels as the codes on which the encoder was trained.

The number of labels depends on how strongly a parameter imprints on the catalogs. With 100 labels, $A_{\mathrm{AGN1}}$ reaches $R^2 = 0.93$ in Swift-EAGLE, 0.88 in SIMBA, 0.28 in IllustrisTNG and about zero in Astrid. For the three kinetic-mode parameters, the readout finds no signal with 50 or fewer labeled catalogs ($|r| \le0.07$). Their $r$ rises steeply between 50 and 200 catalogs, and their $R^2$ is still increasing at $n = 700$.

For a new code, the cosmological directions of the embedding carry over without labels, and about 100 labeled catalogs recover its feedback parameters with $r \gtrsim0.85$.

Response to parameters and to initial conditions

The geometry of the embedding reflects the physical effect of each parameter. The 1P and CV sets of CAMELS probe this geometry around the fiducial model. Here we compute $z$ from equation 14 with $\mathcal{T}$ set to the 2,100 training catalogs of the three suites and $\sigma_{\mathcal{T}} = 1$, and keep its 2,048 leading components. For each parameter $a$, we take the displacement $\boldsymbol{d}_a = \boldsymbol{z}(\theta_a^{\max}) - \boldsymbol{z}_0$ from the fiducial 1P catalog to the one with the largest value of $a$, and project the 1P and CV catalogs onto $\boldsymbol{u}_a = \boldsymbol{d}_a/\lVert\boldsymbol{d}_a\rVert$ (Figure 8). The length $\lVert\boldsymbol{d}_a\rVert$ measures the response distance, and the standard deviation $\sigma_a$ of the CV catalogs along $\boldsymbol{u}_a$ measures the scatter from initial conditions.

Response and cosmic variance. — Apart from the three kinetic-mode amplitudes discussed below, every parameter mover a catalog well beyond this scatter, and its five 1P catalogs follow the parameter value monotonically along $\boldsymbol{u}_a$, with a single exception at the two lowest values of $\Omega_{\mathrm{m}}$ in SIMBA. The supernova parameters shift the catalog by 13 to 59 times $\sigma_a$, $\Omega_{\mathrm{m}}$ by 9 to 19 times and $\sigma_8$ by 4 to 7 times. The AGN parameters do so in SIMBA (20 and 4.6 times), and $A_{\mathrm{AGN2}}$ does so in Astrid (39 times). For $A_{\mathrm{AGN1}}$ and $A_{\mathrm{AGN2}}$ in IllustrisTNG and for $A_{\mathrm{AGN1}}$ in Astrid, the four lower values leave the catalog within about $\sigma_a$ of the fiducial catalog, and only the largest value moves it, by $\lVert\boldsymbol{d}_a\rVert= 1.7$, 2.0 and 1.9, or 4.0, 3.2 and 1.6 times $\sigma_a$. These displacements are at least four times shorter than those of all other parameters. The 1P response alone thus singles out the three parameters with the lowest readout $R^2$.

Physical content of the response. — Figure 9 shows how the catalogs change along these displacements. For one parameter in one suite per panel, we interpolate 11 equally spaced points on the straight line in the embedding from the 1P catalog with the lowest value of the parameter to the one with the highest value, and represent each point by the five nearest of all 1,000 LH catalogs of the same suite. Because the retrieved catalogs are simulations, their statistics show what each point of the line represents. Along every line, the statistics shift from one end to the other in the direction expected for the parameter.

Increasing $\Omega_{\mathrm{m}}$ from 0.1 to 0.5 raises the abundance of subhalos in IllustrisTNG with dark-matter masses of $10^{11}$–$10^{12}\ h^{-1}\,M_{\odot}$ by factors of 8 to 9. Increasing $\sigma_8$ from 0.6 to 1.0 raises the abundance of SIMBA subhalos by a factor of 1.4 at $V_{\max} \approx110\ \mathrm{km}\,\mathrm{s}^{-1}$ and of 12 at 350 $\mathrm{km}\,\mathrm{s}^{-1}$. The wind energy $A_{\mathrm{SN1}}$ lowers the median ratio of stellar to halo mass of IllustrisTNG central galaxies by a factor of about 40 at $M_{200\mathrm{c}} = 10^{11}$–$10^{11.5}\ h^{-1}\,M_{\odot}$ and of 9 at $10^{12}$–$10^{12.25}\ h^{-1}\,M_{\odot}$, and the wind speed $A_{\mathrm{SN2}}$ lowers the abundance of Astrid galaxies with $M_\star= 10^{9.2}$–$10^{9.8}\ h^{-1}\,M_{\odot}$ by factors of 20 to 28. The AGN parameters act on the gas and on the black holes. $A_{\mathrm{AGN1}}$ lowers the median baryon fraction of SIMBA FoF halos at all masses shown, by 10% to 60% and least near $M_{\mathrm{FoF}} = 10^{12}\ h^{-1}\,M_{\odot}$. The thermal-mode amplitude $A_{\mathrm{AGN2}}$ lowers the median black-hole mass of Astrid galaxies at fixed stellar mass by a factor of about 7 at $M_\star= 10^{8}$–$10^{8.5}\ h^{-1}\,M_{\odot}$ and by factors of 14 to 18 at $10^{9.75}$–$10^{10.75}\ h^{-1}\,M_{\odot}$. We choose these statistics for their strong response to the parameters, and the same retrieval applies to any statistic that the catalogs provide. A systematic study of how observables respond to each parameter is left for future work.

Directions of the response. — Figure 10 compares the directions in which the parameters move the catalogs across the Latin hypercube. For each suite, we regress the covariance vectors $\boldsymbol{\xi}$ of all 1,000 LH catalogs jointly on the standardized targets of equation (15), and the coefficient vector of each parameter is its response direction. Noise in the estimated directions biases their norms upwards, and it also biases the inner products between directions of the same suite. We therefore fit every direction on 200 random splits of the catalogs of its suite into two halves, and take the norms and the within-suite inner products from the mean product of the two half-sample estimates.

The cosmological directions are shared across codes. The $\Omega_{\mathrm{m}}$ directions of the three suites have cosines of 0.44 to 0.67

with one another and the $\sigma_8$ directions 0.52 to 0.55, and within each code the two are partly aligned (0.48 to 0.66). This is expected, because gravity and initial conditions act in the same way in all codes and both parameters raise the abundance of massive halos. For the astrophysical parameters, we compare the directions with the responses of galaxies, black holes and gas in the 1P simulations.

Within a code, the alignments agree with these responses. In SIMBA, the momentum flux $A_{\mathrm{AGN1}}$ and the jet speed $A_{\mathrm{AGN2}}$ align (0.52), and the wind speed $A_{\mathrm{SN2}}$ points against the jet speed ($-0.63$). Both AGN parameters spread the gas of the SIMBA 1P set farther from the dark matter next to which it started, and faster winds reduce this spread (Fig. 4 of Gebhardt et al. 2024), because stellar feedback limits the growth of black holes and thereby the output of the jets (Fig. 4 of Delgado et al. 2023; Gebhardt et al. 2024). In IllustrisTNG, the energy $A_{\mathrm{AGN1}}$ and the ejection speed $A_{\mathrm{AGN2}}$ of the kinetic mode align (0.60), and both lower the baryon fractions of halos (Fig. 2 of Delgado Ni et al. [70]. In Astrid, the wind energy $A_{\mathrm{SN1}}$ and the wind speed $A_{\mathrm{SN2}}$ align (0.64), and both reduce the abundance of galaxies (Figure 5 of Ni et al. [70]).

Across codes, parameters that scale similar processes align. The $A_{\mathrm{SN1}}$ directions have cosines of 0.33 to 0.50, and a larger $A_{\mathrm{SN1}}$ reduces the abundance of galaxies at all masses in all three codes (Figure 5 of Ni et al. [70]). The largest of these cosines joins IllustrisTNG and Astrid, where $A_{\mathrm{SN1}}$ scales the same wind energy per unit star-formation rate. The ejection speed of the kinetic mode in IllustrisTNG aligns with the jet speed of SIMBA (0.46), and both lower the baryon fractions of massive halos (Figure 2 of Delgado et al. 2023). The kinetic-mode energy $A_{\mathrm{AGN1}}$ of Astrid aligns with the same ejection speed (0.57). Astrid and IllustrisTNG both inject kinetic feedback locally and without collimation around massive black holes at low accretion rates (Ni et al. [70]), and in both codes the confidence of the readout follows the number of black holes massive enough to drive this mode (§ 3.6).

The thermal-mode amplitude $A_{\mathrm{AGN2}}$ of Astrid points against the kinetic directions of Astrid and IllustrisTNG ($-0.53$, $-0.31$ and $-0.43$) and against the jet speed of SIMBA ($-0.23$). A larger thermal-mode amplitude limits the formation of the massive black holes that turn on the kinetic mode and raises the star-formation rate and the abundance of galaxies in Astrid (Figure 4 and Figure 5 of Ni et al. [70]). This is expected from self-regulated black-hole growth, in which more efficient feedback leaves black holes less massive (Booth & Schaye 2009). § 3.7 follows this effect to the numbers of massive galaxies and black holes.

With responses of 47 and 44 times the scatter from initial conditions, the aligned supernova parameters of Astrid are both recovered with $R^2 > 0.99$. The kinetic-mode parameters of IllustrisTNG move the catalog by about a tenth of the displacement for $\Omega_{\mathrm{m}}$. The direction of $A_{\mathrm{AGN1}}$ also overlaps with those of $\Omega_{\mathrm{m}}$ (0.56 and 0.54) and with the displacements between CV catalogs (0.31), so changes of cosmology or of the initial conditions can mimic it.

Initial conditions. — The bottom row of Figure 10 compares each response direction with the displacements between CV catalogs, which are the directions along which a change of initial conditions moves a catalog. The root-mean-square cosine is 0.32 to 0.47 for $\sigma_8$ and 0.22 to 0.37 for $\Omega_{\mathrm{m}}$, and it is at most 0.14 for eight of the twelve feedback directions. In a 25 $h^{-1}\,\mathrm{Mpc}$ volume, the number of rare massive halos varies from one realization to the next, and this cosmic variance mimics a change of cosmology. Among the feedback parameters, $A_{\mathrm{SN2}}$ and $A_{\mathrm{AGN2}}$ in SIMBA and $A_{\mathrm{AGN1}}$ in IllustrisTNG overlap with the initial conditions about as much as $\sigma_8$ does in SIMBA (0.34, 0.32 and 0.31). The response to $A_{\mathrm{SN1}}$ points in directions that cosmic variance hardly reaches (at most 0.13). The precision of a readout thus follows from the length of the response to a parameter and from the overlap of its direction with those of the other parameters and of the initial conditions.

Confidence score and black-hole content

The mean confidence over the test catalogs agrees with the accuracy of the readout. It agrees with the test $R^2$ to within 0.03 in 16 of the 18 combinations Table 2. The two exceptions are $A_{\mathrm{AGN2}}$ in IllustrisTNG and in SIMBA, where it is 0.73 and 0.91 against $R^2 = 0.63$ and 0.84. Because $p_{\mathrm{conf}}$ needs only the embedding, it anticipates how well each parameter is measured in a set of catalogs whose true parameters are unknown.

The score also ranks individual catalogs. We call a catalog high, middle or low confidence when $p_{\mathrm{conf}} \ge0.7$, $0.3 \le p_{\mathrm{conf}} < 0.7$ or $p_{\mathrm{conf}} < 0.3$. Figure 11 compares the readout in the three tiers for the three parameters with the lowest overall $R^2$ ($A_{\mathrm{AGN1}}$ and $A_{\mathrm{AGN2}}$ in IllustrisTNG and $A_{\mathrm{AGN1}}$ in Astrid). For these parameters, 72.5%, 69.5% and 49% of the test catalogs have high confidence. From the low to the high tier, $R^2$ rises and the RMSE falls in all three. The high-confidence catalogs are measured well, with $R^2 = 0.808$, 0.788 and 0.858 against 0.762, 0.627 and 0.615 for all test catalogs, and the RMSE drops from 0.176, 0.105 and 0.230 to 0.157, 0.081 and 0.145. In the low-confidence catalogs, the readout carries little information ($R^2 = 0.134$, −0.627 and −0.090).

For the well-measured AGN parameters, 97.5% and 97% of the SIMBA catalogs have high confidence for $A_{\mathrm{AGN1}}$ and $A_{\mathrm{AGN2}}$, and all Astrid catalogs have high confidence for $A_{\mathrm{AGN2}}$. The score thus picks out the catalogs whose AGN parameters are well measured.

The confidence follows the black-hole population of each catalog (Figure 12). Kinetic or jet feedback becomes important in black holes above $10^{8}\ M_{\odot}$ in IllustrisTNG (Weinberger et al. 2017), $4.5 \times10^{7}\ M_{\odot}$ in SIMBA (Ni et al. [70]) and $5 \times10^{8}\ h^{-1}\ M_{\odot}$ in Astrid (Ni et al. [70]), and we count the subhalos whose black-hole mass exceeds these values. The median test catalog contains 103 such hosts in SIMBA, 82 in IllustrisTNG and 18 in Astrid, in the same order as the $R^{2}$ of $A_{\mathrm{AGN1}}$ (0.94, 0.76 and 0.62). In IllustrisTNG, the mean $p_{\mathrm{conf}}$ of $A_{\mathrm{AGN1}}$ and $A_{\mathrm{AGN2}}$ rises from 0.48 and 0.41 in catalogs with about 20 such hosts to 0.89 and 0.90 in catalogs with about 170. In SIMBA, catalogs with about 30 hosts already reach $p_{\mathrm{conf}} = 0.88$ for $A_{\mathrm{AGN1}}$ and 0.79 for $A_{\mathrm{AGN2}}$, and above 100 hosts both exceed 0.9. In Astrid, $p_{\mathrm{conf}}$ for the kinetic-mode amplitude $A_{\mathrm{AGN1}}$ rises from 0.03 in catalogs without such black holes to 0.53 in catalogs with about 10 and 0.84 in catalogs with about 75. The thermal-mode amplitude $A_{\mathrm{AGN2}}$ of Astrid has $p_{\mathrm{conf}} > 0.97$ in every catalog.

The information on kinetic and jet feedback in the embedding grows with the number of black holes massive enough to drive it. The error model is fitted to the readout residuals of the training catalogs and is never given the black-hole counts, and its confidence tracks this number from catalog to catalog.

Statistical effect along the principal axes

The directions studied so far are defined by the parameters. Here we use directions that the embedding defines without labels. For each suite, we compute $\mathbf{z}$ with $T$ the 1,000 LH catalogs of the suite (equation (14)), so that its component $z_k$ is the coordinate of a catalog along the $k$-th principal axis of these catalogs. Along each axis, we take the points one standard deviation on either side of the mean and retrieve the five LH catalogs nearest to each point, as for Figure 9. The arrows in Figure 13 give the change between the two sets of catalogs in $\log_{10}$ of the mean numbers of galaxies with $M_{\star} > 10^{10.5}\ h^{-1}\ M_{\odot}$ and of black holes with $M_{\mathrm{BH}} > 10^{8}\ h^{-1}\ M_{\odot}$. The figure shows the two axes of each suite with the longest arrows. The semi-axes of the ellipses are the standard deviations of the same changes between pairs of CV simulations of the suite, each represented by its five nearest LH catalogs, and measure the scatter from initial conditions.

The parameters enter only afterwards. Over the 1,000 LH catalogs, we fit the coordinate along each axis by least squares with a linear function of the logarithms of the six parameters $\theta_j$,

which describes the axis by a power law of the parameters. Parameter $j$ accounts for the share $f_{kj}=b_{kj}^{2}\operatorname{Var}(\log_{10}\theta_j)/\operatorname{Var}(z_k)$ of the variance of $z_k$, and the leading parameter $\theta_{d_k}$ is the one with the largest share. Its exponent is one, and we orient each axis so that $b_{kd_k}>0$. The labels of Figure 13 list the terms in decreasing order of $f_{kj}$ until they reach 95% of $\sum_j f_{kj}$, with at most three terms, and each arrow points towards increasing power law. For the seven axes shown, the six-parameter fit accounts for at least 54% of the variance of $z_k$.

Along most axes, massive galaxies and black holes change in the same direction. In IllustrisTNG and Astrid, the axes that follow the supernova parameters lower both. The leading axis of Astrid follows mainly the wind speed $A_{\mathrm{SN2}}$ and reduces the number of massive galaxies by a factor of 11 and that of massive black holes by a factor of 2.9. In IllustrisTNG, the second axis follows $A_{\rm SN1}\Omega_{\rm m}^{-0.60}\sigma_8^{-1.8}$ and reduces them by factors of 3.2 and 4.6, and the third axis follows $A_{\rm SN2}\sigma_8^{-1.2}$ and reduces them by factors of 2.5 and 1.3.

In SIMBA, the third axis follows mainly $A_{\rm SN2}$ and reduces the number of massive black holes by a factor of 2.0 and that of massive galaxies by 1.2, in line with the suppression of black-hole growth by faster winds (§ 3.5). The leading axis of SIMBA and the second axis of Astrid follow $\Omega_{\rm m}$ and $\sigma_8$. They raise the numbers of massive black holes by factors of 2.6 and 3.7 and those of massive galaxies by 1.5 and 2.0.

One axis moves the two populations apart. Among the axes of all three suites whose six-parameter fit accounts for at least half of the variance of $z_k$, only the fifth axis of Astrid changes the two numbers in opposite directions, each by more than three times the corresponding semi-axis of the ellipse. Figure 13 shows it as a dashed arrow. This axis follows $A_{\rm AGN2}\sigma_8^{1.5}$, with the thermal-mode amplitude as its leading term. Along it, massive galaxies become 1.7 times more numerous and massive black holes 2.8 times rarer. The same effect appears as the anti-alignment of the thermal and kinetic AGN directions (§ 3.5) and as the lower black-hole masses along $A_{\rm AGN2}$ in Figure 9. In the picture of self-regulated growth of § 3.5, stronger thermal feedback leaves fewer black holes massive enough to enter the kinetic mode, and fewer massive galaxies are then quenched by it.

The principal axes of the embedding follow combinations of parameters that change the numbers of massive galaxies and black holes. Among the axes shown, the number of massive galaxies changes most along the supernova axes in IllustrisTNG and Astrid and along the cosmological axis in SIMBA.

Discussion

This section discusses what the embedding shares across galaxy-formation models (§ 4.1), where its information on the parameters comes from (§ 4.2), how it can be used (§ 4.3) and how it can be extended to observations and to broader simulation suites (§ 4.4).

A shared representation across galaxy-formation models

The embedding keeps the catalogs of the three codes in separate regions (Figure 2), and a change of cosmology moves them in partly aligned directions (Figure 10). One representation holds both the differences between the codes and a shared response to cosmology, and this shared response carries over to Swift-EAGLE, which the encoder never saw (§ 3.4). Supervised networks infer $\Omega_{\rm m}$ robustly across codes from the positions and velocities of galaxies or halos ([109] (Shao et al. 2023)), and internal properties can tie them more closely to the code they were trained on ((Villescusa-Navarro et al. 2022c; Echeverri-Rojas et al. 2023; Shao et al. 2023)). The zero-shot readout of Celia draws on the internal properties as well and still follows $\Omega_{\rm m}$ in the new code. The relations that the self-supervised objective learns between the dark-matter and baryonic properties of subhalos appear to be partly shared among galaxy-formation models.

Each feedback amplitude rescales a process of one galaxy-formation model relative to the fiducial value of that model, so the same name can stand for different physics in different codes (Table 1). Readouts of these amplitudes use labels from the new code, and a small labeled set suffices (Figure 7). The encoder itself stays unchanged, so adding a code requires a small labeled set and no retraining.

Parameters that scale similar processes in different codes also move the catalogs in similar directions (§ 3.5), and most of these alignments match the responses of galaxies, black holes and gas in the 1P sets ([70] (Delgado et al. 2023; Gebhardt et al. 2024)). They are candidates for feedback coordinates that are defined across codes. An encoder with one part shared by the codes and another specific to each code would make such coordinates explicit and could carry more of the readout calibration across codes. Mercader-Perez et al. (2026) learn such a split from overlapping observations of the same galaxies, separating the physical content of the galaxies from the effects of each survey. In simulations, the code-specific part would describe the galaxy-formation model of each code and would itself be a subject of study. Simulations that share initial conditions and cosmology and differ only in their galaxy-formation code would provide the training pairs, and the 1P sets of IllustrisTNG, SIMBA and Astrid already share their initial random seed ([70] (Villescusa-Navarro et al. 2021)). Once calibrated against quantities that all codes predict, such as the baryon fractions of halos or the abundance of massive black holes, the shared feedback directions would become targets for inference across codes.

Sources of parameter information

Masked reconstruction and the prediction of the teacher’s object representations both ask the encoder to infer the hidden properties of a subhalo from its visible properties and from the objects around it (§ 2.4). Relations such as those between the dark-matter and stellar masses of a subhalo between a satellite and its host change with cosmology and with the galaxy-formation physics of the simulation. The visible objects of a view carry the imprint of these parameters, so the objective rewards an encoder that uses this imprint to predict the hidden properties. The covariance of the object representations collects these relations over a whole catalog and retains the parameters, which the encoder never received (§ 3.2).

Much of the information on $\Omega_{\mathrm{m}}$ comes from the internal properties of the subhalos (§ 3.3). The maximum circular velocity traces the dark matter of a galaxy, and the objective asks the encoder to predict the dark-matter properties of a subhalo from its baryonic ones and the reverse. CAMELS fixes $\Omega_{\mathrm{b}}$, so $\Omega_{\mathrm{m}}$ sets the dark-matter content of galaxies relative to their baryons (Villaescusa-Navarro et al. 2022c), and the baryon fractions of halos carry $\Omega_{\mathrm{m}}$ directly. Part of the precision on $\Omega_{\mathrm{m}}$ therefore reflects the design of the simulations.

With positions alone, the readout stays below the GNNs of Villanueva-Domingo & Villaescusa-Navarro (2022). This input lies outside the training conditions of the encoder. Its training views hide random groups of properties and random objects, and most objects keep some of their properties, so the encoder rarely sees a catalog in which every object is reduced to its position. The GNNs are built on graphs of paired relations and learn the correlations between positions directly. Baldwin et al. (2026) find that GNNs encoding galaxy–galaxy relations make the best use of galaxy positions, and that Deep Sets draw the information on $\Omega_{\mathrm{m}}$ mainly from the line-of-sight velocities. Training views in which every object keeps only its position, or a reconstruction target for the spatial distribution of the objects, would ask the encoder for this clustering information as well.

The kinetic-mode amplitudes of IllustrisTNG and Astrid are the hardest parameters to measure. They act directly on the few subhalos whose black holes are massive enough to drive kinetic feedback, and the confidence of their readouts grows with the number of such hosts (§ 3.6). The covariance of equation (13) weights every subhalo equally, and mass-weighted pooling or attention pooling over the most massive halos would give these hosts more weight. At fixed parameters, the number of hosts grows in proportion to the volume, and a box with four times the CAMELS volume would contain about as many hosts as the Astrid catalogs with the highest confidence. Simulations of this size would measure how the information on kinetic feedback grows with the number of hosts. The label-efficiency curves of these parameters still rise at the largest training set (Figure 7), so more labeled simulations would raise their precision as well.

Uses of the embedding

In astronomy, foundation models such as AstroCLIP (Parker et al. 2024) and AION-1 (Parker et al. 2025) learn embeddings of galaxy images and spectra from which simple heads recover physical properties, redshifts and classifications. Celia plays the same role for galaxy catalogs. It is trained once without labels, and each downstream task fits only a readout on the frozen embedding of the whole catalog.

Before any readout is fitted, the geometry of the embedding shows which parameters move the catalogs of a suite beyond the scatter from initial conditions, and hence which parameters the suite can constrain at its volume. In CAMELS, the kinetic-mode amplitudes of IllustrisTNG and Astrid move a catalog least, and the direction of $A_{\mathrm{AGN1}}$ in IllustrisTNG overlaps with those of the cosmological parameters and of the initial conditions (§ 3.5). Once a readout and its error model are fitted, the confidence score extends this diagnosis to single catalogs. Its mean follows the readout $R^2$ because the error model is calibrated on the training catalogs, and its variation from catalog to catalog picks out the catalogs whose AGN parameters are well measured (§ 3.6).

The principal axes of the embedding and the retrieval of the nearest LH catalogs let one explore a simulation suite without choosing a summary statistic in advance. The axes are found from the catalogs alone, and the statistics and the parameters enter only afterwards, so the same axes can be read with any quantity that the catalogs provide. In Astrid, this exploration finds an axis led by the thermal-mode amplitude along which massive black holes and massive galaxies change in opposite directions (§ 3.7). This matches the response of the Astrid 1P set (Ni et al. [70]) and the self-regulated growth of black holes (Booth & Schaye 2009). Axes that follow the cosmological parameters recur in SIMBA and Astrid, and such recurring axes could guide the choice of parameters for new simulation suites. Retrieval also turns any direction in the embedding into a sequence of simulated catalogs (Figure 9).

The embedding can also enter inference and emulation directly. Its readout predictions are compressed summaries of the kind used to make likelihood-free inference tractable (Alsing et al. 2018), and a neural posterior estimator trained on the embedding would give the joint posterior of several parameters (Cranmer et al. 2020) and may recover more of the AGN information than a linear map. The 1P sets trace the response of the embedding to each parameter with the random seed of the initial conditions held fixed, and the CV sets measure the scatter that the initial conditions add (§ 3.5). A conditional density model of $z$ given the parameters can be trained to reproduce both. Such a model would emulate the embeddings of new simulations without running them, and the linear readouts would apply to its samples unchanged.

Observations and broader simulation suites

The catalogs used here are simulated at $z = 0$, with real-space positions, full three-dimensional velocities and properties that are not directly observable, such as dark-matter masses and $V_{\max}$. Applying the encoder to galaxy surveys requires inputs restricted to observables in redshift space, with survey selection and measurement noise. The frozen encoder already embeds restricted inputs without retraining, and with line-of-sight velocities or stellar properties its readouts reach a precision comparable to that of supervised networks (§ 3.3). Training views with survey selection and redshift-space distortions would extend the encoder to such data, and catalogs at higher redshift can enter the training as additional views.

The encoder was trained on the six-parameter Latin hypercubes of three codes. The SB28 set of IllustrisTNG varies 28 parameters at once, and the 1P-28 sets of IllustrisTNG and SIMBA vary the additional parameters one at a time (Ni et al. [70]). These sets would test whether the shared directions persist as more physical processes vary. Additional codes in the pretraining would show whether the zero-shot transfer improves with the number of codes. A generative decoder from the embedding to catalogs or their statistics would extend the retrieval of § 3.5 and § 3.7 beyond the simulated catalogs.

Conclusions

We have presented Celia, a self-supervised encoder trained without labels on the subhalo catalogs of the IllustrisTNG, SIMBA and Astrid suites of CAMELS and their N-body counterparts. During training, the encoder fills in subhalo properties hidden from a masked view and predicts the representations that an EMA teacher assigns to the same subhalos in the complete view. The frozen encoder then turns each catalog into one vector. Linear readouts of this vector estimate the simulation parameters with error bars. The same vectors define a response direction for each parameter, and the nearest simulated catalogs translate directions in the embedding into catalog statistics. Our main findings are as follows.

1. The embedding retains the cosmological and supernova parameters almost completely. In all three suites, linear readouts reach $R^2 = 0.999$ for $\Omega_{\mathrm{m}}$, 0.986 to 0.991 for $\sigma_8$ and at least 0.97 for both supernova amplitudes. The precision for the AGN amplitudes depends on how many objects they act on. In SIMBA they reach $R^2 = 0.94$ and 0.84, the thermal mode of Astrid reaches 0.996, and the kinetic-mode amplitudes of IllustrisTNG and Astrid reach 0.62 to 0.76. With the restricted inputs of earlier supervised studies, the readouts are comparable to their graph neural networks once the positions are complemented by line-of-sight velocities or stellar properties.

2. Celia transfers to Swift-EAGLE without retraining. The readout fitted on the three training suites follows the Swift-EAGLE values of $\Omega_{\mathrm{m}}$ with $R^2 = 0.825$ and $r = 0.985$, and those of $\sigma_8$ with $r = 0.79$. The feedback amplitudes refer to the model of each code and need labels from the new code. Readouts fitted on 20 labeled Swift-EAGLE catalogs give $r \ge0.94$ for $\Omega_{\mathrm{m}}$, $\sigma_8$ and $A_{\mathrm{SN1}}$, and readouts fitted on 100 give $r \ge0.96$ for every parameter except $A_{\mathrm{AGN2}}$ ($r = 0.85$), similar to the number of labels that the training codes need.

3. The geometry of the embedding separates the effect of each parameter from the scatter of the initial conditions. The cosmological and supernova parameters move a catalog by 4 to 59 times this scatter, and the kinetic-mode amplitudes move it least. Feedback parameters that scale similar processes in different codes point in similar directions, and most of these alignments match the responses of galaxies, black holes and gas in the 1P simulations (Delgado et al. 2023; Ni et al. 2023; Gebhardt et al. 2024). Simulations retrieved along straight lines in the embedding show the changes in their statistics that each parameter is expected to cause.

4. The confidence score of each prediction comes from the embedding alone, and for the AGN parameters it closely follows the massive black-hole content of the catalog. For the kinetic- and jet-mode amplitudes, it rises with the number of black holes massive enough to drive that feedback. For the kinetic mode of Astrid, its mean grows from 0.03 in catalogs without such black holes to 0.84 in catalogs with about 75.

5. The principal axes of the embedding need no labels. Each follows a combination of parameters, and moving along it changes the numbers of massive galaxies and black holes. In Astrid, the axis led by the thermal-mode amplitude raises the number of massive galaxies by a factor of 1.7 and lowers that of massive black holes by a factor of 2.8, matching the response of the Astrid 1P set to stronger thermal feedback [70]. In SIMBA and Astrid, the axes that follow $\Omega_{\mathrm{m}}$ and $\sigma_8$ raise the numbers of both populations.

Celia provides a shared coordinate system for simulated universes. Catalogs from different galaxy-formation models can be compared in it, and its frozen encoder places the catalogs of a new code there without retraining. Inference, diagnostics and the exploration of a simulation suite all operate on one vector per catalog. Larger simulation volumes, survey-like training views and more galaxy-formation models would extend this coordinate system toward observed galaxy catalogs.

Acknowledgments

This work was supported by the Keystone Astronomy and AI (KAAI) program, funded by the Simons Foundation through its Targeted Grants to Institutes program (SFI-MPS-T-Institutes-00021638). We are grateful for discussions with Francisco Villaescusa-Navarro and with participants of the KAAI program, including Yixi Zhao, Sean Lewis, Fatemeh Hafianzadeh and Patrick LaChance. We also acknowledge funding from the Block Center for Technology and Society at Carnegie Mellon University. TDM and RACC acknowledge funding from NASA ATP 80NSSC18K101 and NASA ATP NNX17AK56G. TDM additionally acknowledges funding from NASA ATP 80NSSC20K0519, NSF PHY-2020295, and NASA Theory grant 80NSSC22K072.

The computational results presented have been achieved using Vera at the Pittsburgh Supercomputing Center (PSC). ChatGPT (OpenAI) and Claude (Anthropic) were used to help format certain captions and figures and to polish the language; all text generated by these tools was carefully reviewed, edited, and verified for accuracy by the authors, who assume full responsibility for the outcome.

Input fields and training settings

Input fields. — Tables 3 and 4 list the 53 input features of a subhalo and their transformations before standardization. Masses, radii, rates and characteristic velocities such as $V_{\max}$ enter as base-10 logarithms in catalog units. These units are $10^{10}\ h^{-1}\ \mathrm{M}_{\odot}$ for masses, $h^{-1}\ \mathrm{kpc}$ for lengths, $\mathrm{km}\ \mathrm{s}^{-1}$ for velocities and $\mathrm{M}_{\odot}\ \mathrm{yr}^{-1}$ for star-formation rates. Positions are fractional box coordinates $x=(r\ \bmod\ L)/L$, and all separations use the minimum image in the periodic box. Positions and velocities relative to the host are normalized by the $R_{200\mathrm{c}}$ and $V_{200\mathrm{c}}$ of the host, and the center offsets of the subhalo and of its host by the total half-mass radius and by $R_{200\mathrm{c}}$, respectively. The host gas, stellar and black-hole masses enter as fractions of the host FoF mass, and the host dark-matter fraction as one minus the summed mass fractions of all other particle types.

All features other than the positions are standardized with means and standard deviations computed on the training catalogs. The three components of the velocity and of the host-relative position and velocity share one scale and have zero mean. The Swift-EAGLE catalogs are standardized with the same statistics. The positivity indicators enter the encoder standardized in the same way and serve as reconstruction targets with their original values 0 and 1.

The inputs exclude the particle counts of subhalos and groups, because together with the masses they give the particle mass and hence $\Omega_{\mathrm{m}}$. They also exclude the number of subhalos in a group and the catalog radii $R_{200\mathrm{c}}$ and $R_{500\mathrm{c}}$, which follow from the unperturbed masses. The subhalo total mass is not an input either (§ 2.4.1). The spin enters only through its magnitude, since its direction is statistically isotropic.

Mass perturbation. — The perturbation of equation (2) shifts the dark-matter mass of each subhalo by an offset $\delta_M$ and the three host masses $M_{\mathrm{FoF}}$, $M_{200\mathrm{c}}$ and $M_{500\mathrm{c}}$ by a common offset $\delta_{M,\mathrm{h}}$. The host scales

are computed from the perturbed masses, and so are the host fractions and the quantities normalized by these scales. No combination of inputs therefore returns an unperturbed mass.

Training settings. — We optimize the student with AdamW with $\beta_1=0.9$, $\beta_2=0.95$ and weight decay 0.05 and clip the gradient norm at 10. The weight decay does not act on one-dimensional parameters such as biases and normalization weights. The forward pass uses bfloat16 mixed precision, and the reconstruction loss is computed in single precision. An epoch has 63 batches. The learning rate rises linearly to $2\times10^{-4}$ over the first 200 steps, stays constant until step 10,000 and then decays linearly to zero at step 18,000. The EMA momentum is $m=0.996$ during the warmup and the stable phase and rises to 1 during the decay,

for $10^4\leq t\leq1.8\times10^4$.

Encoder and masking details

This appendix specifies the input groups and masking families (Appendix § 8.§ 1), the spatial modules (Appendix § 8.§ 2), the propagation of masks between related fields (Appendix § 8.§ 3) and the reconstruction targets (Appendix § 8.§ 4).

Input groups and masking families

The last columns of Table 3 and Table 4 assign every field other than the positions to one of the 21 input groups of equation (3). Several groups join a subhalo field with the corresponding host field. The stellar, gas, star-formation and black-hole groups also contain the host stellar, gas and black-hole mass fractions and the host star-formation rate, so step 1 of the masking hides these host fields as well.

In terms of these groups, the four baryonic families are {stellar, stellar radius, stellar metallicity}, {gas, gas metallicity}, {star formation, SF-gas metallicity} and {black hole, BH accretion}. The seven non-baryonic families are {velocity, host-rel. velocity}, {FoF mass, $M_{200\mathrm{c}}$, $M_{500\mathrm{c}}$}, {DM mass}, {half-mass radii}, {dynamics}, {center offset} and {host center offset}. The host-relative group holds the host-relative positions and the central flag. This group and the host DM fraction belong to no family.

Spatial modules

The pair-attention layer of equation (6) has eight heads of dimension 32. One MLP takes the three separation measures as input and returns the eight biases $B_h$. Its output layer is initialized to zero, so that the layer starts as unbiased attention. Each spatial module adds its output to the tokens and is followed by a residual feed-forward block with pre-layer normalization, a hidden width of 1,024 and GELU activation.

Mask propagation

Choosing whole hosts in step 3 hides every subhalo with probability 0.2, whatever the size of its host. Step 4 hides fields that could be computed from visible ones. The host dark-matter fraction equals one minus the summed fractions of the other particle types (Table 4), so it is hidden whenever the host gas, stellar or black-hole fraction is hidden.

The host-relative positions $\Delta\vec r=\Delta\vec r/R'_{200\mathrm{c}}$ and velocities $\Delta\vec v=\Delta\vec v/V'_{200\mathrm{c}}$ carry the scales of the host. The positions are always visible, and for two members $a$ and $b$ of the same host they give

and the velocities give $V'_{200\mathrm{c}}$ in the same way whenever the host-relative velocities are visible. Either scale returns $M'_{200\mathrm{c}}$ through equation (21). Whenever the host masses are hidden, the host-relative positions and velocities, the host center offset and the central flag are therefore hidden for all members of the host. This rule acts on single groups, so it can hide the host-relative velocities and leave the velocities of the same family visible.

Reconstruction targets

Binary targets that a visible indicator already determines are not scored. These are a positive gas mass when the star-formation rate is visible and positive, a positive stellar mass or host stellar fraction when the gas metallicity or the star-formation rate within $2R_{1/2,\star}$ is visible and positive, and a positive host gas fraction when the star-formation rate of the subhalo or of its host is visible and positive.

Readout details

With $d_{\mathcal{T}}$ retained principal components of variances $\lambda_1,\ldots,\lambda_{d_{\mathcal{T}}}$ over $\mathcal{T}$, the scalar of equation (14) is

It keeps the relative variances of the components. The projection and the scalar are fitted on $\mathcal{T}$ and applied unchanged to the test catalogs.

The regularization strength $\alpha$ is chosen separately for each target from 33 logarithmically spaced values between $10^{-3}$ and $10^5$. The training catalogs are split by their LH index into three folds. Every suite uses the same assignment of LH indices to folds, so catalogs with the same LH index fall into the same fold of a pooled readout. Test catalogs are predicted by the readout fitted on all training catalogs, with $\alpha$ minimizing the cross-validated squared error over the three folds. Each training catalog also receives an out-of-fold prediction $\hat y^{\mathrm{oof}}$ from a readout fitted on the other two folds, with $\alpha$ chosen by a three-fold cross-validation within these catalogs.

The error model is a ridge regression with an unpenalized intercept, and its regularization strength is chosen on the same grid and folds as for the readout. For the multipliers, the error model of each fold computes its floor $\eta$ and the variance $V_0$ on the other two folds.

Calibration of the error bars

For each test catalog and parameter, the probability integral transform (PIT) $u$ is the fraction of training catalogs whose signed normalized out-of-fold residual $r_n/\sigma_n^{\mathrm{ooF}}$ lies below $(y-\hat{y})/\hat{\sigma}(z)$ (Section 2.6). For calibrated error bars, $u$ is uniformly distributed over the test catalogs. Figure 14 shows the histograms of $u$ for the six parameters in the three suites. A $\chi^2$ test of the ten bin counts against the uniform expectation gives $p>0.05$ in 13 of the 18 combinations. The five exceptions are $A_{\mathrm{SN2}}$, $A_{\mathrm{AGN1}}$ and $A_{\mathrm{AGN2}}$ in SIMBA and $\Omega_m$ and $A_{\mathrm{AGN2}}$ in Astrid. The largest departure occurs for $A_{\mathrm{AGN2}}$ in SIMBA, where 37 of the 200 true values lie above the 90th percentile of their predictive distributions, against 20 expected. Over the 18 combinations, the 68% and 90% intervals contain 60% to 78% and 83% to 97% of the true values.

This paper was built using the Open Journal of Astrophysics \LaTeX template. The OJA is a journal which provides fast and easy peer review for new papers in the astro-ph section of the arXiv, making the reviewing process simpler for authors and referees alike. Learn more at http://astro.theoj.org.