Circumbinary disks are found in a variety of astrophysical scenarios, spanning binary star formation to accreting supermassive black hole binaries. The interaction with a circumbinary disk can yield opposite effects on the binary orbit leading to circularization, or exciting the eccentricity, widening the orbit or shrinking it and facilitating mergers. We present a new formalism for the long-term evolution of the disk-binary interaction based on the results of recent suites of hydrodynamic simulations, which resolve the complex geometry of the gas in the vicinity of the binary and fully account for the gravitational and accretion forces. We release a python package, \texttt{spindler}, that implements our model. We show that, unless the mass reservoir feeding the disk is comparable to the mass of the binary, accretion onto the binary depletes the disk mass before inducing a significant change in orbital separation or mass ratio. This finding implies that, in most scenarios, interaction with a circumbinary disk is not an efficient mechanism to shrink the orbit of the binary. However, as long as the mass of the disk is at least a few percent of the mass of the binary, the interaction can excite the eccentricity up to an equilibrium value, and induce a statistical correlation between mass ratio and eccentricity. We consider the applicability of our model to a variety of astrophysical scenarios: during star formation, in evolved stellar binaries, triples and in supermassive black hole binaries. We discuss the theoretical and observational implications of our predictions.

Evolved cool stars of various masses are major cosmic engines, delivering substantial mechanical and radiative feedback to the interstellar medium through strong stellar winds and supernova ejecta. These stars play a pivotal role in enriching the interstellar medium with vital chemical elements that constitute the essential building blocks for forming subsequent generations of stars, planets, and potentially even life. Within the complex tapestry of processes occurring in the atmospheres of these cool and luminous stars, convection takes center stage. Convection is a non-local, complex phenomenon marked by non-linear interactions across diverse length scales within a multi-dimensional framework. For these particular stars, characterized by their considerable luminosities and extensive scale heights, convection transitions to a global scale. This transition is facilitated by transmitting radiative energy through the non-uniform outer layers of their atmospheres. To fully understand this phenomenon, the application of global comprehensive 3D radiation-hydrodynamics simulations of stellar convection is of paramount importance. We present two state-of-the-art numerical codes: CO5BOLD and Athena++. Furthermore, we provide a view on their applications as: pivotal roles in enabling a comprehensive investigation into the dynamic processes linked to convection; and critical tools for accurately modeling the emissions produced during shock breakouts in Type II-P Supernovae.

The equipartition model is widely used to estimate magnetic field strength from synchrotron intensity in radio galaxies, yet the validity of its underlying assumptions remains uncertain. Using an Arepo simulation which incorporates a two-moment cosmic ray (CR) transport scheme and a multiphase interstellar medium, we compare magnetic fields inferred from synthetic synchrotron emission maps with the true fields in the simulation. Starting from the derivation of the equipartition formula, we find that the deviation between the equipartition magnetic field and the true magnetic field depends only weakly on the ratio of the magnetic to the CR energy density. In practice, for both face-on and edge-on projections, the equipartition model slightly overestimates the total synchrotron-weighted magnetic field with mean offsets of 32% (0.17 dex) and 36% (0.2 dex), even though the energy equipartition does not hold locally. Beyond these average offsets, a clear trend emerges in edge-on projections that the model underestimates the field in the disk and overestimates it in the halo. Our results demonstrate that the validity of the equipartition model depends only weakly on the strict fulfillment of energy equipartition, and that the equipartition model remains a practical method for estimating magnetic field strengths in face-on projection maps based on our CR-magnetohydrodynamics simulation.

The CGM around unobscured AGN has received much attention in recent years. Comparatively, nebulae associated with obscured AGN are less studied. Here, we simulate the Ly$\alpha$, H$\alpha$, and HeII nebulae around the two types of AGN at $z=2-3$ with ten massive systems from the FIRE simulations based on the unified model to show their differences and to test if they can be used to constrain the AGN model. We post-process the data with the CLOUDY and the Ly$\alpha$ radiative transfer code, RASCAS. Overall, we find that the Ly$\alpha$ nebulae around the unobscured AGN (type-I nebulae) and obscured AGN (type-II nebulae) do not exhibit significant differences in the luminosity, area, and HeII/Ly$\alpha$ when the simulated cutout is set to the halo virial radius. Whereas, the type-II nebulae exhibit less symmetric morphologies, flatter surface brightness profiles, and larger emission line widths (at $R\geq 10$ kpc) than those of the type-I nebulae. These nebulae properties exhibit complicated correlations with the AGN, indicating that nebulae observations can be applied to constrain the AGN engine. However, independent observations on nebulae in the mentioned emissions are insufficient to test the unified model as a priori in observations is not possible to know the direction and opening angle of the ionization cone. We prompt that the joint observations of Ly$\alpha$ nebulae and radio jets can help to reveal the ionization cone to probe the unified model. Our calculations suggest that this method requires $\geq 75$ type-II Ly$\alpha$ nebulae with current instruments to reach a confidence level of $\geq 95\%$.

The polarized light of the cosmic microwave background is sensitive to new physics that violates parity symmetry. For example, the interaction of photons with the fields of elusive dark matter and dark energy could cause a uniform rotation of the plane of linear polarization across the sky, an effect known as cosmic birefringence. We extract the cosmological rotation angle, $\beta$, using Bayesian analysis of parity-violating correlations, $EB$ and $TB$, of polarization data from the Atacama Cosmology Telescope (ACT) Data Release 6. We use prior probabilities for instrumental miscalibration angles derived from the optics model for the ACT telescope and instruments, and marginalize over a residual intensity-to-polarization leakage. We measure $\beta = 0.215^\circ\pm 0.074^\circ$ (68\% confidence level), which excludes $\beta=0$ with a statistical significance of $2.9\sigma$. Although there remain systematics in the ACT data that are not understood and do not allow us to draw strong cosmological conclusions, this result is consistent with previous independent results from the \wmap\ and \planck\ missions. It is suggestive that independent data sets and analyses using different methodologies have yielded the same sign and comparable magnitudes for $\beta$.

We create the first large-scale mock spectroscopic survey of gas absorption sightlines traversing the interstellar medium (ISM), circumgalactic medium (CGM), and intergalactic medium (IGM) surrounding galaxies of virtual Universes. That is, we create mock, or synthetic, absorption spectra by drawing lines-of-sight through cosmological hydrodynamical simulations, using a new mesh-free Voronoi ray-tracing algorithm. The result is the Synthetic Absorption Line Spectral Almanac (SALSA), which is publicly released on a feature-rich online science platform (this http URL). It spans a range of ions, transitions, instruments, observational characteristics, assumptions, redshifts, and simulations. These include, but are not limited to: (ions) HI, OI, CI, MgI, MgII, FeII, SiII, CaII, ZnII, SiIII, SiIV, NV, CII, CIV, OVI; (instruments) SDSS-BOSS, KECK-HIRES, UVES, COS, DESI, 4MOST, WEAVE, XSHOOTER; (model choices) with/without dust depletion, noise, quasar continua, foregrounds; (redshift) from z=0 to z~6; (ancillary data) integrated equivalent widths, column densities, distances and properties of nearby galaxies; (simulations) IllustrisTNG including TNG50, TNG-Cluster, EAGLE, and SIMBA. This scope is not fixed, and will grow and evolve with community interest and requests over time -- suggestions are welcome. The resulting dataset is generic and broadly applicable, enabling diverse science goals such as: (i) studies of the underlying physical gas structures giving rise to particular absorption signatures, (ii) galaxy-absorber and halo-absorber correlations, (iii) virtual surveys and survey strategy optimization, (iv) stacking experiments and the identification of faint absorption features, (v) assessment of data reduction methods and completeness calculations, (vi) inference of physical properties from observables, and (vii) apples-to-apples comparisons between simulations and data.

The Lyman-alpha (Ly$\alpha$) forest at $z \sim 5$ offers a primary probe to constrain the history of the Epoch of Reionization (EoR), retaining thermal and ionization signatures imprinted by the reionization process. In this work, we present a new inference framework based on JAX that combines forward-modeled Ly$\alpha$ forest observables with differentiable neural emulators and Hamiltonian Monte Carlo (HMC). We construct a dataset of 501 low-resolution simulations generated with user-defined reionization histories and compute a set of 1D Ly$\alpha$ power spectra and model-dependent covariance matrices. We then train two independent neural emulators that achieve sub-percent errors across relevant scales and combine them with HMC to efficiently perform parameter estimation. We validate this framework by applying it to a suite of mock observations, demonstrating that the true parameters are reliably recovered. While this work is limited by the low resolution of the simulations used, our results highlight the potential of this method for inferring the reionization history from high-redshift Ly$\alpha$ forest measurements. Future improvements in our reionization models will further enhance its ability to extract constraints from observational datasets.

Weak gravitational lensing is a powerful tool for precision tests of cosmology. As the expected deflection angles are small, predictions based on non-linear N-body simulations are commonly computed with the Born approximation. Here we examine this assumption using ${\small DORIAN}$, a newly developed full-sky ray-tracing scheme applied to high-resolution mass-shell outputs of the two largest simulations in the MillenniumTNG suite, each with a 3000 Mpc box containing almost 1.1 trillion cold dark matter particles in addition to 16.7 billion particles representing massive neutrinos. We examine simple two-point statistics like the angular power spectrum of the convergence field, as well as statistics sensitive to higher order correlations such as peak and minimum statistics, void statistics, and Minkowski functionals of the convergence maps. Overall, we find only small differences between the Born approximation and a full ray-tracing treatment. While these are negligibly small at power-spectrum level, some higher order statistics show more sizable effects; ray-tracing is necessary to achieve percent level precision. At the resolution reached here, full-sky maps with 0.8 billion pixels and an angular resolution of 0.43 arcmin, we find that interpolation accuracy can introduce appreciable errors in ray-tracing results. We therefore implemented an interpolation method based on nonuniform fast Fourier transforms (NUFFT) along with more traditional methods. Bilinear interpolation introduces significant smoothing, while nearest grid point sampling agrees well with NUFFT, at least for our fiducial source redshift, $z_s=1.0$, and for the 1 arcmin smoothing we use for higher-order statistics.

It is often understood that supernova (SN) feedback in galaxies is responsible for regulating star formation and generating gaseous outflows. However, a detailed look at their effect on the local interstellar medium (ISM) on small mass scales in simulations shows that these processes proceed in clearly distinct channels. We demonstrate this finding in two independent simulations with solar-mass resolution, LYRA and RIGEL, of an isolated dwarf galaxy. Focusing on the immediate environment surrounding SNe, our findings suggest that the large-scale effect of a given SN on the galaxy is best predicted by its immediate local density. Outflows are driven by SNe in diffuse regions expanding to their cooling radii on large (\sim kpc) scales, while dense star-forming regions are disrupted in a localized (\sim pc) manner. However, these separate feedback channels are only distinguishable at very high numerical resolutions capable of following scales <<10^3 \msun. On larger scales, ISM densities are greatly mis-estimated, and differences between local environments of SNe become severely washed out. We demonstrate the practical implications of this effect by comparing with a mid-resolution simulation (Mptcl. \sim 200 \msun) of the same dwarf using the SMUGGLE model. The coarse-resolution simulation cannot self-consistently determine whether a given SN is responsible for generating outflows or suppressing star formation, suggesting that emergent galaxy physics such as star formation regulation through hot-phase outflows is fundamentally unresolvable by subgrid stellar feedback models, without appealing directly to simulations with highly resolved ISM.

Recent Pulsar Timing Arrays (PTAs) results provided strong evidence for a stochastic gravitational wave background (sGWB), consistent with a population of merging massive black holes (MBHs) at z&lt;1. Meanwhile, JWST observations at z&gt;5 suggest a higher number density of accreting MBHs than previously estimated. Together with constraints from local MBHs and high-$z$ quasars, these findings offer a unique opportunity to test MBH seeding and early growth models. We explore this using ${\tt L-Galaxies}\textit{BH}$, a new extension of the galaxy formation model ${\tt L-Galaxies}$, developed to explicitly model all stages of MBH evolution, including seeding, accretion, and binary dynamics. To take advantage of both the high resolution of the ${\tt MillenniumII}$ and the large volume of the ${\tt Millennium}$ simulations, we run ${\tt L-Galaxies}\textit{BH}$ on the former and use its outputs as initial conditions for the latter, via our $\textit{grafting}$ method. We find that reproducing the number density of high-$z$ active MBHs observed by JWST requires either a heavy seed formation rate significantly higher than that predicted by current models ($\gtrsim 0.01 Mpc^{-3}$ at $z \sim 10$), or widespread formation of light seeds undergoing multiple phases of super-Eddington accretion. Furthermore, matching the amplitude of the PTA sGWB signal requires nearly all galaxies with stellar masses M_{*}&gt; 10^9 M_\odot to host central MBHs by $z\sim0$. Given the extreme heavy seed densities required to satisfy both PTA and JWST constraints, our results favor a scenario in which MBHs originate from light seeds that grow rapidly and efficiently in the early universe. This work demonstrates the power of combining multi-messenger data with physical models to probe the origins and evolution of MBHs across cosmic time.

We present James Webb Space Telescope (JWST) NIRSpec 1.7--5.5 micron observations of SN~2024ggi at +285.51 and +385.27 days post-explosion. The late-time nebular spectra are dominated by emission lines from various ionization states of H, Ca, Ar, C, Mg, Ni, Co, and Fe. We also detect strong CO emission in both the first overtone and fundamental vibrational bands. Most atomic features exhibit asymmetric line profiles, indicating an aspherical explosion. Using observed fluxes combined with non-LTE radiative-transfer simulations, we develop a data-driven method that resolves the complex molecular-emission region, constrains its 3D structure, and reproduces high-fidelity spectral profiles. We find that, CO is mostly formed prior to +285d past explosion. The subsequent evolution is dominated by the evaporation of CO with CO mass varying from M(CO) of 8.7E-3 to 1.3E-3 Mo, and with instabilities growing from almost homogeneous to highly clumped (density contrast f_c of 1.2 to 2). The minimum velocity of CO only slightly decreases between epochs (v_1 of 1200 and 1100 km/sec), with the reference temperature dropping from T_1 of 2400 and 1900K.

Luminous Red Novae (LRNe) are transients associated with mass ejection during stellar mergers and common envelope evolution (CEE). LRNe have the potential to illuminate the poorly understood phases of binary evolution leading up to the CEE, during the mass ejection phase, and in the immediate aftermath. However, the mechanism responsible for powering LRN light curves and the origin of their observed diversity remain open questions. Here, we perform two-dimensional moving-mesh radiation-hydrodynamic simulations of LRNe that take into account hydrogen and helium recombination and relevant opacities. We study a typical high-mass stellar merger, which dynamically ejects 2 $M_\odot$ with a characteristic velocity of 410 km/s. This ejecta collides with 2.7 $M_\odot$ of equatorially concentrated circumbinary material (CBM) left behind from a prior phase of non-conservative runaway mass transfer. We find that the resulting light curve is composed of a short, blue peak followed by a redder, predominantly shock-powered plateau with luminosities reaching up to $10^{41}$ erg/s and durations up to 200 days. These luminosities are significantly higher, and the durations much longer, than those produced by a simple spherical ejection of the same mass. They also depend in a complex way on the radial distribution of the CBM and the viewing angle. The shock is embedded in the ejecta and its observational signatures during the optically-thick phase are largely hidden. Our results are broadly compatible with observations of the brightest extragalactic LRNe and pave the way for the transformation of LRNe into powerful probes of binary evolution.

Little Red Dots (LRDs) are compact sources at $z>5$ discovered through JWST spectroscopy. Their spectra exhibit broad Balmer emission lines ($\gtrsim1000\rm~km~s^{-1}$), alongside absorption features and a pronounced Balmer break -- evidence for a dense, neutral hydrogen medium with the $n=2$ state. When interpreted as arising from AGN broad-line regions, inferred black hole masses from local scaling relations exceed expectations given their stellar masses, challenging models of early black hole--galaxy co-evolution. However, radiative transfer effects in dense media may also impact the formation of hydrogen emission lines. We model three scattering processes shaping hydrogen line profiles: resonance scattering by hydrogen in the $n=2$ state, Raman scattering of UV radiation by ground-state hydrogen, and Thomson scattering by free electrons. Using 3D Monte Carlo radiative transfer simulations with multi-branching resonance transitions, we examine their imprint on line shapes and ratios. Resonance scattering produces strong deviations from Case B flux ratios, clear differences between H$\alpha$ and H$\beta$, and encodes gas kinematics in line profiles but cannot broaden H$\beta$ due to conversion to Pa$\alpha$. While Raman scattering can yield broad wings, scattering of UV continuum is disfavored given the absence of strong FWHM variations across transitions. Raman scattering of higher Lyman-series emission can produce H$\alpha$/H$\beta$ wing width ratios of $\gtrsim1.28$, agreeing with observations. Thomson scattering can reproduce the observed $\gtrsim1000~\rm km\, s^{-1}$ wings under plausible conditions, e.g., $T_{\rm e} \sim 10^4\rm \, K$ and $N_{\rm e}\sim10^{24}\rm~cm^{-2}$ -- and lead to black hole mass overestimates by factors $\gtrsim10$. Our results provide a framework for interpreting hydrogen lines in LRDs and similar systems.

Small-scale winds driven from accretion discs surrounding active galactic nuclei (AGN) are expected to launch kpc-scale outflows into their host galaxies. However, the ways in which the structure of the interstellar medium (ISM) affects the multiphase content and impact of the outflow remains uncertain. We present a series of numerical experiments featuring a realistic small-scale AGN wind with velocity $5\times 10^3-10^4\ \rm{km/s}$ interacting with an isolated galaxy disc with a manually-controlled clumpy ISM, followed at sub-pc resolution. Our simulations are performed with AREPO and probe a wide range of AGN luminosities ($L=10^{43-47}\ \rm{erg/s}$) and ISM substructures. In homogeneous discs, the AGN wind sweeps up an outflowing, cooling shell, where the emerging cold phase dominates the mass and kinetic energy budgets, reaching a momentum flux $\dot{p} \approx 7\ L/c$. However, when the ISM is clumpy, outflow properties are profoundly different. They contain small, long-lived (&gt; 5\ \rm{Myr}), cold (T&lt;10^{4.5}\ \rm{K}) cloudlets entrained in the faster, hot outflow phase, which are only present in the outflow if radiative cooling is included in the simulation. While the cold phase dominates the mass of the outflow, most of the kinetic luminosity is now carried by a tenuous, hot phase with T &gt; 10^7 \ \rm K. While the hot phases reaches momentum fluxes $\dot{p} \approx (1 - 5)\ L/c$, energy-driven bubbles couple to the cold phase inefficiently, producing modest momentum fluxes \dot{p} &lt; L/c in the fast-outflowing cold gas. These low momentum fluxes could lead to the outflows being misclassified as momentum-driven using common observational diagnostics. We also show predictions for scaling relations between outflow properties and AGN luminosity and discuss the challenges in constraining outflow driving mechanisms and kinetic coupling efficiencies using observed quantities.

The 21-cm line of hydrogen is the most promising probe of the Dark Ages and Cosmic Dawn. We combine hydrodynamical simulations with a large-scale grid in order to calculate the effect of non-linear structure formation on the large-scale 21-cm power spectrum, focusing on redshifts $z=20-40$. As the clumping effect arises from small-scale density fluctuations, it offers a unique opportunity to probe the standard cold dark matter model in a new regime and thus potentially investigate the properties of dark matter. To this end, we also study a warm dark matter $-$ like model with a Gaussian cutoff on a scale of 50 kpc. We find that clumping has a significant impact on the large-scale 21-cm power spectrum. For example, for the Dark Ages case at $z=30$ and wavenumber $k=0.05$ Mpc$^{-1}$, small-scale clustering enhances the 21-cm power spectrum by $13\%$. Once Lyman-$\alpha$ coupling kicks in due to the first stars, the 21-cm signal strengthens, and the effect of clumping grows; it suppresses the observable power spectrum at $z=20$ by a factor of two, while the cutoff model has less than half the clumping impact. The clumping effect is significantly higher than the sensitivity of the planned Square Kilometre Array (SKA) AA$^\star$ configuration, by up to a factor of 20 for standard cold dark matter, though detection will require separation from foregrounds and from astrophysical contributions to the 21-cm power spectrum.

Forward models of the galaxy density field enable simulation based inference as well as field level inference of galaxy clustering. However, these analysis techniques require forward models that are both computationally fast and robust to modeling uncertainties in the relation between galaxies and matter. Both requirements can be addressed with the Effective Field Theory of Large Scale Structure. Here, we focus on the physical and numerical convergence of the LEFTfield model. Based on the perturbative nature of the forward model, we derive an analytic understanding of the leading numerical errors, and we compare our estimates to high-resolution and N-body references. This allows us to derive a set of best-practice recommendations for the numerical accuracy parameters, which are completely specified by the desired order of the perturbative solution and the cut-off scale. We verify these recommendations by an extended set of parameter recovery tests from fully nonlinear mock data and find very consistent results. A single evaluation of the forward model takes seconds, making cosmological analyses of galaxy clustering data based on forward models computationally feasible.

We present measurements of black hole masses and Eddington ratios for a sample of 38 bright (M$_{1450}$ < -24.4 mag) quasars at 5.8 < z < 7.5, derived from VLT/X-shooter near-IR spectroscopy of their broad CIV and MgII emission lines. The black hole masses (on average M$_{BH}$ ~ 4.6 x 10$^9$ M$_{\odot}$) and accretion rates (with Eddington ratios ranging between 0.1 and 1.0) are broadly consistent with that of similarly luminous 0.3 < z < 2.3 quasars, but there is evidence for a mild increase in the median Eddington ratio going towards z > 6. Combined with deep ALMA observations of the [CII] 158 $\mu$m line from the quasar host galaxies and VLT/MUSE investigations of the extended Ly$\alpha$ halos, this study provides fundamental clues to models of the formation and growth of the first massive galaxies and black holes. Compared to local scaling relations, z > 5.7 black holes appear to be over-massive with respect to their host galaxies, and their accretion properties do not change with host galaxy morphology. Under the assumption that the kinematics of the T ~ 10$^4$ K gas, traced by the extended Ly$\alpha$ halos, are dominated by the gravitational potential of the dark matter halo, we report a similar relation between the black hole mass and circular velocity to the one reported for z ~ 0 galaxies. These results paint a picture where the first supermassive black holes reside in massive halos at z > 6 and lead the first stages of galaxy formation by rapidly growing in mass with a duty cycle of order unity. However, this duty cycle needs to drastically drop towards lower redshifts, while the host galaxies continue forming stars at a rate of hundreds of solar masses per year, sustained by the large reservoirs of cool gas surrounding them.

This is the second paper in the HOWLS (higher-order weak lensing statistics) series exploring the usage of non-Gaussian statistics for cosmology inference within \textit{Euclid}. With respect to our first paper, we develop a full tomographic analysis based on realistic photometric redshifts which allows us to derive Fisher forecasts in the ($\sigma_8$, $w_0$) plane for a \textit{Euclid}-like data release 1 (DR1) setup. We find that the 5 higher-order statistics (HOSs) that satisfy the Gaussian likelihood assumption of the Fisher formalism (1-point probability distribution function, $\ell$1-norm, peak counts, Minkowski functionals, and Betti numbers) each outperform the shear 2-point correlation functions by a factor $2.5$ on the $w_0$ forecasts, with only marginal improvement when used in combination with 2-point estimators, suggesting that every HOS is able to retrieve both the non-Gaussian and Gaussian information of the matter density field. The similar performance of the different estimators\inlinecomment{, with a slight preference for Minkowski functionals and 1-point probability distribution function,} is explained by a homogeneous use of multi-scale and tomographic information, optimized to lower computational costs. These results hold for the $3$ mass mapping techniques of the \textit{Euclid} pipeline: aperture mass, Kaiser--Squires, and Kaiser--Squires plus, and are unaffected by the application of realistic star masks. Finally, we explore the use of HOSs with the Bernardeau--Nishimichi--Taruya (BNT) nulling scheme approach, finding promising results towards applying physical scale cuts to HOSs.

The abundance of dark matter haloes is a key cosmological probe in forthcoming galaxy surveys. The theoretical understanding of the halo mass function (HMF) is limited by our incomplete knowledge of the origin of non-universality and its cosmological parameter dependence. We present a deep learning model which compresses the linear matter power spectrum into three independent factors which are necessary and sufficient to describe the $z=0$ HMF from the state-of-the-art AEMULUS emulator to sub-per cent accuracy in a $w$CDM$+N_\mathrm{eff}$ parameter space. Additional information about growth history does not improve the accuracy of HMF predictions if the matter power spectrum is already provided as input, because required aspects of the former can be inferred from the latter. The three factors carry information about the universal and non-universal aspects of the HMF, which we interrogate via the information-theoretic measure of mutual information. We find that non-universality is captured by recent growth history after matter-dark-energy equality and $N_\mathrm{eff}$ for $M\sim 10^{13} \, \mathrm{M_\odot}\, h^{-1}$ haloes, and by $\Omega_{\rm m}$ for $M\sim 10^{15} \, \mathrm{M_\odot}\, h^{-1}$. The compact representation learnt by our model can inform the design of emulator training sets to achieve high emulator accuracy with fewer simulations.